
(FPCore (a b) :precision binary64 (- (+ (pow (+ (* a a) (* b b)) 2.0) (* 4.0 (+ (* (* a a) (- 1.0 a)) (* (* b b) (+ 3.0 a))))) 1.0))
double code(double a, double b) {
return (pow(((a * a) + (b * b)), 2.0) + (4.0 * (((a * a) * (1.0 - a)) + ((b * b) * (3.0 + a))))) - 1.0;
}
real(8) function code(a, b)
real(8), intent (in) :: a
real(8), intent (in) :: b
code = ((((a * a) + (b * b)) ** 2.0d0) + (4.0d0 * (((a * a) * (1.0d0 - a)) + ((b * b) * (3.0d0 + a))))) - 1.0d0
end function
public static double code(double a, double b) {
return (Math.pow(((a * a) + (b * b)), 2.0) + (4.0 * (((a * a) * (1.0 - a)) + ((b * b) * (3.0 + a))))) - 1.0;
}
def code(a, b): return (math.pow(((a * a) + (b * b)), 2.0) + (4.0 * (((a * a) * (1.0 - a)) + ((b * b) * (3.0 + a))))) - 1.0
function code(a, b) return Float64(Float64((Float64(Float64(a * a) + Float64(b * b)) ^ 2.0) + Float64(4.0 * Float64(Float64(Float64(a * a) * Float64(1.0 - a)) + Float64(Float64(b * b) * Float64(3.0 + a))))) - 1.0) end
function tmp = code(a, b) tmp = ((((a * a) + (b * b)) ^ 2.0) + (4.0 * (((a * a) * (1.0 - a)) + ((b * b) * (3.0 + a))))) - 1.0; end
code[a_, b_] := N[(N[(N[Power[N[(N[(a * a), $MachinePrecision] + N[(b * b), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] + N[(4.0 * N[(N[(N[(a * a), $MachinePrecision] * N[(1.0 - a), $MachinePrecision]), $MachinePrecision] + N[(N[(b * b), $MachinePrecision] * N[(3.0 + a), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] - 1.0), $MachinePrecision]
\begin{array}{l}
\\
\left({\left(a \cdot a + b \cdot b\right)}^{2} + 4 \cdot \left(\left(a \cdot a\right) \cdot \left(1 - a\right) + \left(b \cdot b\right) \cdot \left(3 + a\right)\right)\right) - 1
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 12 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (a b) :precision binary64 (- (+ (pow (+ (* a a) (* b b)) 2.0) (* 4.0 (+ (* (* a a) (- 1.0 a)) (* (* b b) (+ 3.0 a))))) 1.0))
double code(double a, double b) {
return (pow(((a * a) + (b * b)), 2.0) + (4.0 * (((a * a) * (1.0 - a)) + ((b * b) * (3.0 + a))))) - 1.0;
}
real(8) function code(a, b)
real(8), intent (in) :: a
real(8), intent (in) :: b
code = ((((a * a) + (b * b)) ** 2.0d0) + (4.0d0 * (((a * a) * (1.0d0 - a)) + ((b * b) * (3.0d0 + a))))) - 1.0d0
end function
public static double code(double a, double b) {
return (Math.pow(((a * a) + (b * b)), 2.0) + (4.0 * (((a * a) * (1.0 - a)) + ((b * b) * (3.0 + a))))) - 1.0;
}
def code(a, b): return (math.pow(((a * a) + (b * b)), 2.0) + (4.0 * (((a * a) * (1.0 - a)) + ((b * b) * (3.0 + a))))) - 1.0
function code(a, b) return Float64(Float64((Float64(Float64(a * a) + Float64(b * b)) ^ 2.0) + Float64(4.0 * Float64(Float64(Float64(a * a) * Float64(1.0 - a)) + Float64(Float64(b * b) * Float64(3.0 + a))))) - 1.0) end
function tmp = code(a, b) tmp = ((((a * a) + (b * b)) ^ 2.0) + (4.0 * (((a * a) * (1.0 - a)) + ((b * b) * (3.0 + a))))) - 1.0; end
code[a_, b_] := N[(N[(N[Power[N[(N[(a * a), $MachinePrecision] + N[(b * b), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] + N[(4.0 * N[(N[(N[(a * a), $MachinePrecision] * N[(1.0 - a), $MachinePrecision]), $MachinePrecision] + N[(N[(b * b), $MachinePrecision] * N[(3.0 + a), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] - 1.0), $MachinePrecision]
\begin{array}{l}
\\
\left({\left(a \cdot a + b \cdot b\right)}^{2} + 4 \cdot \left(\left(a \cdot a\right) \cdot \left(1 - a\right) + \left(b \cdot b\right) \cdot \left(3 + a\right)\right)\right) - 1
\end{array}
(FPCore (a b)
:precision binary64
(let* ((t_0 (pow (sqrt (fma b b (* a a))) -1.0))
(t_1 (fma a a (* b b)))
(t_2 (* 4.0 (+ (* (* a a) (- 1.0 a)) (* (* b b) (+ a 3.0))))))
(if (<= (+ (pow (+ (* a a) (* b b)) 2.0) t_2) 4e+278)
(+ (+ t_2 (/ t_1 (/ 1.0 t_1))) -1.0)
(+ -1.0 (+ (/ t_1 (* t_0 t_0)) (* 4.0 (* (* b b) 3.0)))))))
double code(double a, double b) {
double t_0 = pow(sqrt(fma(b, b, (a * a))), -1.0);
double t_1 = fma(a, a, (b * b));
double t_2 = 4.0 * (((a * a) * (1.0 - a)) + ((b * b) * (a + 3.0)));
double tmp;
if ((pow(((a * a) + (b * b)), 2.0) + t_2) <= 4e+278) {
tmp = (t_2 + (t_1 / (1.0 / t_1))) + -1.0;
} else {
tmp = -1.0 + ((t_1 / (t_0 * t_0)) + (4.0 * ((b * b) * 3.0)));
}
return tmp;
}
function code(a, b) t_0 = sqrt(fma(b, b, Float64(a * a))) ^ -1.0 t_1 = fma(a, a, Float64(b * b)) t_2 = Float64(4.0 * Float64(Float64(Float64(a * a) * Float64(1.0 - a)) + Float64(Float64(b * b) * Float64(a + 3.0)))) tmp = 0.0 if (Float64((Float64(Float64(a * a) + Float64(b * b)) ^ 2.0) + t_2) <= 4e+278) tmp = Float64(Float64(t_2 + Float64(t_1 / Float64(1.0 / t_1))) + -1.0); else tmp = Float64(-1.0 + Float64(Float64(t_1 / Float64(t_0 * t_0)) + Float64(4.0 * Float64(Float64(b * b) * 3.0)))); end return tmp end
code[a_, b_] := Block[{t$95$0 = N[Power[N[Sqrt[N[(b * b + N[(a * a), $MachinePrecision]), $MachinePrecision]], $MachinePrecision], -1.0], $MachinePrecision]}, Block[{t$95$1 = N[(a * a + N[(b * b), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$2 = N[(4.0 * N[(N[(N[(a * a), $MachinePrecision] * N[(1.0 - a), $MachinePrecision]), $MachinePrecision] + N[(N[(b * b), $MachinePrecision] * N[(a + 3.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[N[(N[Power[N[(N[(a * a), $MachinePrecision] + N[(b * b), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] + t$95$2), $MachinePrecision], 4e+278], N[(N[(t$95$2 + N[(t$95$1 / N[(1.0 / t$95$1), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + -1.0), $MachinePrecision], N[(-1.0 + N[(N[(t$95$1 / N[(t$95$0 * t$95$0), $MachinePrecision]), $MachinePrecision] + N[(4.0 * N[(N[(b * b), $MachinePrecision] * 3.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := {\left(\sqrt{\mathsf{fma}\left(b, b, a \cdot a\right)}\right)}^{-1}\\
t_1 := \mathsf{fma}\left(a, a, b \cdot b\right)\\
t_2 := 4 \cdot \left(\left(a \cdot a\right) \cdot \left(1 - a\right) + \left(b \cdot b\right) \cdot \left(a + 3\right)\right)\\
\mathbf{if}\;{\left(a \cdot a + b \cdot b\right)}^{2} + t\_2 \leq 4 \cdot 10^{+278}:\\
\;\;\;\;\left(t\_2 + \frac{t\_1}{\frac{1}{t\_1}}\right) + -1\\
\mathbf{else}:\\
\;\;\;\;-1 + \left(\frac{t\_1}{t\_0 \cdot t\_0} + 4 \cdot \left(\left(b \cdot b\right) \cdot 3\right)\right)\\
\end{array}
\end{array}
if (+.f64 (pow.f64 (+.f64 (*.f64 a a) (*.f64 b b)) #s(literal 2 binary64)) (*.f64 #s(literal 4 binary64) (+.f64 (*.f64 (*.f64 a a) (-.f64 #s(literal 1 binary64) a)) (*.f64 (*.f64 b b) (+.f64 #s(literal 3 binary64) a))))) < 3.99999999999999985e278Initial program 99.7%
lift-pow.f64N/A
unpow2N/A
lift-+.f64N/A
flip3-+N/A
clear-numN/A
un-div-invN/A
lower-/.f64N/A
lift-+.f64N/A
lift-*.f64N/A
lower-fma.f64N/A
clear-numN/A
Applied rewrites99.8%
if 3.99999999999999985e278 < (+.f64 (pow.f64 (+.f64 (*.f64 a a) (*.f64 b b)) #s(literal 2 binary64)) (*.f64 #s(literal 4 binary64) (+.f64 (*.f64 (*.f64 a a) (-.f64 #s(literal 1 binary64) a)) (*.f64 (*.f64 b b) (+.f64 #s(literal 3 binary64) a))))) Initial program 53.5%
lift-pow.f64N/A
unpow2N/A
lift-+.f64N/A
flip3-+N/A
clear-numN/A
un-div-invN/A
lower-/.f64N/A
lift-+.f64N/A
lift-*.f64N/A
lower-fma.f64N/A
clear-numN/A
Applied rewrites53.5%
Taylor expanded in a around 0
*-commutativeN/A
lower-*.f64N/A
unpow2N/A
lower-*.f64100.0
Applied rewrites100.0%
lift-/.f64N/A
inv-powN/A
unpow1N/A
sqr-powN/A
unpow-prod-downN/A
lower-*.f64N/A
lower-pow.f64N/A
metadata-evalN/A
unpow1/2N/A
lower-sqrt.f64N/A
lift-fma.f64N/A
lift-*.f64N/A
+-commutativeN/A
lift-*.f64N/A
lower-fma.f64N/A
lower-pow.f64N/A
Applied rewrites100.0%
Final simplification99.9%
(FPCore (a b)
:precision binary64
(let* ((t_0 (fma a a (* b b)))
(t_1 (/ t_0 (/ 1.0 t_0)))
(t_2 (* 4.0 (+ (* (* a a) (- 1.0 a)) (* (* b b) (+ a 3.0))))))
(if (<= (+ (pow (+ (* a a) (* b b)) 2.0) t_2) 2e+62)
(+ (+ t_2 t_1) -1.0)
(+ -1.0 (+ t_1 (* 4.0 (* (* b b) 3.0)))))))
double code(double a, double b) {
double t_0 = fma(a, a, (b * b));
double t_1 = t_0 / (1.0 / t_0);
double t_2 = 4.0 * (((a * a) * (1.0 - a)) + ((b * b) * (a + 3.0)));
double tmp;
if ((pow(((a * a) + (b * b)), 2.0) + t_2) <= 2e+62) {
tmp = (t_2 + t_1) + -1.0;
} else {
tmp = -1.0 + (t_1 + (4.0 * ((b * b) * 3.0)));
}
return tmp;
}
function code(a, b) t_0 = fma(a, a, Float64(b * b)) t_1 = Float64(t_0 / Float64(1.0 / t_0)) t_2 = Float64(4.0 * Float64(Float64(Float64(a * a) * Float64(1.0 - a)) + Float64(Float64(b * b) * Float64(a + 3.0)))) tmp = 0.0 if (Float64((Float64(Float64(a * a) + Float64(b * b)) ^ 2.0) + t_2) <= 2e+62) tmp = Float64(Float64(t_2 + t_1) + -1.0); else tmp = Float64(-1.0 + Float64(t_1 + Float64(4.0 * Float64(Float64(b * b) * 3.0)))); end return tmp end
code[a_, b_] := Block[{t$95$0 = N[(a * a + N[(b * b), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(t$95$0 / N[(1.0 / t$95$0), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$2 = N[(4.0 * N[(N[(N[(a * a), $MachinePrecision] * N[(1.0 - a), $MachinePrecision]), $MachinePrecision] + N[(N[(b * b), $MachinePrecision] * N[(a + 3.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[N[(N[Power[N[(N[(a * a), $MachinePrecision] + N[(b * b), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] + t$95$2), $MachinePrecision], 2e+62], N[(N[(t$95$2 + t$95$1), $MachinePrecision] + -1.0), $MachinePrecision], N[(-1.0 + N[(t$95$1 + N[(4.0 * N[(N[(b * b), $MachinePrecision] * 3.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \mathsf{fma}\left(a, a, b \cdot b\right)\\
t_1 := \frac{t\_0}{\frac{1}{t\_0}}\\
t_2 := 4 \cdot \left(\left(a \cdot a\right) \cdot \left(1 - a\right) + \left(b \cdot b\right) \cdot \left(a + 3\right)\right)\\
\mathbf{if}\;{\left(a \cdot a + b \cdot b\right)}^{2} + t\_2 \leq 2 \cdot 10^{+62}:\\
\;\;\;\;\left(t\_2 + t\_1\right) + -1\\
\mathbf{else}:\\
\;\;\;\;-1 + \left(t\_1 + 4 \cdot \left(\left(b \cdot b\right) \cdot 3\right)\right)\\
\end{array}
\end{array}
if (+.f64 (pow.f64 (+.f64 (*.f64 a a) (*.f64 b b)) #s(literal 2 binary64)) (*.f64 #s(literal 4 binary64) (+.f64 (*.f64 (*.f64 a a) (-.f64 #s(literal 1 binary64) a)) (*.f64 (*.f64 b b) (+.f64 #s(literal 3 binary64) a))))) < 2.00000000000000007e62Initial program 99.9%
lift-pow.f64N/A
unpow2N/A
lift-+.f64N/A
flip3-+N/A
clear-numN/A
un-div-invN/A
lower-/.f64N/A
lift-+.f64N/A
lift-*.f64N/A
lower-fma.f64N/A
clear-numN/A
Applied rewrites99.9%
if 2.00000000000000007e62 < (+.f64 (pow.f64 (+.f64 (*.f64 a a) (*.f64 b b)) #s(literal 2 binary64)) (*.f64 #s(literal 4 binary64) (+.f64 (*.f64 (*.f64 a a) (-.f64 #s(literal 1 binary64) a)) (*.f64 (*.f64 b b) (+.f64 #s(literal 3 binary64) a))))) Initial program 63.6%
lift-pow.f64N/A
unpow2N/A
lift-+.f64N/A
flip3-+N/A
clear-numN/A
un-div-invN/A
lower-/.f64N/A
lift-+.f64N/A
lift-*.f64N/A
lower-fma.f64N/A
clear-numN/A
Applied rewrites63.6%
Taylor expanded in a around 0
*-commutativeN/A
lower-*.f64N/A
unpow2N/A
lower-*.f6499.9
Applied rewrites99.9%
Final simplification99.9%
herbie shell --seed 2024228
(FPCore (a b)
:name "Bouland and Aaronson, Equation (24)"
:precision binary64
(- (+ (pow (+ (* a a) (* b b)) 2.0) (* 4.0 (+ (* (* a a) (- 1.0 a)) (* (* b b) (+ 3.0 a))))) 1.0))