Data.Number.Erf:$cinvnormcdf from erf-2.0.0.0, B

Average Error: 0.0 → 0.1

Time: 8.5s

Precision: binary64

Cost: 704

Math

\[x - \frac{y}{1 + \frac{x \cdot y}{2}} \]

↓

\[x + \frac{1}{x \cdot -0.5 + \frac{-1}{y}} \]

FPCore

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

↓

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

C

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

↓

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

Fortran

real(8) function code(x, y)
    real(8), intent (in) :: x
    real(8), intent (in) :: y
    code = x - (y / (1.0d0 + ((x * y) / 2.0d0)))
end function

↓

real(8) function code(x, y)
    real(8), intent (in) :: x
    real(8), intent (in) :: y
    code = x + (1.0d0 / ((x * (-0.5d0)) + ((-1.0d0) / y)))
end function

Java

public static double code(double x, double y) {
	return x - (y / (1.0 + ((x * y) / 2.0)));
}

↓

public static double code(double x, double y) {
	return x + (1.0 / ((x * -0.5) + (-1.0 / y)));
}

Python

def code(x, y):
	return x - (y / (1.0 + ((x * y) / 2.0)))

↓

def code(x, y):
	return x + (1.0 / ((x * -0.5) + (-1.0 / y)))

Julia

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

↓

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

MATLAB

function tmp = code(x, y)
	tmp = x - (y / (1.0 + ((x * y) / 2.0)));
end

↓

function tmp = code(x, y)
	tmp = x + (1.0 / ((x * -0.5) + (-1.0 / y)));
end

Wolfram

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

↓

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

Derivation

1. Initial program 0.0

   \[x - \frac{y}{1 + \frac{x \cdot y}{2}} \]

2. Simplified 0.1

   \[\leadsto \color{blue}{x + \frac{1}{\mathsf{fma}\left(x, -0.5, \frac{-1}{y}\right)}} \]
   Proof
   (+.f64 x (/.f64 1 (fma.f64 x -1/2 (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (fma.f64 x (Rewrite<= metadata-eval (/.f64 1 -2)) (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (fma.f64 x (/.f64 1 (Rewrite<= metadata-eval (/.f64 2 -1))) (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (Rewrite<= fma-def_binary64 (+.f64 (*.f64 x (/.f64 1 (/.f64 2 -1))) (/.f64 -1 y))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (+.f64 (Rewrite<= *-commutative_binary64 (*.f64 (/.f64 1 (/.f64 2 -1)) x)) (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (+.f64 (*.f64 (/.f64 1 (Rewrite=> metadata-eval -2)) x) (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (+.f64 (*.f64 (Rewrite=> metadata-eval -1/2) x) (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (+.f64 (*.f64 (Rewrite<= metadata-eval (/.f64 -1 2)) x) (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (+.f64 (Rewrite<= associate-/r/_binary64 (/.f64 -1 (/.f64 2 x))) (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (+.f64 (Rewrite<= associate-/l*_binary64 (/.f64 (*.f64 -1 x) 2)) (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (+.f64 (Rewrite<= associate-*r/_binary64 (*.f64 -1 (/.f64 x 2))) (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (Rewrite=> fma-def_binary64 (fma.f64 -1 (/.f64 x 2) (/.f64 -1 y))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (fma.f64 (Rewrite<= metadata-eval (neg.f64 1)) (/.f64 x 2) (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (fma.f64 (neg.f64 (Rewrite<= *-inverses_binary64 (/.f64 y y))) (/.f64 x 2) (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (fma.f64 (Rewrite<= distribute-frac-neg_binary64 (/.f64 (neg.f64 y) y)) (/.f64 x 2) (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (fma.f64 (/.f64 (neg.f64 y) y) (Rewrite<= /-rgt-identity_binary64 (/.f64 (/.f64 x 2) 1)) (/.f64 -1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (fma.f64 (/.f64 (neg.f64 y) y) (/.f64 (/.f64 x 2) 1) (/.f64 (Rewrite<= metadata-eval (neg.f64 1)) y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (fma.f64 (/.f64 (neg.f64 y) y) (/.f64 (/.f64 x 2) 1) (Rewrite<= distribute-neg-frac_binary64 (neg.f64 (/.f64 1 y)))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (Rewrite<= fma-neg_binary64 (-.f64 (*.f64 (/.f64 (neg.f64 y) y) (/.f64 (/.f64 x 2) 1)) (/.f64 1 y))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (-.f64 (Rewrite<= times-frac_binary64 (/.f64 (*.f64 (neg.f64 y) (/.f64 x 2)) (*.f64 y 1))) (/.f64 1 y)))): 2 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (-.f64 (/.f64 (Rewrite<= distribute-lft-neg-in_binary64 (neg.f64 (*.f64 y (/.f64 x 2)))) (*.f64 y 1)) (/.f64 1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (-.f64 (/.f64 (neg.f64 (Rewrite<= *-commutative_binary64 (*.f64 (/.f64 x 2) y))) (*.f64 y 1)) (/.f64 1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (-.f64 (/.f64 (neg.f64 (Rewrite<= associate-/r/_binary64 (/.f64 x (/.f64 2 y)))) (*.f64 y 1)) (/.f64 1 y)))): 3 points increase in error, 1 points decrease in error
   (+.f64 x (/.f64 1 (-.f64 (/.f64 (Rewrite<= distribute-frac-neg_binary64 (/.f64 (neg.f64 x) (/.f64 2 y))) (*.f64 y 1)) (/.f64 1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (-.f64 (/.f64 (/.f64 (neg.f64 x) (/.f64 2 y)) (Rewrite=> *-rgt-identity_binary64 y)) (/.f64 1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 1 (Rewrite<= div-sub_binary64 (/.f64 (-.f64 (/.f64 (neg.f64 x) (/.f64 2 y)) 1) y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (Rewrite<= associate-/l*_binary64 (/.f64 (*.f64 1 y) (-.f64 (/.f64 (neg.f64 x) (/.f64 2 y)) 1)))): 3 points increase in error, 7 points decrease in error
   (+.f64 x (/.f64 (Rewrite=> *-lft-identity_binary64 y) (-.f64 (/.f64 (neg.f64 x) (/.f64 2 y)) 1))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 y (Rewrite=> sub-neg_binary64 (+.f64 (/.f64 (neg.f64 x) (/.f64 2 y)) (neg.f64 1))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 y (+.f64 (Rewrite=> distribute-frac-neg_binary64 (neg.f64 (/.f64 x (/.f64 2 y)))) (neg.f64 1)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 y (+.f64 (neg.f64 (Rewrite<= associate-/l*_binary64 (/.f64 (*.f64 x y) 2))) (neg.f64 1)))): 4 points increase in error, 4 points decrease in error
   (+.f64 x (/.f64 y (Rewrite<= distribute-neg-in_binary64 (neg.f64 (+.f64 (/.f64 (*.f64 x y) 2) 1))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 y (neg.f64 (Rewrite<= +-commutative_binary64 (+.f64 1 (/.f64 (*.f64 x y) 2)))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 y (Rewrite=> neg-mul-1_binary64 (*.f64 -1 (+.f64 1 (/.f64 (*.f64 x y) 2)))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (Rewrite<= associate-/l/_binary64 (/.f64 (/.f64 y (+.f64 1 (/.f64 (*.f64 x y) 2))) -1))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 (/.f64 y (+.f64 1 (/.f64 (*.f64 x y) 2))) (Rewrite<= metadata-eval (/.f64 1 -1)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (Rewrite<= associate-/l*_binary64 (/.f64 (*.f64 (/.f64 y (+.f64 1 (/.f64 (*.f64 x y) 2))) -1) 1))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 (Rewrite<= *-commutative_binary64 (*.f64 -1 (/.f64 y (+.f64 1 (/.f64 (*.f64 x y) 2))))) 1)): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 (Rewrite<= neg-mul-1_binary64 (neg.f64 (/.f64 y (+.f64 1 (/.f64 (*.f64 x y) 2))))) 1)): 0 points increase in error, 0 points decrease in error
   (+.f64 x (Rewrite=> /-rgt-identity_binary64 (neg.f64 (/.f64 y (+.f64 1 (/.f64 (*.f64 x y) 2)))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (Rewrite=> neg-mul-1_binary64 (*.f64 -1 (/.f64 y (+.f64 1 (/.f64 (*.f64 x y) 2)))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (Rewrite=> *-commutative_binary64 (*.f64 (/.f64 y (+.f64 1 (/.f64 (*.f64 x y) 2))) -1))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (Rewrite<= *-commutative_binary64 (*.f64 -1 (/.f64 y (+.f64 1 (/.f64 (*.f64 x y) 2)))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (Rewrite<= neg-mul-1_binary64 (neg.f64 (/.f64 y (+.f64 1 (/.f64 (*.f64 x y) 2)))))): 0 points increase in error, 0 points decrease in error
   (Rewrite<= sub-neg_binary64 (-.f64 x (/.f64 y (+.f64 1 (/.f64 (*.f64 x y) 2))))): 0 points increase in error, 0 points decrease in error

3. Taylor expanded in x around 0 0.1

   \[\leadsto x + \frac{1}{\color{blue}{-0.5 \cdot x - \frac{1}{y}}} \]

4. Final simplification 0.1

   \[\leadsto x + \frac{1}{x \cdot -0.5 + \frac{-1}{y}} \]

Alternatives

Alternative 1
Error	12.9
Cost	656

\[\begin{array}{l} \mathbf{if}\;x \leq -2.2581364671467762 \cdot 10^{-7}:\\ \;\;\;\;x\\ \mathbf{elif}\;x \leq -6.3317407257963676 \cdot 10^{-49}:\\ \;\;\;\;\frac{-2}{x}\\ \mathbf{elif}\;x \leq -1.6064033313862016 \cdot 10^{-148}:\\ \;\;\;\;x\\ \mathbf{elif}\;x \leq 1.086368859771886 \cdot 10^{-100}:\\ \;\;\;\;-y\\ \mathbf{else}:\\ \;\;\;\;x\\ \end{array} \]

Alternative 2
Error	8.0
Cost	588

\[\begin{array}{l} \mathbf{if}\;x \leq -2.2581364671467762 \cdot 10^{-7}:\\ \;\;\;\;x\\ \mathbf{elif}\;x \leq -6.3317407257963676 \cdot 10^{-49}:\\ \;\;\;\;\frac{-2}{x}\\ \mathbf{elif}\;x \leq 1.256492122538902 \cdot 10^{-6}:\\ \;\;\;\;x - y\\ \mathbf{else}:\\ \;\;\;\;x\\ \end{array} \]

Alternative 3
Error	5.6
Cost	584

\[\begin{array}{l} t_0 := x + \frac{-2}{x}\\ \mathbf{if}\;y \leq -3.339674474706417 \cdot 10^{+92}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;y \leq 4.356904807041854 \cdot 10^{+115}:\\ \;\;\;\;x - y\\ \mathbf{else}:\\ \;\;\;\;t_0\\ \end{array} \]

Alternative 4
Error	12.9
Cost	392

\[\begin{array}{l} \mathbf{if}\;x \leq -1.6064033313862016 \cdot 10^{-148}:\\ \;\;\;\;x\\ \mathbf{elif}\;x \leq 1.086368859771886 \cdot 10^{-100}:\\ \;\;\;\;-y\\ \mathbf{else}:\\ \;\;\;\;x\\ \end{array} \]

Alternative 5
Error	23.5
Cost	64

\[x \]

Reproduce

herbie shell --seed 2022318 
(FPCore (x y)
  :name "Data.Number.Erf:$cinvnormcdf from erf-2.0.0.0, B"
  :precision binary64
  (- x (/ y (+ 1.0 (/ (* x y) 2.0)))))