
(FPCore (A B C F)
:precision binary64
(let* ((t_0 (- (pow B 2.0) (* (* 4.0 A) C))))
(/
(-
(sqrt
(*
(* 2.0 (* t_0 F))
(+ (+ A C) (sqrt (+ (pow (- A C) 2.0) (pow B 2.0)))))))
t_0)))
double code(double A, double B, double C, double F) {
double t_0 = pow(B, 2.0) - ((4.0 * A) * C);
return -sqrt(((2.0 * (t_0 * F)) * ((A + C) + sqrt((pow((A - C), 2.0) + pow(B, 2.0)))))) / t_0;
}
real(8) function code(a, b, c, f)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: c
real(8), intent (in) :: f
real(8) :: t_0
t_0 = (b ** 2.0d0) - ((4.0d0 * a) * c)
code = -sqrt(((2.0d0 * (t_0 * f)) * ((a + c) + sqrt((((a - c) ** 2.0d0) + (b ** 2.0d0)))))) / t_0
end function
public static double code(double A, double B, double C, double F) {
double t_0 = Math.pow(B, 2.0) - ((4.0 * A) * C);
return -Math.sqrt(((2.0 * (t_0 * F)) * ((A + C) + Math.sqrt((Math.pow((A - C), 2.0) + Math.pow(B, 2.0)))))) / t_0;
}
def code(A, B, C, F): t_0 = math.pow(B, 2.0) - ((4.0 * A) * C) return -math.sqrt(((2.0 * (t_0 * F)) * ((A + C) + math.sqrt((math.pow((A - C), 2.0) + math.pow(B, 2.0)))))) / t_0
function code(A, B, C, F) t_0 = Float64((B ^ 2.0) - Float64(Float64(4.0 * A) * C)) return Float64(Float64(-sqrt(Float64(Float64(2.0 * Float64(t_0 * F)) * Float64(Float64(A + C) + sqrt(Float64((Float64(A - C) ^ 2.0) + (B ^ 2.0))))))) / t_0) end
function tmp = code(A, B, C, F) t_0 = (B ^ 2.0) - ((4.0 * A) * C); tmp = -sqrt(((2.0 * (t_0 * F)) * ((A + C) + sqrt((((A - C) ^ 2.0) + (B ^ 2.0)))))) / t_0; end
code[A_, B_, C_, F_] := Block[{t$95$0 = N[(N[Power[B, 2.0], $MachinePrecision] - N[(N[(4.0 * A), $MachinePrecision] * C), $MachinePrecision]), $MachinePrecision]}, N[((-N[Sqrt[N[(N[(2.0 * N[(t$95$0 * F), $MachinePrecision]), $MachinePrecision] * N[(N[(A + C), $MachinePrecision] + N[Sqrt[N[(N[Power[N[(A - C), $MachinePrecision], 2.0], $MachinePrecision] + N[Power[B, 2.0], $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]) / t$95$0), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := {B}^{2} - \left(4 \cdot A\right) \cdot C\\
\frac{-\sqrt{\left(2 \cdot \left(t\_0 \cdot F\right)\right) \cdot \left(\left(A + C\right) + \sqrt{{\left(A - C\right)}^{2} + {B}^{2}}\right)}}{t\_0}
\end{array}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 9 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (A B C F)
:precision binary64
(let* ((t_0 (- (pow B 2.0) (* (* 4.0 A) C))))
(/
(-
(sqrt
(*
(* 2.0 (* t_0 F))
(+ (+ A C) (sqrt (+ (pow (- A C) 2.0) (pow B 2.0)))))))
t_0)))
double code(double A, double B, double C, double F) {
double t_0 = pow(B, 2.0) - ((4.0 * A) * C);
return -sqrt(((2.0 * (t_0 * F)) * ((A + C) + sqrt((pow((A - C), 2.0) + pow(B, 2.0)))))) / t_0;
}
real(8) function code(a, b, c, f)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: c
real(8), intent (in) :: f
real(8) :: t_0
t_0 = (b ** 2.0d0) - ((4.0d0 * a) * c)
code = -sqrt(((2.0d0 * (t_0 * f)) * ((a + c) + sqrt((((a - c) ** 2.0d0) + (b ** 2.0d0)))))) / t_0
end function
public static double code(double A, double B, double C, double F) {
double t_0 = Math.pow(B, 2.0) - ((4.0 * A) * C);
return -Math.sqrt(((2.0 * (t_0 * F)) * ((A + C) + Math.sqrt((Math.pow((A - C), 2.0) + Math.pow(B, 2.0)))))) / t_0;
}
def code(A, B, C, F): t_0 = math.pow(B, 2.0) - ((4.0 * A) * C) return -math.sqrt(((2.0 * (t_0 * F)) * ((A + C) + math.sqrt((math.pow((A - C), 2.0) + math.pow(B, 2.0)))))) / t_0
function code(A, B, C, F) t_0 = Float64((B ^ 2.0) - Float64(Float64(4.0 * A) * C)) return Float64(Float64(-sqrt(Float64(Float64(2.0 * Float64(t_0 * F)) * Float64(Float64(A + C) + sqrt(Float64((Float64(A - C) ^ 2.0) + (B ^ 2.0))))))) / t_0) end
function tmp = code(A, B, C, F) t_0 = (B ^ 2.0) - ((4.0 * A) * C); tmp = -sqrt(((2.0 * (t_0 * F)) * ((A + C) + sqrt((((A - C) ^ 2.0) + (B ^ 2.0)))))) / t_0; end
code[A_, B_, C_, F_] := Block[{t$95$0 = N[(N[Power[B, 2.0], $MachinePrecision] - N[(N[(4.0 * A), $MachinePrecision] * C), $MachinePrecision]), $MachinePrecision]}, N[((-N[Sqrt[N[(N[(2.0 * N[(t$95$0 * F), $MachinePrecision]), $MachinePrecision] * N[(N[(A + C), $MachinePrecision] + N[Sqrt[N[(N[Power[N[(A - C), $MachinePrecision], 2.0], $MachinePrecision] + N[Power[B, 2.0], $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]) / t$95$0), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := {B}^{2} - \left(4 \cdot A\right) \cdot C\\
\frac{-\sqrt{\left(2 \cdot \left(t\_0 \cdot F\right)\right) \cdot \left(\left(A + C\right) + \sqrt{{\left(A - C\right)}^{2} + {B}^{2}}\right)}}{t\_0}
\end{array}
\end{array}
B_m = (fabs.f64 B)
NOTE: A, B_m, C, and F should be sorted in increasing order before calling this function.
(FPCore (A B_m C F)
:precision binary64
(let* ((t_0 (fma B_m B_m (* -4.0 (* A C)))))
(if (<= (pow B_m 2.0) 5e-257)
(/ (sqrt (* (* F (* 2.0 t_0)) (* 2.0 C))) (- t_0))
(if (<= (pow B_m 2.0) 5e+104)
(*
(- (sqrt F))
(* (sqrt (/ (* (* B_m B_m) -0.5) A)) (/ (sqrt 2.0) B_m)))
(- (/ (sqrt (* 2.0 F)) (sqrt B_m)))))))B_m = fabs(B);
assert(A < B_m && B_m < C && C < F);
double code(double A, double B_m, double C, double F) {
double t_0 = fma(B_m, B_m, (-4.0 * (A * C)));
double tmp;
if (pow(B_m, 2.0) <= 5e-257) {
tmp = sqrt(((F * (2.0 * t_0)) * (2.0 * C))) / -t_0;
} else if (pow(B_m, 2.0) <= 5e+104) {
tmp = -sqrt(F) * (sqrt((((B_m * B_m) * -0.5) / A)) * (sqrt(2.0) / B_m));
} else {
tmp = -(sqrt((2.0 * F)) / sqrt(B_m));
}
return tmp;
}
B_m = abs(B) A, B_m, C, F = sort([A, B_m, C, F]) function code(A, B_m, C, F) t_0 = fma(B_m, B_m, Float64(-4.0 * Float64(A * C))) tmp = 0.0 if ((B_m ^ 2.0) <= 5e-257) tmp = Float64(sqrt(Float64(Float64(F * Float64(2.0 * t_0)) * Float64(2.0 * C))) / Float64(-t_0)); elseif ((B_m ^ 2.0) <= 5e+104) tmp = Float64(Float64(-sqrt(F)) * Float64(sqrt(Float64(Float64(Float64(B_m * B_m) * -0.5) / A)) * Float64(sqrt(2.0) / B_m))); else tmp = Float64(-Float64(sqrt(Float64(2.0 * F)) / sqrt(B_m))); end return tmp end
B_m = N[Abs[B], $MachinePrecision]
NOTE: A, B_m, C, and F should be sorted in increasing order before calling this function.
code[A_, B$95$m_, C_, F_] := Block[{t$95$0 = N[(B$95$m * B$95$m + N[(-4.0 * N[(A * C), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[N[Power[B$95$m, 2.0], $MachinePrecision], 5e-257], N[(N[Sqrt[N[(N[(F * N[(2.0 * t$95$0), $MachinePrecision]), $MachinePrecision] * N[(2.0 * C), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] / (-t$95$0)), $MachinePrecision], If[LessEqual[N[Power[B$95$m, 2.0], $MachinePrecision], 5e+104], N[((-N[Sqrt[F], $MachinePrecision]) * N[(N[Sqrt[N[(N[(N[(B$95$m * B$95$m), $MachinePrecision] * -0.5), $MachinePrecision] / A), $MachinePrecision]], $MachinePrecision] * N[(N[Sqrt[2.0], $MachinePrecision] / B$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], (-N[(N[Sqrt[N[(2.0 * F), $MachinePrecision]], $MachinePrecision] / N[Sqrt[B$95$m], $MachinePrecision]), $MachinePrecision])]]]
\begin{array}{l}
B_m = \left|B\right|
\\
[A, B_m, C, F] = \mathsf{sort}([A, B_m, C, F])\\
\\
\begin{array}{l}
t_0 := \mathsf{fma}\left(B\_m, B\_m, -4 \cdot \left(A \cdot C\right)\right)\\
\mathbf{if}\;{B\_m}^{2} \leq 5 \cdot 10^{-257}:\\
\;\;\;\;\frac{\sqrt{\left(F \cdot \left(2 \cdot t\_0\right)\right) \cdot \left(2 \cdot C\right)}}{-t\_0}\\
\mathbf{elif}\;{B\_m}^{2} \leq 5 \cdot 10^{+104}:\\
\;\;\;\;\left(-\sqrt{F}\right) \cdot \left(\sqrt{\frac{\left(B\_m \cdot B\_m\right) \cdot -0.5}{A}} \cdot \frac{\sqrt{2}}{B\_m}\right)\\
\mathbf{else}:\\
\;\;\;\;-\frac{\sqrt{2 \cdot F}}{\sqrt{B\_m}}\\
\end{array}
\end{array}
if (pow.f64 B #s(literal 2 binary64)) < 4.99999999999999989e-257Initial program 24.5%
Taylor expanded in A around 0
lower-+.f64N/A
lower-sqrt.f64N/A
+-commutativeN/A
unpow2N/A
lower-fma.f64N/A
unpow2N/A
lower-*.f6430.3
Applied rewrites30.3%
Applied rewrites30.3%
Taylor expanded in A around -inf
lower-*.f6454.5
Applied rewrites54.5%
if 4.99999999999999989e-257 < (pow.f64 B #s(literal 2 binary64)) < 4.9999999999999997e104Initial program 25.3%
Taylor expanded in C around 0
mul-1-negN/A
*-commutativeN/A
distribute-rgt-neg-inN/A
mul-1-negN/A
lower-*.f64N/A
Applied rewrites15.3%
Taylor expanded in A around -inf
Applied rewrites30.3%
Applied rewrites39.1%
if 4.9999999999999997e104 < (pow.f64 B #s(literal 2 binary64)) Initial program 12.0%
Taylor expanded in B around inf
mul-1-negN/A
lower-neg.f64N/A
*-commutativeN/A
lower-*.f64N/A
lower-sqrt.f64N/A
lower-sqrt.f64N/A
lower-/.f6449.4
Applied rewrites49.4%
Applied rewrites68.9%
Final simplification55.9%
B_m = (fabs.f64 B)
NOTE: A, B_m, C, and F should be sorted in increasing order before calling this function.
(FPCore (A B_m C F)
:precision binary64
(let* ((t_0 (* (* 4.0 A) C)))
(if (<= (pow B_m 2.0) 5e+104)
(/
(sqrt
(*
(* 2.0 (* (- (pow B_m 2.0) t_0) F))
(fma (/ (* B_m B_m) A) -0.5 (* 2.0 C))))
(- t_0 (pow B_m 2.0)))
(- (/ (sqrt (* 2.0 F)) (sqrt B_m))))))B_m = fabs(B);
assert(A < B_m && B_m < C && C < F);
double code(double A, double B_m, double C, double F) {
double t_0 = (4.0 * A) * C;
double tmp;
if (pow(B_m, 2.0) <= 5e+104) {
tmp = sqrt(((2.0 * ((pow(B_m, 2.0) - t_0) * F)) * fma(((B_m * B_m) / A), -0.5, (2.0 * C)))) / (t_0 - pow(B_m, 2.0));
} else {
tmp = -(sqrt((2.0 * F)) / sqrt(B_m));
}
return tmp;
}
B_m = abs(B) A, B_m, C, F = sort([A, B_m, C, F]) function code(A, B_m, C, F) t_0 = Float64(Float64(4.0 * A) * C) tmp = 0.0 if ((B_m ^ 2.0) <= 5e+104) tmp = Float64(sqrt(Float64(Float64(2.0 * Float64(Float64((B_m ^ 2.0) - t_0) * F)) * fma(Float64(Float64(B_m * B_m) / A), -0.5, Float64(2.0 * C)))) / Float64(t_0 - (B_m ^ 2.0))); else tmp = Float64(-Float64(sqrt(Float64(2.0 * F)) / sqrt(B_m))); end return tmp end
B_m = N[Abs[B], $MachinePrecision]
NOTE: A, B_m, C, and F should be sorted in increasing order before calling this function.
code[A_, B$95$m_, C_, F_] := Block[{t$95$0 = N[(N[(4.0 * A), $MachinePrecision] * C), $MachinePrecision]}, If[LessEqual[N[Power[B$95$m, 2.0], $MachinePrecision], 5e+104], N[(N[Sqrt[N[(N[(2.0 * N[(N[(N[Power[B$95$m, 2.0], $MachinePrecision] - t$95$0), $MachinePrecision] * F), $MachinePrecision]), $MachinePrecision] * N[(N[(N[(B$95$m * B$95$m), $MachinePrecision] / A), $MachinePrecision] * -0.5 + N[(2.0 * C), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] / N[(t$95$0 - N[Power[B$95$m, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], (-N[(N[Sqrt[N[(2.0 * F), $MachinePrecision]], $MachinePrecision] / N[Sqrt[B$95$m], $MachinePrecision]), $MachinePrecision])]]
\begin{array}{l}
B_m = \left|B\right|
\\
[A, B_m, C, F] = \mathsf{sort}([A, B_m, C, F])\\
\\
\begin{array}{l}
t_0 := \left(4 \cdot A\right) \cdot C\\
\mathbf{if}\;{B\_m}^{2} \leq 5 \cdot 10^{+104}:\\
\;\;\;\;\frac{\sqrt{\left(2 \cdot \left(\left({B\_m}^{2} - t\_0\right) \cdot F\right)\right) \cdot \mathsf{fma}\left(\frac{B\_m \cdot B\_m}{A}, -0.5, 2 \cdot C\right)}}{t\_0 - {B\_m}^{2}}\\
\mathbf{else}:\\
\;\;\;\;-\frac{\sqrt{2 \cdot F}}{\sqrt{B\_m}}\\
\end{array}
\end{array}
if (pow.f64 B #s(literal 2 binary64)) < 4.9999999999999997e104Initial program 24.7%
Taylor expanded in A around -inf
*-commutativeN/A
lower-fma.f64N/A
lower-/.f64N/A
unpow2N/A
lower-*.f64N/A
lower-*.f6441.3
Applied rewrites41.3%
if 4.9999999999999997e104 < (pow.f64 B #s(literal 2 binary64)) Initial program 9.8%
Taylor expanded in B around inf
mul-1-negN/A
lower-neg.f64N/A
*-commutativeN/A
lower-*.f64N/A
lower-sqrt.f64N/A
lower-sqrt.f64N/A
lower-/.f6445.8
Applied rewrites45.8%
Applied rewrites64.2%
Final simplification50.8%
herbie shell --seed 2024229
(FPCore (A B C F)
:name "ABCF->ab-angle a"
:precision binary64
(/ (- (sqrt (* (* 2.0 (* (- (pow B 2.0) (* (* 4.0 A) C)) F)) (+ (+ A C) (sqrt (+ (pow (- A C) 2.0) (pow B 2.0))))))) (- (pow B 2.0) (* (* 4.0 A) C))))