Numeric.SpecFunctions:invIncompleteGamma from math-functions-0.1.5.2, B

Average Error: 28.83% → 0.34%

Time: 45.8s

Precision: binary64

Math

\[1 - \log \left(1 - \frac{x - y}{1 - y}\right) \]

↓

\[\begin{array}{l} t_0 := \mathsf{log1p}\left(\frac{x - y}{y - 1}\right)\\ t_1 := \frac{\frac{1 - x}{x - 1}}{y}\\ t_2 := \log \left(\frac{1}{y}\right) + \log \left(x - 1\right)\\ t_3 := 1 - t_0\\ \mathbf{if}\;y \leq -410000000000:\\ \;\;\;\;\begin{array}{l} \mathbf{if}\;1 - x > 0 \land \frac{-1}{y} > 0:\\ \;\;\;\;1 + \left(t_1 - \log \left(\frac{x - 1}{y}\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\left(t_1 - -1\right) - \left(\log \left(\frac{-1}{y}\right) + \log \left(1 - x\right)\right)\\ \end{array}\\ \mathbf{elif}\;y \leq 1.9 \cdot 10^{+44}:\\ \;\;\;\;t_3\\ \mathbf{elif}\;t_0 + 1 \ne 0:\\ \;\;\;\;\frac{1 - {t_2}^{2}}{1 + t_2}\\ \mathbf{else}:\\ \;\;\;\;t_3\\ \end{array} \]

FPCore

(FPCore (x y) :precision binary64 (- 1.0 (log (- 1.0 (/ (- x y) (- 1.0 y))))))

↓

(FPCore (x y)
 :precision binary64
 (let* ((t_0 (log1p (/ (- x y) (- y 1.0))))
        (t_1 (/ (/ (- 1.0 x) (- x 1.0)) y))
        (t_2 (+ (log (/ 1.0 y)) (log (- x 1.0))))
        (t_3 (- 1.0 t_0)))
   (if (<= y -410000000000.0)
     (if (and (> (- 1.0 x) 0.0) (> (/ -1.0 y) 0.0))
       (+ 1.0 (- t_1 (log (/ (- x 1.0) y))))
       (- (- t_1 -1.0) (+ (log (/ -1.0 y)) (log (- 1.0 x)))))
     (if (<= y 1.9e+44)
       t_3
       (if (!= (+ t_0 1.0) 0.0) (/ (- 1.0 (pow t_2 2.0)) (+ 1.0 t_2)) t_3)))))

C

double code(double x, double y) {
	return 1.0 - log((1.0 - ((x - y) / (1.0 - y))));
}

↓

double code(double x, double y) {
	double t_0 = log1p(((x - y) / (y - 1.0)));
	double t_1 = ((1.0 - x) / (x - 1.0)) / y;
	double t_2 = log((1.0 / y)) + log((x - 1.0));
	double t_3 = 1.0 - t_0;
	double tmp_1;
	if (y <= -410000000000.0) {
		double tmp_2;
		if (((1.0 - x) > 0.0) && ((-1.0 / y) > 0.0)) {
			tmp_2 = 1.0 + (t_1 - log(((x - 1.0) / y)));
		} else {
			tmp_2 = (t_1 - -1.0) - (log((-1.0 / y)) + log((1.0 - x)));
		}
		tmp_1 = tmp_2;
	} else if (y <= 1.9e+44) {
		tmp_1 = t_3;
	} else if ((t_0 + 1.0) != 0.0) {
		tmp_1 = (1.0 - pow(t_2, 2.0)) / (1.0 + t_2);
	} else {
		tmp_1 = t_3;
	}
	return tmp_1;
}

Java

public static double code(double x, double y) {
	return 1.0 - Math.log((1.0 - ((x - y) / (1.0 - y))));
}

↓

public static double code(double x, double y) {
	double t_0 = Math.log1p(((x - y) / (y - 1.0)));
	double t_1 = ((1.0 - x) / (x - 1.0)) / y;
	double t_2 = Math.log((1.0 / y)) + Math.log((x - 1.0));
	double t_3 = 1.0 - t_0;
	double tmp_1;
	if (y <= -410000000000.0) {
		double tmp_2;
		if (((1.0 - x) > 0.0) && ((-1.0 / y) > 0.0)) {
			tmp_2 = 1.0 + (t_1 - Math.log(((x - 1.0) / y)));
		} else {
			tmp_2 = (t_1 - -1.0) - (Math.log((-1.0 / y)) + Math.log((1.0 - x)));
		}
		tmp_1 = tmp_2;
	} else if (y <= 1.9e+44) {
		tmp_1 = t_3;
	} else if ((t_0 + 1.0) != 0.0) {
		tmp_1 = (1.0 - Math.pow(t_2, 2.0)) / (1.0 + t_2);
	} else {
		tmp_1 = t_3;
	}
	return tmp_1;
}

Python

def code(x, y):
	return 1.0 - math.log((1.0 - ((x - y) / (1.0 - y))))

↓

def code(x, y):
	t_0 = math.log1p(((x - y) / (y - 1.0)))
	t_1 = ((1.0 - x) / (x - 1.0)) / y
	t_2 = math.log((1.0 / y)) + math.log((x - 1.0))
	t_3 = 1.0 - t_0
	tmp_1 = 0
	if y <= -410000000000.0:
		tmp_2 = 0
		if ((1.0 - x) > 0.0) and ((-1.0 / y) > 0.0):
			tmp_2 = 1.0 + (t_1 - math.log(((x - 1.0) / y)))
		else:
			tmp_2 = (t_1 - -1.0) - (math.log((-1.0 / y)) + math.log((1.0 - x)))
		tmp_1 = tmp_2
	elif y <= 1.9e+44:
		tmp_1 = t_3
	elif (t_0 + 1.0) != 0.0:
		tmp_1 = (1.0 - math.pow(t_2, 2.0)) / (1.0 + t_2)
	else:
		tmp_1 = t_3
	return tmp_1

Julia

function code(x, y)
	return Float64(1.0 - log(Float64(1.0 - Float64(Float64(x - y) / Float64(1.0 - y)))))
end

↓

function code(x, y)
	t_0 = log1p(Float64(Float64(x - y) / Float64(y - 1.0)))
	t_1 = Float64(Float64(Float64(1.0 - x) / Float64(x - 1.0)) / y)
	t_2 = Float64(log(Float64(1.0 / y)) + log(Float64(x - 1.0)))
	t_3 = Float64(1.0 - t_0)
	tmp_1 = 0.0
	if (y <= -410000000000.0)
		tmp_2 = 0.0
		if ((Float64(1.0 - x) > 0.0) && (Float64(-1.0 / y) > 0.0))
			tmp_2 = Float64(1.0 + Float64(t_1 - log(Float64(Float64(x - 1.0) / y))));
		else
			tmp_2 = Float64(Float64(t_1 - -1.0) - Float64(log(Float64(-1.0 / y)) + log(Float64(1.0 - x))));
		end
		tmp_1 = tmp_2;
	elseif (y <= 1.9e+44)
		tmp_1 = t_3;
	elseif (Float64(t_0 + 1.0) != 0.0)
		tmp_1 = Float64(Float64(1.0 - (t_2 ^ 2.0)) / Float64(1.0 + t_2));
	else
		tmp_1 = t_3;
	end
	return tmp_1
end

Wolfram

code[x_, y_] := N[(1.0 - N[Log[N[(1.0 - N[(N[(x - y), $MachinePrecision] / N[(1.0 - y), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]

↓

code[x_, y_] := Block[{t$95$0 = N[Log[1 + N[(N[(x - y), $MachinePrecision] / N[(y - 1.0), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]}, Block[{t$95$1 = N[(N[(N[(1.0 - x), $MachinePrecision] / N[(x - 1.0), $MachinePrecision]), $MachinePrecision] / y), $MachinePrecision]}, Block[{t$95$2 = N[(N[Log[N[(1.0 / y), $MachinePrecision]], $MachinePrecision] + N[Log[N[(x - 1.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]}, Block[{t$95$3 = N[(1.0 - t$95$0), $MachinePrecision]}, If[LessEqual[y, -410000000000.0], If[And[Greater[N[(1.0 - x), $MachinePrecision], 0.0], Greater[N[(-1.0 / y), $MachinePrecision], 0.0]], N[(1.0 + N[(t$95$1 - N[Log[N[(N[(x - 1.0), $MachinePrecision] / y), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(t$95$1 - -1.0), $MachinePrecision] - N[(N[Log[N[(-1.0 / y), $MachinePrecision]], $MachinePrecision] + N[Log[N[(1.0 - x), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], If[LessEqual[y, 1.9e+44], t$95$3, If[Unequal[N[(t$95$0 + 1.0), $MachinePrecision], 0.0], N[(N[(1.0 - N[Power[t$95$2, 2.0], $MachinePrecision]), $MachinePrecision] / N[(1.0 + t$95$2), $MachinePrecision]), $MachinePrecision], t$95$3]]]]]]]

Target

Original	28.83%
Target	0.16%
Herbie	0.34%

\[\begin{array}{l} \mathbf{if}\;y < -81284752.61947241:\\ \;\;\;\;1 - \log \left(\frac{x}{y \cdot y} - \left(\frac{1}{y} - \frac{x}{y}\right)\right)\\ \mathbf{elif}\;y < 3.0094271212461764 \cdot 10^{+25}:\\ \;\;\;\;\log \left(\frac{e^{1}}{1 - \frac{x - y}{1 - y}}\right)\\ \mathbf{else}:\\ \;\;\;\;1 - \log \left(\frac{x}{y \cdot y} - \left(\frac{1}{y} - \frac{x}{y}\right)\right)\\ \end{array} \]

Derivation

1. Split input into 3 regimes

2. if y < -4.1e11

   1. Initial program 83.34

      \[1 - \log \left(1 - \frac{x - y}{1 - y}\right) \]

   2. Simplified 83.34

      \[\leadsto \color{blue}{1 - \mathsf{log1p}\left(\frac{x - y}{y - 1}\right)} \]
      Proof

   3. Taylor expanded in y around -inf 0.5

      \[\leadsto \color{blue}{\left(1 + -1 \cdot \frac{\frac{x}{x - 1} - \frac{1}{x - 1}}{y}\right) - \left(\log \left(-1 \cdot \left(x - 1\right)\right) + \log \left(\frac{-1}{y}\right)\right)} \]

   4. Simplified 0.5

      \[\leadsto \color{blue}{\left(1 + \left(-\frac{\frac{x}{x - 1} - \frac{1}{x - 1}}{y}\right)\right) - \left(\log \left(-\left(x - 1\right)\right) + \log \left(\frac{-1}{y}\right)\right)} \]
      Proof

   5. Applied egg-rr 0.02

      \[\leadsto \color{blue}{\begin{array}{l} \color{blue}{\mathbf{if}\;1 - x > 0 \land \frac{-1}{y} > 0:\\ \;\;\;\;1 + \left(\frac{\frac{1 - x}{x - 1}}{y} - \log \left(\frac{x - 1}{y}\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\left(\frac{\frac{1 - x}{x - 1}}{y} - -1\right) - \left(\log \left(\frac{-1}{y}\right) + \log \left(1 - x\right)\right)\\ } \end{array}} \]

   if -4.1e11 < y < 1.9000000000000001e44

   1. Initial program 0.36

      \[1 - \log \left(1 - \frac{x - y}{1 - y}\right) \]

   2. Simplified 0.33

      \[\leadsto \color{blue}{1 - \mathsf{log1p}\left(\frac{x - y}{y - 1}\right)} \]
      Proof

   if 1.9000000000000001e44 < y

   1. Initial program 53.74

      \[1 - \log \left(1 - \frac{x - y}{1 - y}\right) \]

   2. Simplified 53.74

      \[\leadsto \color{blue}{1 - \mathsf{log1p}\left(\frac{x - y}{y - 1}\right)} \]
      Proof

   3. Applied egg-rr 53.81

      \[\leadsto \color{blue}{\begin{array}{l} \color{blue}{\mathbf{if}\;\mathsf{log1p}\left(\frac{x - y}{y - 1}\right) + 1 \ne 0:\\ \;\;\;\;\frac{1 - {\left(\mathsf{log1p}\left(\frac{x - y}{y - 1}\right)\right)}^{2}}{\mathsf{log1p}\left(\frac{x - y}{y - 1}\right) + 1}\\ \mathbf{else}:\\ \;\;\;\;1 - \mathsf{log1p}\left(\frac{x - y}{y - 1}\right)\\ } \end{array}} \]

   4. Taylor expanded in y around inf 1.54

      \[\leadsto \begin{array}{l} \mathbf{if}\;\mathsf{log1p}\left(\frac{x - y}{y - 1}\right) + 1 \ne 0:\\ \;\;\;\;\color{blue}{\frac{1 - {\left(\log \left(\frac{1}{y}\right) + \log \left(x - 1\right)\right)}^{2}}{1 + \left(\log \left(\frac{1}{y}\right) + \log \left(x - 1\right)\right)}}\\ \mathbf{else}:\\ \;\;\;\;1 - \mathsf{log1p}\left(\frac{x - y}{y - 1}\right)\\ \end{array} \]
3. Recombined 3 regimes into one program.

Reproduce

herbie shell --seed 2023136 
(FPCore (x y)
  :name "Numeric.SpecFunctions:invIncompleteGamma from math-functions-0.1.5.2, B"
  :precision binary64

  :herbie-target
  (if (< y -81284752.61947241) (- 1.0 (log (- (/ x (* y y)) (- (/ 1.0 y) (/ x y))))) (if (< y 3.0094271212461764e+25) (log (/ (exp 1.0) (- 1.0 (/ (- x y) (- 1.0 y))))) (- 1.0 (log (- (/ x (* y y)) (- (/ 1.0 y) (/ x y)))))))

  (- 1.0 (log (- 1.0 (/ (- x y) (- 1.0 y))))))