
(FPCore (a b) :precision binary64 (- (+ (pow (+ (* a a) (* b b)) 2.0) (* 4.0 (+ (* (* a a) (+ 1.0 a)) (* (* b b) (- 1.0 (* 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) * (1.0 - (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) * (1.0d0 - (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) * (1.0 - (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) * (1.0 - (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(1.0 - 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) * (1.0 - (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[(1.0 - N[(3.0 * a), $MachinePrecision]), $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(1 - 3 \cdot a\right)\right)\right) - 1
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 9 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) (- 1.0 (* 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) * (1.0 - (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) * (1.0d0 - (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) * (1.0 - (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) * (1.0 - (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(1.0 - 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) * (1.0 - (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[(1.0 - N[(3.0 * a), $MachinePrecision]), $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(1 - 3 \cdot a\right)\right)\right) - 1
\end{array}
(FPCore (a b) :precision binary64 (fma (* (fma (+ 4.0 a) a 4.0) a) a (fma (* (fma (fma 2.0 a -12.0) a (fma b b 4.0)) b) b -1.0)))
double code(double a, double b) {
return fma((fma((4.0 + a), a, 4.0) * a), a, fma((fma(fma(2.0, a, -12.0), a, fma(b, b, 4.0)) * b), b, -1.0));
}
function code(a, b) return fma(Float64(fma(Float64(4.0 + a), a, 4.0) * a), a, fma(Float64(fma(fma(2.0, a, -12.0), a, fma(b, b, 4.0)) * b), b, -1.0)) end
code[a_, b_] := N[(N[(N[(N[(4.0 + a), $MachinePrecision] * a + 4.0), $MachinePrecision] * a), $MachinePrecision] * a + N[(N[(N[(N[(2.0 * a + -12.0), $MachinePrecision] * a + N[(b * b + 4.0), $MachinePrecision]), $MachinePrecision] * b), $MachinePrecision] * b + -1.0), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\mathsf{fma}\left(\mathsf{fma}\left(4 + a, a, 4\right) \cdot a, a, \mathsf{fma}\left(\mathsf{fma}\left(\mathsf{fma}\left(2, a, -12\right), a, \mathsf{fma}\left(b, b, 4\right)\right) \cdot b, b, -1\right)\right)
\end{array}
Initial program 69.4%
Taylor expanded in b around 0
+-commutativeN/A
associate-+l+N/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
+-commutativeN/A
lower-fma.f64N/A
Applied rewrites79.2%
Taylor expanded in b around 0
Applied rewrites99.9%
(FPCore (a b) :precision binary64 (fma (* (fma b b 4.0) b) b (+ -1.0 (* (* (fma (+ 4.0 a) a 4.0) a) a))))
double code(double a, double b) {
return fma((fma(b, b, 4.0) * b), b, (-1.0 + ((fma((4.0 + a), a, 4.0) * a) * a)));
}
function code(a, b) return fma(Float64(fma(b, b, 4.0) * b), b, Float64(-1.0 + Float64(Float64(fma(Float64(4.0 + a), a, 4.0) * a) * a))) end
code[a_, b_] := N[(N[(N[(b * b + 4.0), $MachinePrecision] * b), $MachinePrecision] * b + N[(-1.0 + N[(N[(N[(N[(4.0 + a), $MachinePrecision] * a + 4.0), $MachinePrecision] * a), $MachinePrecision] * a), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\mathsf{fma}\left(\mathsf{fma}\left(b, b, 4\right) \cdot b, b, -1 + \left(\mathsf{fma}\left(4 + a, a, 4\right) \cdot a\right) \cdot a\right)
\end{array}
Initial program 69.4%
Taylor expanded in b around 0
+-commutativeN/A
associate-+l+N/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
+-commutativeN/A
lower-fma.f64N/A
Applied rewrites79.2%
Taylor expanded in b around 0
Applied rewrites99.9%
Taylor expanded in a around 0
Applied rewrites99.7%
Applied rewrites99.7%
(FPCore (a b) :precision binary64 (fma (* (fma (+ 4.0 a) a 4.0) a) a (fma (* (fma b b 4.0) b) b -1.0)))
double code(double a, double b) {
return fma((fma((4.0 + a), a, 4.0) * a), a, fma((fma(b, b, 4.0) * b), b, -1.0));
}
function code(a, b) return fma(Float64(fma(Float64(4.0 + a), a, 4.0) * a), a, fma(Float64(fma(b, b, 4.0) * b), b, -1.0)) end
code[a_, b_] := N[(N[(N[(N[(4.0 + a), $MachinePrecision] * a + 4.0), $MachinePrecision] * a), $MachinePrecision] * a + N[(N[(N[(b * b + 4.0), $MachinePrecision] * b), $MachinePrecision] * b + -1.0), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\mathsf{fma}\left(\mathsf{fma}\left(4 + a, a, 4\right) \cdot a, a, \mathsf{fma}\left(\mathsf{fma}\left(b, b, 4\right) \cdot b, b, -1\right)\right)
\end{array}
Initial program 69.4%
Taylor expanded in b around 0
+-commutativeN/A
associate-+l+N/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
+-commutativeN/A
lower-fma.f64N/A
Applied rewrites79.2%
Taylor expanded in b around 0
Applied rewrites99.9%
Taylor expanded in a around 0
Applied rewrites99.7%
(FPCore (a b) :precision binary64 (if (or (<= a -4.4e+23) (not (<= a 9.5e+17))) (- (* (* a a) (* a a)) 1.0) (fma (* (fma b b 4.0) b) b -1.0)))
double code(double a, double b) {
double tmp;
if ((a <= -4.4e+23) || !(a <= 9.5e+17)) {
tmp = ((a * a) * (a * a)) - 1.0;
} else {
tmp = fma((fma(b, b, 4.0) * b), b, -1.0);
}
return tmp;
}
function code(a, b) tmp = 0.0 if ((a <= -4.4e+23) || !(a <= 9.5e+17)) tmp = Float64(Float64(Float64(a * a) * Float64(a * a)) - 1.0); else tmp = fma(Float64(fma(b, b, 4.0) * b), b, -1.0); end return tmp end
code[a_, b_] := If[Or[LessEqual[a, -4.4e+23], N[Not[LessEqual[a, 9.5e+17]], $MachinePrecision]], N[(N[(N[(a * a), $MachinePrecision] * N[(a * a), $MachinePrecision]), $MachinePrecision] - 1.0), $MachinePrecision], N[(N[(N[(b * b + 4.0), $MachinePrecision] * b), $MachinePrecision] * b + -1.0), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;a \leq -4.4 \cdot 10^{+23} \lor \neg \left(a \leq 9.5 \cdot 10^{+17}\right):\\
\;\;\;\;\left(a \cdot a\right) \cdot \left(a \cdot a\right) - 1\\
\mathbf{else}:\\
\;\;\;\;\mathsf{fma}\left(\mathsf{fma}\left(b, b, 4\right) \cdot b, b, -1\right)\\
\end{array}
\end{array}
if a < -4.40000000000000017e23 or 9.5e17 < a Initial program 34.9%
Taylor expanded in a around inf
lower-pow.f6497.2
Applied rewrites97.2%
Applied rewrites97.1%
if -4.40000000000000017e23 < a < 9.5e17Initial program 99.9%
Taylor expanded in b around 0
+-commutativeN/A
associate-+l+N/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
+-commutativeN/A
lower-fma.f64N/A
Applied rewrites99.9%
Taylor expanded in a around 0
sub-negN/A
metadata-evalN/A
pow-sqrN/A
distribute-rgt-inN/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
metadata-evalN/A
lower-fma.f64N/A
*-commutativeN/A
lower-*.f64N/A
+-commutativeN/A
unpow2N/A
lower-fma.f6498.8
Applied rewrites98.8%
Final simplification98.0%
(FPCore (a b)
:precision binary64
(if (<= a -4.4e+23)
(- (* (* (* a a) a) a) 1.0)
(if (<= a 9.5e+17)
(fma (* (fma b b 4.0) b) b -1.0)
(- (* (* a a) (* a a)) 1.0))))
double code(double a, double b) {
double tmp;
if (a <= -4.4e+23) {
tmp = (((a * a) * a) * a) - 1.0;
} else if (a <= 9.5e+17) {
tmp = fma((fma(b, b, 4.0) * b), b, -1.0);
} else {
tmp = ((a * a) * (a * a)) - 1.0;
}
return tmp;
}
function code(a, b) tmp = 0.0 if (a <= -4.4e+23) tmp = Float64(Float64(Float64(Float64(a * a) * a) * a) - 1.0); elseif (a <= 9.5e+17) tmp = fma(Float64(fma(b, b, 4.0) * b), b, -1.0); else tmp = Float64(Float64(Float64(a * a) * Float64(a * a)) - 1.0); end return tmp end
code[a_, b_] := If[LessEqual[a, -4.4e+23], N[(N[(N[(N[(a * a), $MachinePrecision] * a), $MachinePrecision] * a), $MachinePrecision] - 1.0), $MachinePrecision], If[LessEqual[a, 9.5e+17], N[(N[(N[(b * b + 4.0), $MachinePrecision] * b), $MachinePrecision] * b + -1.0), $MachinePrecision], N[(N[(N[(a * a), $MachinePrecision] * N[(a * a), $MachinePrecision]), $MachinePrecision] - 1.0), $MachinePrecision]]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;a \leq -4.4 \cdot 10^{+23}:\\
\;\;\;\;\left(\left(a \cdot a\right) \cdot a\right) \cdot a - 1\\
\mathbf{elif}\;a \leq 9.5 \cdot 10^{+17}:\\
\;\;\;\;\mathsf{fma}\left(\mathsf{fma}\left(b, b, 4\right) \cdot b, b, -1\right)\\
\mathbf{else}:\\
\;\;\;\;\left(a \cdot a\right) \cdot \left(a \cdot a\right) - 1\\
\end{array}
\end{array}
if a < -4.40000000000000017e23Initial program 20.8%
Taylor expanded in b around 0
+-commutativeN/A
associate-+l+N/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
+-commutativeN/A
lower-fma.f64N/A
Applied rewrites20.9%
Taylor expanded in b around 0
*-commutativeN/A
associate-*l*N/A
+-commutativeN/A
distribute-rgt1-inN/A
*-commutativeN/A
metadata-evalN/A
pow-sqrN/A
distribute-lft-inN/A
associate-+r+N/A
unpow2N/A
distribute-rgt-inN/A
*-commutativeN/A
unpow2N/A
associate-*r*N/A
lower-*.f64N/A
Applied rewrites97.6%
Taylor expanded in a around inf
Applied rewrites97.6%
if -4.40000000000000017e23 < a < 9.5e17Initial program 99.9%
Taylor expanded in b around 0
+-commutativeN/A
associate-+l+N/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
+-commutativeN/A
lower-fma.f64N/A
Applied rewrites99.9%
Taylor expanded in a around 0
sub-negN/A
metadata-evalN/A
pow-sqrN/A
distribute-rgt-inN/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
metadata-evalN/A
lower-fma.f64N/A
*-commutativeN/A
lower-*.f64N/A
+-commutativeN/A
unpow2N/A
lower-fma.f6498.8
Applied rewrites98.8%
if 9.5e17 < a Initial program 52.7%
Taylor expanded in a around inf
lower-pow.f6496.5
Applied rewrites96.5%
Applied rewrites96.4%
(FPCore (a b) :precision binary64 (if (or (<= a -8.5e+139) (not (<= a 6.6e+153))) (- (* (* a a) 4.0) 1.0) (fma (* (fma b b 4.0) b) b -1.0)))
double code(double a, double b) {
double tmp;
if ((a <= -8.5e+139) || !(a <= 6.6e+153)) {
tmp = ((a * a) * 4.0) - 1.0;
} else {
tmp = fma((fma(b, b, 4.0) * b), b, -1.0);
}
return tmp;
}
function code(a, b) tmp = 0.0 if ((a <= -8.5e+139) || !(a <= 6.6e+153)) tmp = Float64(Float64(Float64(a * a) * 4.0) - 1.0); else tmp = fma(Float64(fma(b, b, 4.0) * b), b, -1.0); end return tmp end
code[a_, b_] := If[Or[LessEqual[a, -8.5e+139], N[Not[LessEqual[a, 6.6e+153]], $MachinePrecision]], N[(N[(N[(a * a), $MachinePrecision] * 4.0), $MachinePrecision] - 1.0), $MachinePrecision], N[(N[(N[(b * b + 4.0), $MachinePrecision] * b), $MachinePrecision] * b + -1.0), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;a \leq -8.5 \cdot 10^{+139} \lor \neg \left(a \leq 6.6 \cdot 10^{+153}\right):\\
\;\;\;\;\left(a \cdot a\right) \cdot 4 - 1\\
\mathbf{else}:\\
\;\;\;\;\mathsf{fma}\left(\mathsf{fma}\left(b, b, 4\right) \cdot b, b, -1\right)\\
\end{array}
\end{array}
if a < -8.5e139 or 6.59999999999999989e153 < a Initial program 24.0%
Taylor expanded in b around 0
*-commutativeN/A
lower-fma.f64N/A
unpow2N/A
associate-*r*N/A
*-commutativeN/A
lower-*.f64N/A
+-commutativeN/A
distribute-lft-inN/A
*-rgt-identityN/A
lower-fma.f64N/A
lower-pow.f6445.3
Applied rewrites45.3%
Taylor expanded in a around 0
Applied rewrites96.4%
if -8.5e139 < a < 6.59999999999999989e153Initial program 88.2%
Taylor expanded in b around 0
+-commutativeN/A
associate-+l+N/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
+-commutativeN/A
lower-fma.f64N/A
Applied rewrites93.3%
Taylor expanded in a around 0
sub-negN/A
metadata-evalN/A
pow-sqrN/A
distribute-rgt-inN/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
metadata-evalN/A
lower-fma.f64N/A
*-commutativeN/A
lower-*.f64N/A
+-commutativeN/A
unpow2N/A
lower-fma.f6485.0
Applied rewrites85.0%
Final simplification88.3%
(FPCore (a b) :precision binary64 (if (<= (* b b) 5e+256) (- (* (* a a) 4.0) 1.0) (fma (* 4.0 b) b -1.0)))
double code(double a, double b) {
double tmp;
if ((b * b) <= 5e+256) {
tmp = ((a * a) * 4.0) - 1.0;
} else {
tmp = fma((4.0 * b), b, -1.0);
}
return tmp;
}
function code(a, b) tmp = 0.0 if (Float64(b * b) <= 5e+256) tmp = Float64(Float64(Float64(a * a) * 4.0) - 1.0); else tmp = fma(Float64(4.0 * b), b, -1.0); end return tmp end
code[a_, b_] := If[LessEqual[N[(b * b), $MachinePrecision], 5e+256], N[(N[(N[(a * a), $MachinePrecision] * 4.0), $MachinePrecision] - 1.0), $MachinePrecision], N[(N[(4.0 * b), $MachinePrecision] * b + -1.0), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;b \cdot b \leq 5 \cdot 10^{+256}:\\
\;\;\;\;\left(a \cdot a\right) \cdot 4 - 1\\
\mathbf{else}:\\
\;\;\;\;\mathsf{fma}\left(4 \cdot b, b, -1\right)\\
\end{array}
\end{array}
if (*.f64 b b) < 5.00000000000000015e256Initial program 73.9%
Taylor expanded in b around 0
*-commutativeN/A
lower-fma.f64N/A
unpow2N/A
associate-*r*N/A
*-commutativeN/A
lower-*.f64N/A
+-commutativeN/A
distribute-lft-inN/A
*-rgt-identityN/A
lower-fma.f64N/A
lower-pow.f6459.9
Applied rewrites59.9%
Taylor expanded in a around 0
Applied rewrites63.4%
if 5.00000000000000015e256 < (*.f64 b b) Initial program 57.7%
Taylor expanded in b around 0
+-commutativeN/A
associate-+l+N/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
+-commutativeN/A
lower-fma.f64N/A
Applied rewrites80.3%
Taylor expanded in a around 0
sub-negN/A
metadata-evalN/A
pow-sqrN/A
distribute-rgt-inN/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
metadata-evalN/A
lower-fma.f64N/A
*-commutativeN/A
lower-*.f64N/A
+-commutativeN/A
unpow2N/A
lower-fma.f64100.0
Applied rewrites100.0%
Taylor expanded in b around 0
Applied rewrites91.1%
(FPCore (a b) :precision binary64 (fma (* 4.0 b) b -1.0))
double code(double a, double b) {
return fma((4.0 * b), b, -1.0);
}
function code(a, b) return fma(Float64(4.0 * b), b, -1.0) end
code[a_, b_] := N[(N[(4.0 * b), $MachinePrecision] * b + -1.0), $MachinePrecision]
\begin{array}{l}
\\
\mathsf{fma}\left(4 \cdot b, b, -1\right)
\end{array}
Initial program 69.4%
Taylor expanded in b around 0
+-commutativeN/A
associate-+l+N/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
+-commutativeN/A
lower-fma.f64N/A
Applied rewrites79.2%
Taylor expanded in a around 0
sub-negN/A
metadata-evalN/A
pow-sqrN/A
distribute-rgt-inN/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
metadata-evalN/A
lower-fma.f64N/A
*-commutativeN/A
lower-*.f64N/A
+-commutativeN/A
unpow2N/A
lower-fma.f6471.0
Applied rewrites71.0%
Taylor expanded in b around 0
Applied rewrites50.9%
(FPCore (a b) :precision binary64 -1.0)
double code(double a, double b) {
return -1.0;
}
real(8) function code(a, b)
real(8), intent (in) :: a
real(8), intent (in) :: b
code = -1.0d0
end function
public static double code(double a, double b) {
return -1.0;
}
def code(a, b): return -1.0
function code(a, b) return -1.0 end
function tmp = code(a, b) tmp = -1.0; end
code[a_, b_] := -1.0
\begin{array}{l}
\\
-1
\end{array}
Initial program 69.4%
Taylor expanded in b around 0
+-commutativeN/A
associate-+l+N/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
+-commutativeN/A
lower-fma.f64N/A
Applied rewrites79.2%
Taylor expanded in a around 0
sub-negN/A
metadata-evalN/A
pow-sqrN/A
distribute-rgt-inN/A
unpow2N/A
associate-*l*N/A
*-commutativeN/A
metadata-evalN/A
lower-fma.f64N/A
*-commutativeN/A
lower-*.f64N/A
+-commutativeN/A
unpow2N/A
lower-fma.f6471.0
Applied rewrites71.0%
Taylor expanded in b around 0
Applied rewrites24.5%
herbie shell --seed 2024321
(FPCore (a b)
:name "Bouland and Aaronson, Equation (25)"
:precision binary64
(- (+ (pow (+ (* a a) (* b b)) 2.0) (* 4.0 (+ (* (* a a) (+ 1.0 a)) (* (* b b) (- 1.0 (* 3.0 a)))))) 1.0))