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

Average Error: 0.0 → 0.1

Time: 8.0s

Precision: binary64

Cost: 704

Math

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

↓

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

FPCore

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

↓

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

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 / ((1.0 / y) + (x * 0.5)));
}

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) / ((1.0d0 / y) + (x * 0.5d0)))
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 / ((1.0 / y) + (x * 0.5)));
}

Python

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

↓

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

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(1.0 / y) + Float64(x * 0.5))))
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 / ((1.0 / y) + (x * 0.5)));
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[(1.0 / y), $MachinePrecision] + N[(x * 0.5), $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 (Rewrite<= fma-def_binary64 (+.f64 (*.f64 x (/.f64 1 2)) (/.f64 1 y))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (+.f64 (Rewrite<= *-commutative_binary64 (*.f64 (/.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)))): 4 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (+.f64 (Rewrite=> *-lft-identity_binary64 (/.f64 x 2)) (/.f64 1 y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (Rewrite=> +-commutative_binary64 (+.f64 (/.f64 1 y) (/.f64 x 2))))): 0 points increase in error, 4 points decrease in error
   (+.f64 x (/.f64 -1 (+.f64 (/.f64 1 y) (/.f64 (Rewrite<= remove-double-neg_binary64 (neg.f64 (neg.f64 x))) 2)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (+.f64 (/.f64 1 y) (Rewrite=> distribute-frac-neg_binary64 (neg.f64 (/.f64 (neg.f64 x) 2)))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (Rewrite<= sub-neg_binary64 (-.f64 (/.f64 1 y) (/.f64 (neg.f64 x) 2))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (Rewrite<= /-rgt-identity_binary64 (/.f64 (/.f64 (neg.f64 x) 2) 1))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (/.f64 (/.f64 (neg.f64 x) 2) (Rewrite<= *-inverses_binary64 (/.f64 (neg.f64 y) (neg.f64 y))))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (Rewrite<= associate-/l*_binary64 (/.f64 (*.f64 (/.f64 (neg.f64 x) 2) (neg.f64 y)) (neg.f64 y)))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (/.f64 (*.f64 (/.f64 (neg.f64 x) 2) (neg.f64 y)) (Rewrite<= /-rgt-identity_binary64 (/.f64 (neg.f64 y) 1)))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (/.f64 (*.f64 (/.f64 (neg.f64 x) 2) (neg.f64 y)) (/.f64 (Rewrite=> neg-mul-1_binary64 (*.f64 -1 y)) 1))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (/.f64 (*.f64 (/.f64 (neg.f64 x) 2) (neg.f64 y)) (/.f64 (Rewrite=> *-commutative_binary64 (*.f64 y -1)) 1))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (/.f64 (*.f64 (/.f64 (neg.f64 x) 2) (neg.f64 y)) (Rewrite=> associate-/l*_binary64 (/.f64 y (/.f64 1 -1))))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (/.f64 (*.f64 (/.f64 (neg.f64 x) 2) (neg.f64 y)) (/.f64 y (Rewrite=> metadata-eval -1)))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (Rewrite=> associate-/l*_binary64 (/.f64 (/.f64 (neg.f64 x) 2) (/.f64 (/.f64 y -1) (neg.f64 y))))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (/.f64 (/.f64 (neg.f64 x) 2) (Rewrite<= associate-/r*_binary64 (/.f64 y (*.f64 -1 (neg.f64 y)))))))): 4 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (/.f64 (/.f64 (neg.f64 x) 2) (/.f64 y (Rewrite=> mul-1-neg_binary64 (neg.f64 (neg.f64 y)))))))): 0 points increase in error, 4 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (/.f64 (/.f64 (neg.f64 x) 2) (/.f64 y (Rewrite=> remove-double-neg_binary64 y)))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (Rewrite<= associate-/l*_binary64 (/.f64 (*.f64 (/.f64 (neg.f64 x) 2) y) y))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (-.f64 (/.f64 1 y) (/.f64 (Rewrite<= associate-/r/_binary64 (/.f64 (neg.f64 x) (/.f64 2 y))) y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (Rewrite<= div-sub_binary64 (/.f64 (-.f64 1 (/.f64 (neg.f64 x) (/.f64 2 y))) y)))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (/.f64 (Rewrite<= unsub-neg_binary64 (+.f64 1 (neg.f64 (/.f64 (neg.f64 x) (/.f64 2 y))))) y))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (/.f64 (+.f64 1 (Rewrite<= distribute-frac-neg_binary64 (/.f64 (neg.f64 (neg.f64 x)) (/.f64 2 y)))) y))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (/.f64 (+.f64 1 (/.f64 (Rewrite=> remove-double-neg_binary64 x) (/.f64 2 y))) y))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (/.f64 -1 (/.f64 (+.f64 1 (Rewrite<= associate-/l*_binary64 (/.f64 (*.f64 x y) 2))) y))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (Rewrite<= associate-/l*_binary64 (/.f64 (*.f64 -1 y) (+.f64 1 (/.f64 (*.f64 x y) 2))))): 0 points increase in error, 0 points decrease in error
   (+.f64 x (Rewrite<= associate-*r/_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. Applied egg-rr 0.1

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

4. Final simplification 0.1

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

Alternatives

Alternative 1
Error	5.8
Cost	585

\[\begin{array}{l} \mathbf{if}\;y \leq -5.5 \cdot 10^{+143} \lor \neg \left(y \leq 3.1 \cdot 10^{+88}\right):\\ \;\;\;\;x + \frac{-2}{x}\\ \mathbf{else}:\\ \;\;\;\;x - y\\ \end{array} \]

Alternative 2
Error	8.0
Cost	456

\[\begin{array}{l} \mathbf{if}\;x \leq -2.9 \cdot 10^{-6}:\\ \;\;\;\;x\\ \mathbf{elif}\;x \leq 1.4:\\ \;\;\;\;x - y\\ \mathbf{else}:\\ \;\;\;\;x\\ \end{array} \]

Alternative 3
Error	23.7
Cost	64

\[x \]

Reproduce

herbie shell --seed 2022343 
(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)))))