
(FPCore (a b c) :precision binary64 (/ (+ (- b) (sqrt (- (* b b) (* (* 3.0 a) c)))) (* 3.0 a)))
double code(double a, double b, double c) {
return (-b + sqrt(((b * b) - ((3.0 * a) * c)))) / (3.0 * a);
}
real(8) function code(a, b, c)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: c
code = (-b + sqrt(((b * b) - ((3.0d0 * a) * c)))) / (3.0d0 * a)
end function
public static double code(double a, double b, double c) {
return (-b + Math.sqrt(((b * b) - ((3.0 * a) * c)))) / (3.0 * a);
}
def code(a, b, c): return (-b + math.sqrt(((b * b) - ((3.0 * a) * c)))) / (3.0 * a)
function code(a, b, c) return Float64(Float64(Float64(-b) + sqrt(Float64(Float64(b * b) - Float64(Float64(3.0 * a) * c)))) / Float64(3.0 * a)) end
function tmp = code(a, b, c) tmp = (-b + sqrt(((b * b) - ((3.0 * a) * c)))) / (3.0 * a); end
code[a_, b_, c_] := N[(N[((-b) + N[Sqrt[N[(N[(b * b), $MachinePrecision] - N[(N[(3.0 * a), $MachinePrecision] * c), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] / N[(3.0 * a), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{\left(-b\right) + \sqrt{b \cdot b - \left(3 \cdot a\right) \cdot c}}{3 \cdot a}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 14 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (a b c) :precision binary64 (/ (+ (- b) (sqrt (- (* b b) (* (* 3.0 a) c)))) (* 3.0 a)))
double code(double a, double b, double c) {
return (-b + sqrt(((b * b) - ((3.0 * a) * c)))) / (3.0 * a);
}
real(8) function code(a, b, c)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: c
code = (-b + sqrt(((b * b) - ((3.0d0 * a) * c)))) / (3.0d0 * a)
end function
public static double code(double a, double b, double c) {
return (-b + Math.sqrt(((b * b) - ((3.0 * a) * c)))) / (3.0 * a);
}
def code(a, b, c): return (-b + math.sqrt(((b * b) - ((3.0 * a) * c)))) / (3.0 * a)
function code(a, b, c) return Float64(Float64(Float64(-b) + sqrt(Float64(Float64(b * b) - Float64(Float64(3.0 * a) * c)))) / Float64(3.0 * a)) end
function tmp = code(a, b, c) tmp = (-b + sqrt(((b * b) - ((3.0 * a) * c)))) / (3.0 * a); end
code[a_, b_, c_] := N[(N[((-b) + N[Sqrt[N[(N[(b * b), $MachinePrecision] - N[(N[(3.0 * a), $MachinePrecision] * c), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] / N[(3.0 * a), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{\left(-b\right) + \sqrt{b \cdot b - \left(3 \cdot a\right) \cdot c}}{3 \cdot a}
\end{array}
(FPCore (a b c)
:precision binary64
(let* ((t_0 (* b (* b b)))
(t_1 (fma c (* a -3.0) (* b b)))
(t_2 (sqrt t_1))
(t_3 (* t_1 t_2)))
(if (<= b 0.118)
(/
(/
(fma t_0 (* t_0 t_0) (* (* t_1 (* t_1 t_1)) (- t_3)))
(fma b (* (* b b) t_0) (+ (* t_3 t_3) (* t_0 t_3))))
(* (* a -3.0) (fma b b (fma b t_2 t_1))))
(/
1.0
(fma
a
(fma
a
(*
-3.0
(fma a (/ (* (* c c) -0.5625) (pow b 5.0)) (* -0.375 (/ c t_0))))
(/ 1.5 b))
(* -2.0 (/ b c)))))))
double code(double a, double b, double c) {
double t_0 = b * (b * b);
double t_1 = fma(c, (a * -3.0), (b * b));
double t_2 = sqrt(t_1);
double t_3 = t_1 * t_2;
double tmp;
if (b <= 0.118) {
tmp = (fma(t_0, (t_0 * t_0), ((t_1 * (t_1 * t_1)) * -t_3)) / fma(b, ((b * b) * t_0), ((t_3 * t_3) + (t_0 * t_3)))) / ((a * -3.0) * fma(b, b, fma(b, t_2, t_1)));
} else {
tmp = 1.0 / fma(a, fma(a, (-3.0 * fma(a, (((c * c) * -0.5625) / pow(b, 5.0)), (-0.375 * (c / t_0)))), (1.5 / b)), (-2.0 * (b / c)));
}
return tmp;
}
function code(a, b, c) t_0 = Float64(b * Float64(b * b)) t_1 = fma(c, Float64(a * -3.0), Float64(b * b)) t_2 = sqrt(t_1) t_3 = Float64(t_1 * t_2) tmp = 0.0 if (b <= 0.118) tmp = Float64(Float64(fma(t_0, Float64(t_0 * t_0), Float64(Float64(t_1 * Float64(t_1 * t_1)) * Float64(-t_3))) / fma(b, Float64(Float64(b * b) * t_0), Float64(Float64(t_3 * t_3) + Float64(t_0 * t_3)))) / Float64(Float64(a * -3.0) * fma(b, b, fma(b, t_2, t_1)))); else tmp = Float64(1.0 / fma(a, fma(a, Float64(-3.0 * fma(a, Float64(Float64(Float64(c * c) * -0.5625) / (b ^ 5.0)), Float64(-0.375 * Float64(c / t_0)))), Float64(1.5 / b)), Float64(-2.0 * Float64(b / c)))); end return tmp end
code[a_, b_, c_] := Block[{t$95$0 = N[(b * N[(b * b), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(c * N[(a * -3.0), $MachinePrecision] + N[(b * b), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$2 = N[Sqrt[t$95$1], $MachinePrecision]}, Block[{t$95$3 = N[(t$95$1 * t$95$2), $MachinePrecision]}, If[LessEqual[b, 0.118], N[(N[(N[(t$95$0 * N[(t$95$0 * t$95$0), $MachinePrecision] + N[(N[(t$95$1 * N[(t$95$1 * t$95$1), $MachinePrecision]), $MachinePrecision] * (-t$95$3)), $MachinePrecision]), $MachinePrecision] / N[(b * N[(N[(b * b), $MachinePrecision] * t$95$0), $MachinePrecision] + N[(N[(t$95$3 * t$95$3), $MachinePrecision] + N[(t$95$0 * t$95$3), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(N[(a * -3.0), $MachinePrecision] * N[(b * b + N[(b * t$95$2 + t$95$1), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(1.0 / N[(a * N[(a * N[(-3.0 * N[(a * N[(N[(N[(c * c), $MachinePrecision] * -0.5625), $MachinePrecision] / N[Power[b, 5.0], $MachinePrecision]), $MachinePrecision] + N[(-0.375 * N[(c / t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + N[(1.5 / b), $MachinePrecision]), $MachinePrecision] + N[(-2.0 * N[(b / c), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := b \cdot \left(b \cdot b\right)\\
t_1 := \mathsf{fma}\left(c, a \cdot -3, b \cdot b\right)\\
t_2 := \sqrt{t\_1}\\
t_3 := t\_1 \cdot t\_2\\
\mathbf{if}\;b \leq 0.118:\\
\;\;\;\;\frac{\frac{\mathsf{fma}\left(t\_0, t\_0 \cdot t\_0, \left(t\_1 \cdot \left(t\_1 \cdot t\_1\right)\right) \cdot \left(-t\_3\right)\right)}{\mathsf{fma}\left(b, \left(b \cdot b\right) \cdot t\_0, t\_3 \cdot t\_3 + t\_0 \cdot t\_3\right)}}{\left(a \cdot -3\right) \cdot \mathsf{fma}\left(b, b, \mathsf{fma}\left(b, t\_2, t\_1\right)\right)}\\
\mathbf{else}:\\
\;\;\;\;\frac{1}{\mathsf{fma}\left(a, \mathsf{fma}\left(a, -3 \cdot \mathsf{fma}\left(a, \frac{\left(c \cdot c\right) \cdot -0.5625}{{b}^{5}}, -0.375 \cdot \frac{c}{t\_0}\right), \frac{1.5}{b}\right), -2 \cdot \frac{b}{c}\right)}\\
\end{array}
\end{array}
if b < 0.11799999999999999Initial program 87.3%
Applied rewrites87.4%
Applied rewrites88.7%
Applied rewrites90.0%
if 0.11799999999999999 < b Initial program 52.5%
Applied rewrites52.5%
Applied rewrites52.5%
Taylor expanded in a around 0
Applied rewrites93.2%
Taylor expanded in c around 0
Applied rewrites93.2%
Final simplification92.9%
(FPCore (a b c)
:precision binary64
(let* ((t_0 (fma c (* a -3.0) (* b b)))
(t_1 (* (* c a) -9.0))
(t_2 (fma (* (* a a) (* c c)) 27.0 (* -0.25 (* t_1 t_1))))
(t_3 (fma (* -27.0 (* a (* a a))) (* c (* c c)) (* -0.5 (* t_1 t_2)))))
(/
(*
b
(-
(/ (* 0.5 (fma 0.25 (* t_2 t_2) (* (* 0.5 t_1) t_3))) (pow b 6.0))
(fma 0.5 t_1 (* 0.5 (+ (/ t_3 (pow b 4.0)) (/ t_2 (* b b)))))))
(* (* a -3.0) (fma b b (fma b (sqrt t_0) t_0))))))
double code(double a, double b, double c) {
double t_0 = fma(c, (a * -3.0), (b * b));
double t_1 = (c * a) * -9.0;
double t_2 = fma(((a * a) * (c * c)), 27.0, (-0.25 * (t_1 * t_1)));
double t_3 = fma((-27.0 * (a * (a * a))), (c * (c * c)), (-0.5 * (t_1 * t_2)));
return (b * (((0.5 * fma(0.25, (t_2 * t_2), ((0.5 * t_1) * t_3))) / pow(b, 6.0)) - fma(0.5, t_1, (0.5 * ((t_3 / pow(b, 4.0)) + (t_2 / (b * b))))))) / ((a * -3.0) * fma(b, b, fma(b, sqrt(t_0), t_0)));
}
function code(a, b, c) t_0 = fma(c, Float64(a * -3.0), Float64(b * b)) t_1 = Float64(Float64(c * a) * -9.0) t_2 = fma(Float64(Float64(a * a) * Float64(c * c)), 27.0, Float64(-0.25 * Float64(t_1 * t_1))) t_3 = fma(Float64(-27.0 * Float64(a * Float64(a * a))), Float64(c * Float64(c * c)), Float64(-0.5 * Float64(t_1 * t_2))) return Float64(Float64(b * Float64(Float64(Float64(0.5 * fma(0.25, Float64(t_2 * t_2), Float64(Float64(0.5 * t_1) * t_3))) / (b ^ 6.0)) - fma(0.5, t_1, Float64(0.5 * Float64(Float64(t_3 / (b ^ 4.0)) + Float64(t_2 / Float64(b * b))))))) / Float64(Float64(a * -3.0) * fma(b, b, fma(b, sqrt(t_0), t_0)))) end
code[a_, b_, c_] := Block[{t$95$0 = N[(c * N[(a * -3.0), $MachinePrecision] + N[(b * b), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(N[(c * a), $MachinePrecision] * -9.0), $MachinePrecision]}, Block[{t$95$2 = N[(N[(N[(a * a), $MachinePrecision] * N[(c * c), $MachinePrecision]), $MachinePrecision] * 27.0 + N[(-0.25 * N[(t$95$1 * t$95$1), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$3 = N[(N[(-27.0 * N[(a * N[(a * a), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(c * N[(c * c), $MachinePrecision]), $MachinePrecision] + N[(-0.5 * N[(t$95$1 * t$95$2), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, N[(N[(b * N[(N[(N[(0.5 * N[(0.25 * N[(t$95$2 * t$95$2), $MachinePrecision] + N[(N[(0.5 * t$95$1), $MachinePrecision] * t$95$3), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[Power[b, 6.0], $MachinePrecision]), $MachinePrecision] - N[(0.5 * t$95$1 + N[(0.5 * N[(N[(t$95$3 / N[Power[b, 4.0], $MachinePrecision]), $MachinePrecision] + N[(t$95$2 / N[(b * b), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(N[(a * -3.0), $MachinePrecision] * N[(b * b + N[(b * N[Sqrt[t$95$0], $MachinePrecision] + t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \mathsf{fma}\left(c, a \cdot -3, b \cdot b\right)\\
t_1 := \left(c \cdot a\right) \cdot -9\\
t_2 := \mathsf{fma}\left(\left(a \cdot a\right) \cdot \left(c \cdot c\right), 27, -0.25 \cdot \left(t\_1 \cdot t\_1\right)\right)\\
t_3 := \mathsf{fma}\left(-27 \cdot \left(a \cdot \left(a \cdot a\right)\right), c \cdot \left(c \cdot c\right), -0.5 \cdot \left(t\_1 \cdot t\_2\right)\right)\\
\frac{b \cdot \left(\frac{0.5 \cdot \mathsf{fma}\left(0.25, t\_2 \cdot t\_2, \left(0.5 \cdot t\_1\right) \cdot t\_3\right)}{{b}^{6}} - \mathsf{fma}\left(0.5, t\_1, 0.5 \cdot \left(\frac{t\_3}{{b}^{4}} + \frac{t\_2}{b \cdot b}\right)\right)\right)}{\left(a \cdot -3\right) \cdot \mathsf{fma}\left(b, b, \mathsf{fma}\left(b, \sqrt{t\_0}, t\_0\right)\right)}
\end{array}
\end{array}
Initial program 54.9%
Applied rewrites54.9%
Applied rewrites55.8%
Taylor expanded in b around inf
Applied rewrites91.5%
Final simplification91.5%
herbie shell --seed 2024223
(FPCore (a b c)
:name "Cubic critical, narrow range"
:precision binary64
:pre (and (and (and (< 1.0536712127723509e-8 a) (< a 94906265.62425156)) (and (< 1.0536712127723509e-8 b) (< b 94906265.62425156))) (and (< 1.0536712127723509e-8 c) (< c 94906265.62425156)))
(/ (+ (- b) (sqrt (- (* b b) (* (* 3.0 a) c)))) (* 3.0 a)))