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

Average Error: 2.1 → 0.8

Time: 23.6s

Precision: binary64

Cost: 33860

Math

\[x \cdot e^{y \cdot \left(\log z - t\right) + a \cdot \left(\log \left(1 - z\right) - b\right)} \]

↓

\[\begin{array}{l} t_1 := y \cdot \left(\log z - t\right) + a \cdot \left(\log \left(1 - z\right) - b\right)\\ \mathbf{if}\;t_1 \leq -1 \cdot 10^{-6}:\\ \;\;\;\;x \cdot e^{t_1}\\ \mathbf{else}:\\ \;\;\;\;x \cdot e^{a \cdot \left(\mathsf{log1p}\left(-z\right) - b\right)}\\ \end{array} \]

FPCore

(FPCore (x y z t a b)
 :precision binary64
 (* x (exp (+ (* y (- (log z) t)) (* a (- (log (- 1.0 z)) b))))))

↓

(FPCore (x y z t a b)
 :precision binary64
 (let* ((t_1 (+ (* y (- (log z) t)) (* a (- (log (- 1.0 z)) b)))))
   (if (<= t_1 -1e-6) (* x (exp t_1)) (* x (exp (* a (- (log1p (- z)) b)))))))

C

double code(double x, double y, double z, double t, double a, double b) {
	return x * exp(((y * (log(z) - t)) + (a * (log((1.0 - z)) - b))));
}

↓

double code(double x, double y, double z, double t, double a, double b) {
	double t_1 = (y * (log(z) - t)) + (a * (log((1.0 - z)) - b));
	double tmp;
	if (t_1 <= -1e-6) {
		tmp = x * exp(t_1);
	} else {
		tmp = x * exp((a * (log1p(-z) - b)));
	}
	return tmp;
}

Java

public static double code(double x, double y, double z, double t, double a, double b) {
	return x * Math.exp(((y * (Math.log(z) - t)) + (a * (Math.log((1.0 - z)) - b))));
}

↓

public static double code(double x, double y, double z, double t, double a, double b) {
	double t_1 = (y * (Math.log(z) - t)) + (a * (Math.log((1.0 - z)) - b));
	double tmp;
	if (t_1 <= -1e-6) {
		tmp = x * Math.exp(t_1);
	} else {
		tmp = x * Math.exp((a * (Math.log1p(-z) - b)));
	}
	return tmp;
}

Python

def code(x, y, z, t, a, b):
	return x * math.exp(((y * (math.log(z) - t)) + (a * (math.log((1.0 - z)) - b))))

↓

def code(x, y, z, t, a, b):
	t_1 = (y * (math.log(z) - t)) + (a * (math.log((1.0 - z)) - b))
	tmp = 0
	if t_1 <= -1e-6:
		tmp = x * math.exp(t_1)
	else:
		tmp = x * math.exp((a * (math.log1p(-z) - b)))
	return tmp

Julia

function code(x, y, z, t, a, b)
	return Float64(x * exp(Float64(Float64(y * Float64(log(z) - t)) + Float64(a * Float64(log(Float64(1.0 - z)) - b)))))
end

↓

function code(x, y, z, t, a, b)
	t_1 = Float64(Float64(y * Float64(log(z) - t)) + Float64(a * Float64(log(Float64(1.0 - z)) - b)))
	tmp = 0.0
	if (t_1 <= -1e-6)
		tmp = Float64(x * exp(t_1));
	else
		tmp = Float64(x * exp(Float64(a * Float64(log1p(Float64(-z)) - b))));
	end
	return tmp
end

Wolfram

code[x_, y_, z_, t_, a_, b_] := N[(x * N[Exp[N[(N[(y * N[(N[Log[z], $MachinePrecision] - t), $MachinePrecision]), $MachinePrecision] + N[(a * N[(N[Log[N[(1.0 - z), $MachinePrecision]], $MachinePrecision] - b), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]

↓

code[x_, y_, z_, t_, a_, b_] := Block[{t$95$1 = N[(N[(y * N[(N[Log[z], $MachinePrecision] - t), $MachinePrecision]), $MachinePrecision] + N[(a * N[(N[Log[N[(1.0 - z), $MachinePrecision]], $MachinePrecision] - b), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[t$95$1, -1e-6], N[(x * N[Exp[t$95$1], $MachinePrecision]), $MachinePrecision], N[(x * N[Exp[N[(a * N[(N[Log[1 + (-z)], $MachinePrecision] - b), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]]

Derivation

1. Split input into 2 regimes

2. if (+.f64 (*.f64 y (-.f64 (log.f64 z) t)) (*.f64 a (-.f64 (log.f64 (-.f64 1 z)) b))) < -9.99999999999999955e-7

   1. Initial program 0.1

      \[x \cdot e^{y \cdot \left(\log z - t\right) + a \cdot \left(\log \left(1 - z\right) - b\right)} \]

   if -9.99999999999999955e-7 < (+.f64 (*.f64 y (-.f64 (log.f64 z) t)) (*.f64 a (-.f64 (log.f64 (-.f64 1 z)) b)))

   1. Initial program 6.3

      \[x \cdot e^{y \cdot \left(\log z - t\right) + a \cdot \left(\log \left(1 - z\right) - b\right)} \]

   2. Simplified 0.8

      \[\leadsto \color{blue}{x \cdot e^{\mathsf{fma}\left(a, \mathsf{log1p}\left(-z\right) - b, y \cdot \left(\log z - t\right)\right)}} \]
      Proof
      (*.f64 x (exp.f64 (fma.f64 a (-.f64 (log1p.f64 (neg.f64 z)) b) (*.f64 y (-.f64 (log.f64 z) t))))): 0 points increase in error, 0 points decrease in error
      (*.f64 x (exp.f64 (fma.f64 a (-.f64 (Rewrite<= log1p-def_binary64 (log.f64 (+.f64 1 (neg.f64 z)))) b) (*.f64 y (-.f64 (log.f64 z) t))))): 9 points increase in error, 0 points decrease in error
      (*.f64 x (exp.f64 (fma.f64 a (-.f64 (log.f64 (Rewrite<= sub-neg_binary64 (-.f64 1 z))) b) (*.f64 y (-.f64 (log.f64 z) t))))): 0 points increase in error, 0 points decrease in error
      (*.f64 x (exp.f64 (Rewrite<= fma-def_binary64 (+.f64 (*.f64 a (-.f64 (log.f64 (-.f64 1 z)) b)) (*.f64 y (-.f64 (log.f64 z) t)))))): 0 points increase in error, 0 points decrease in error
      (*.f64 x (exp.f64 (Rewrite<= +-commutative_binary64 (+.f64 (*.f64 y (-.f64 (log.f64 z) t)) (*.f64 a (-.f64 (log.f64 (-.f64 1 z)) b)))))): 0 points increase in error, 0 points decrease in error

   3. Taylor expanded in a around inf 6.9

      \[\leadsto x \cdot e^{\color{blue}{\left(\log \left(1 - z\right) - b\right) \cdot a}} \]

   4. Simplified 2.2

      \[\leadsto x \cdot e^{\color{blue}{\left(\mathsf{log1p}\left(-z\right) - b\right) \cdot a}} \]
      Proof
      (*.f64 (-.f64 (log1p.f64 (neg.f64 z)) b) a): 0 points increase in error, 0 points decrease in error
      (*.f64 (-.f64 (Rewrite<= log1p-def_binary64 (log.f64 (+.f64 1 (neg.f64 z)))) b) a): 12 points increase in error, 53 points decrease in error
      (*.f64 (-.f64 (log.f64 (Rewrite<= sub-neg_binary64 (-.f64 1 z))) b) a): 0 points increase in error, 0 points decrease in error
3. Recombined 2 regimes into one program.

4. Final simplification 0.8

   \[\leadsto \begin{array}{l} \mathbf{if}\;y \cdot \left(\log z - t\right) + a \cdot \left(\log \left(1 - z\right) - b\right) \leq -1 \cdot 10^{-6}:\\ \;\;\;\;x \cdot e^{y \cdot \left(\log z - t\right) + a \cdot \left(\log \left(1 - z\right) - b\right)}\\ \mathbf{else}:\\ \;\;\;\;x \cdot e^{a \cdot \left(\mathsf{log1p}\left(-z\right) - b\right)}\\ \end{array} \]

Alternatives

Alternative 1
Error	0.3
Cost	26368

\[x \cdot e^{\mathsf{fma}\left(a, \mathsf{log1p}\left(-z\right) - b, y \cdot \left(\log z - t\right)\right)} \]

Alternative 2
Error	6.3
Cost	13576

\[\begin{array}{l} \mathbf{if}\;y \leq -4.6 \cdot 10^{+42}:\\ \;\;\;\;x \cdot e^{y \cdot \left(-t\right)}\\ \mathbf{elif}\;y \leq 5.4 \cdot 10^{-15}:\\ \;\;\;\;x \cdot e^{a \cdot \left(\mathsf{log1p}\left(-z\right) - b\right)}\\ \mathbf{else}:\\ \;\;\;\;x \cdot {z}^{y}\\ \end{array} \]

Alternative 3
Error	8.0
Cost	6984

\[\begin{array}{l} \mathbf{if}\;y \leq -4.6 \cdot 10^{+42}:\\ \;\;\;\;x \cdot e^{y \cdot \left(-t\right)}\\ \mathbf{elif}\;y \leq 8.5 \cdot 10^{-16}:\\ \;\;\;\;\frac{x}{e^{a \cdot b}}\\ \mathbf{else}:\\ \;\;\;\;x \cdot {z}^{y}\\ \end{array} \]

Alternative 4
Error	10.7
Cost	6852

\[\begin{array}{l} \mathbf{if}\;y \leq 8.5 \cdot 10^{-16}:\\ \;\;\;\;\frac{x}{e^{a \cdot b}}\\ \mathbf{else}:\\ \;\;\;\;x \cdot {z}^{y}\\ \end{array} \]

Alternative 5
Error	19.2
Cost	6788

\[\begin{array}{l} \mathbf{if}\;y \leq 5.4 \cdot 10^{-15}:\\ \;\;\;\;\frac{x}{1 + a \cdot b}\\ \mathbf{else}:\\ \;\;\;\;x \cdot {z}^{y}\\ \end{array} \]

Alternative 6
Error	35.0
Cost	712

\[\begin{array}{l} t_1 := \frac{x}{1 + a \cdot b}\\ \mathbf{if}\;b \leq -1.16 \cdot 10^{-144}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;b \leq 0:\\ \;\;\;\;a \cdot \left(x \cdot \left(-b\right)\right)\\ \mathbf{else}:\\ \;\;\;\;t_1\\ \end{array} \]

Alternative 7
Error	27.6
Cost	708

\[\begin{array}{l} \mathbf{if}\;y \leq 8.4 \cdot 10^{-16}:\\ \;\;\;\;\frac{x}{1 + a \cdot b}\\ \mathbf{else}:\\ \;\;\;\;\left(1 + a \cdot \left(x \cdot b\right)\right) + -1\\ \end{array} \]

Alternative 8
Error	39.8
Cost	648

\[\begin{array}{l} t_1 := a \cdot \left(x \cdot \left(-b\right)\right)\\ \mathbf{if}\;y \leq -2.8 \cdot 10^{+14}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;y \leq 5.3 \cdot 10^{-15}:\\ \;\;\;\;x\\ \mathbf{else}:\\ \;\;\;\;t_1\\ \end{array} \]

Alternative 9
Error	44.6
Cost	64

\[x \]

Reproduce

herbie shell --seed 2022317 
(FPCore (x y z t a b)
  :name "Numeric.SpecFunctions:incompleteBetaApprox from math-functions-0.1.5.2, B"
  :precision binary64
  (* x (exp (+ (* y (- (log z) t)) (* a (- (log (- 1.0 z)) b))))))