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

Average Error: 20.0 → 5.2

Time: 5.5s

Precision: binary64

\[ \begin{array}{c}[x, y] = \mathsf{sort}([x, y])\\ \end{array} \]

Math

\[\frac{x \cdot y}{\left(\left(x + y\right) \cdot \left(x + y\right)\right) \cdot \left(\left(x + y\right) + 1\right)} \]

↓

\[\begin{array}{l} \mathbf{if}\;x \leq -3.3 \cdot 10^{+104}:\\ \;\;\;\;\frac{y}{{x}^{2}}\\ \mathbf{elif}\;x \leq -8 \cdot 10^{-24}:\\ \;\;\;\;y \cdot \frac{x}{{\left(x + y\right)}^{3} + {\left(x + y\right)}^{2}}\\ \mathbf{else}:\\ \;\;\;\;\frac{x}{\left(\mathsf{fma}\left(y, y, \frac{x}{\frac{y}{x}}\right) + \left(\mathsf{fma}\left(3, x \cdot x, x \cdot 2\right) + \frac{{x}^{3}}{y}\right)\right) + \mathsf{fma}\left(3, x \cdot y, y\right)}\\ \end{array} \]

FPCore

(FPCore (x y)
 :precision binary64
 (/ (* x y) (* (* (+ x y) (+ x y)) (+ (+ x y) 1.0))))

↓

(FPCore (x y)
 :precision binary64
 (if (<= x -3.3e+104)
   (/ y (pow x 2.0))
   (if (<= x -8e-24)
     (* y (/ x (+ (pow (+ x y) 3.0) (pow (+ x y) 2.0))))
     (/
      x
      (+
       (+
        (fma y y (/ x (/ y x)))
        (+ (fma 3.0 (* x x) (* x 2.0)) (/ (pow x 3.0) y)))
       (fma 3.0 (* x y) y))))))

C

double code(double x, double y) {
	return (x * y) / (((x + y) * (x + y)) * ((x + y) + 1.0));
}

↓

double code(double x, double y) {
	double tmp;
	if (x <= -3.3e+104) {
		tmp = y / pow(x, 2.0);
	} else if (x <= -8e-24) {
		tmp = y * (x / (pow((x + y), 3.0) + pow((x + y), 2.0)));
	} else {
		tmp = x / ((fma(y, y, (x / (y / x))) + (fma(3.0, (x * x), (x * 2.0)) + (pow(x, 3.0) / y))) + fma(3.0, (x * y), y));
	}
	return tmp;
}

Julia

function code(x, y)
	return Float64(Float64(x * y) / Float64(Float64(Float64(x + y) * Float64(x + y)) * Float64(Float64(x + y) + 1.0)))
end

↓

function code(x, y)
	tmp = 0.0
	if (x <= -3.3e+104)
		tmp = Float64(y / (x ^ 2.0));
	elseif (x <= -8e-24)
		tmp = Float64(y * Float64(x / Float64((Float64(x + y) ^ 3.0) + (Float64(x + y) ^ 2.0))));
	else
		tmp = Float64(x / Float64(Float64(fma(y, y, Float64(x / Float64(y / x))) + Float64(fma(3.0, Float64(x * x), Float64(x * 2.0)) + Float64((x ^ 3.0) / y))) + fma(3.0, Float64(x * y), y)));
	end
	return tmp
end

Wolfram

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

↓

code[x_, y_] := If[LessEqual[x, -3.3e+104], N[(y / N[Power[x, 2.0], $MachinePrecision]), $MachinePrecision], If[LessEqual[x, -8e-24], N[(y * N[(x / N[(N[Power[N[(x + y), $MachinePrecision], 3.0], $MachinePrecision] + N[Power[N[(x + y), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(x / N[(N[(N[(y * y + N[(x / N[(y / x), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + N[(N[(3.0 * N[(x * x), $MachinePrecision] + N[(x * 2.0), $MachinePrecision]), $MachinePrecision] + N[(N[Power[x, 3.0], $MachinePrecision] / y), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + N[(3.0 * N[(x * y), $MachinePrecision] + y), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Error

Bits error versus x

Bits error versus y

Target

Original	20.0
Target	0.1
Herbie	5.2

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

Derivation

1. Split input into 3 regimes

2. if x < -3.29999999999999985e104

   1. Initial program 25.6

      \[\frac{x \cdot y}{\left(\left(x + y\right) \cdot \left(x + y\right)\right) \cdot \left(\left(x + y\right) + 1\right)} \]

   2. Simplified 13.6

      \[\leadsto \color{blue}{\frac{x}{\frac{\mathsf{fma}\left(x + y, x + y, {\left(x + y\right)}^{3}\right)}{y}}} \]

   3. Taylor expanded in x around inf 10.9

      \[\leadsto \color{blue}{\frac{y}{{x}^{2}}} \]

   if -3.29999999999999985e104 < x < -7.99999999999999939e-24

   1. Initial program 8.3

      \[\frac{x \cdot y}{\left(\left(x + y\right) \cdot \left(x + y\right)\right) \cdot \left(\left(x + y\right) + 1\right)} \]

   2. Simplified 10.9

      \[\leadsto \color{blue}{\frac{x}{\frac{\mathsf{fma}\left(x + y, x + y, {\left(x + y\right)}^{3}\right)}{y}}} \]

   3. Applied egg-rr 5.1

      \[\leadsto \color{blue}{y \cdot \frac{x}{{\left(x + y\right)}^{3} + {\left(x + y\right)}^{2}}} \]

   if -7.99999999999999939e-24 < x

   1. Initial program 20.4

      \[\frac{x \cdot y}{\left(\left(x + y\right) \cdot \left(x + y\right)\right) \cdot \left(\left(x + y\right) + 1\right)} \]

   2. Simplified 10.9

      \[\leadsto \color{blue}{\frac{x}{\frac{\mathsf{fma}\left(x + y, x + y, {\left(x + y\right)}^{3}\right)}{y}}} \]

   3. Taylor expanded in x around 0 5.5

      \[\leadsto \frac{x}{\color{blue}{3 \cdot \left(y \cdot x\right) + \left(\frac{{x}^{3}}{y} + \left(y + \left({y}^{2} + \left(\frac{{x}^{2}}{y} + \left(3 \cdot {x}^{2} + 2 \cdot x\right)\right)\right)\right)\right)}} \]

   4. Simplified 2.1

      \[\leadsto \frac{x}{\color{blue}{\left(\mathsf{fma}\left(y, y, \frac{x}{\frac{y}{x}}\right) + \left(\mathsf{fma}\left(3, x \cdot x, 2 \cdot x\right) + \frac{{x}^{3}}{y}\right)\right) + \mathsf{fma}\left(3, y \cdot x, y\right)}} \]
3. Recombined 3 regimes into one program.

4. Final simplification 5.2

   \[\leadsto \begin{array}{l} \mathbf{if}\;x \leq -3.3 \cdot 10^{+104}:\\ \;\;\;\;\frac{y}{{x}^{2}}\\ \mathbf{elif}\;x \leq -8 \cdot 10^{-24}:\\ \;\;\;\;y \cdot \frac{x}{{\left(x + y\right)}^{3} + {\left(x + y\right)}^{2}}\\ \mathbf{else}:\\ \;\;\;\;\frac{x}{\left(\mathsf{fma}\left(y, y, \frac{x}{\frac{y}{x}}\right) + \left(\mathsf{fma}\left(3, x \cdot x, x \cdot 2\right) + \frac{{x}^{3}}{y}\right)\right) + \mathsf{fma}\left(3, x \cdot y, y\right)}\\ \end{array} \]

Reproduce

herbie shell --seed 2022162 
(FPCore (x y)
  :name "Numeric.SpecFunctions:incompleteBetaApprox from math-functions-0.1.5.2, A"
  :precision binary64

  :herbie-target
  (/ (/ (/ x (+ (+ y 1.0) x)) (+ y x)) (/ 1.0 (/ y (+ y x))))

  (/ (* x y) (* (* (+ x y) (+ x y)) (+ (+ x y) 1.0))))