
(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 7 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 (/ (* c (* a 3.0)) (* (* -3.0 a) (+ b (sqrt (fma b b (* -3.0 (* a c))))))))
double code(double a, double b, double c) {
return (c * (a * 3.0)) / ((-3.0 * a) * (b + sqrt(fma(b, b, (-3.0 * (a * c))))));
}
function code(a, b, c) return Float64(Float64(c * Float64(a * 3.0)) / Float64(Float64(-3.0 * a) * Float64(b + sqrt(fma(b, b, Float64(-3.0 * Float64(a * c))))))) end
code[a_, b_, c_] := N[(N[(c * N[(a * 3.0), $MachinePrecision]), $MachinePrecision] / N[(N[(-3.0 * a), $MachinePrecision] * N[(b + N[Sqrt[N[(b * b + N[(-3.0 * N[(a * c), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{c \cdot \left(a \cdot 3\right)}{\left(-3 \cdot a\right) \cdot \left(b + \sqrt{\mathsf{fma}\left(b, b, -3 \cdot \left(a \cdot c\right)\right)}\right)}
\end{array}
Initial program 55.9%
Applied rewrites55.9%
Applied rewrites57.3%
Taylor expanded in b around 0
lower-*.f64N/A
lower-*.f6499.1
Applied rewrites99.1%
lift-/.f64N/A
lift-/.f64N/A
associate-/l/N/A
lower-/.f64N/A
Applied rewrites99.3%
(FPCore (a b c) :precision binary64 (* (/ (* c (* a 3.0)) (* a (+ b (sqrt (fma b b (* -3.0 (* a c))))))) -0.3333333333333333))
double code(double a, double b, double c) {
return ((c * (a * 3.0)) / (a * (b + sqrt(fma(b, b, (-3.0 * (a * c))))))) * -0.3333333333333333;
}
function code(a, b, c) return Float64(Float64(Float64(c * Float64(a * 3.0)) / Float64(a * Float64(b + sqrt(fma(b, b, Float64(-3.0 * Float64(a * c))))))) * -0.3333333333333333) end
code[a_, b_, c_] := N[(N[(N[(c * N[(a * 3.0), $MachinePrecision]), $MachinePrecision] / N[(a * N[(b + N[Sqrt[N[(b * b + N[(-3.0 * N[(a * c), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * -0.3333333333333333), $MachinePrecision]
\begin{array}{l}
\\
\frac{c \cdot \left(a \cdot 3\right)}{a \cdot \left(b + \sqrt{\mathsf{fma}\left(b, b, -3 \cdot \left(a \cdot c\right)\right)}\right)} \cdot -0.3333333333333333
\end{array}
Initial program 55.9%
Applied rewrites55.9%
Applied rewrites57.3%
Taylor expanded in b around 0
lower-*.f64N/A
lower-*.f6499.1
Applied rewrites99.1%
lift-/.f64N/A
div-invN/A
lower-*.f64N/A
Applied rewrites99.0%
(FPCore (a b c) :precision binary64 (if (<= b 6.7) (* (/ 0.3333333333333333 a) (- (sqrt (fma b b (* -3.0 (* a c)))) b)) (/ (fma a (/ (* (* c c) -0.375) (* b b)) (* c -0.5)) b)))
double code(double a, double b, double c) {
double tmp;
if (b <= 6.7) {
tmp = (0.3333333333333333 / a) * (sqrt(fma(b, b, (-3.0 * (a * c)))) - b);
} else {
tmp = fma(a, (((c * c) * -0.375) / (b * b)), (c * -0.5)) / b;
}
return tmp;
}
function code(a, b, c) tmp = 0.0 if (b <= 6.7) tmp = Float64(Float64(0.3333333333333333 / a) * Float64(sqrt(fma(b, b, Float64(-3.0 * Float64(a * c)))) - b)); else tmp = Float64(fma(a, Float64(Float64(Float64(c * c) * -0.375) / Float64(b * b)), Float64(c * -0.5)) / b); end return tmp end
code[a_, b_, c_] := If[LessEqual[b, 6.7], N[(N[(0.3333333333333333 / a), $MachinePrecision] * N[(N[Sqrt[N[(b * b + N[(-3.0 * N[(a * c), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] - b), $MachinePrecision]), $MachinePrecision], N[(N[(a * N[(N[(N[(c * c), $MachinePrecision] * -0.375), $MachinePrecision] / N[(b * b), $MachinePrecision]), $MachinePrecision] + N[(c * -0.5), $MachinePrecision]), $MachinePrecision] / b), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;b \leq 6.7:\\
\;\;\;\;\frac{0.3333333333333333}{a} \cdot \left(\sqrt{\mathsf{fma}\left(b, b, -3 \cdot \left(a \cdot c\right)\right)} - b\right)\\
\mathbf{else}:\\
\;\;\;\;\frac{\mathsf{fma}\left(a, \frac{\left(c \cdot c\right) \cdot -0.375}{b \cdot b}, c \cdot -0.5\right)}{b}\\
\end{array}
\end{array}
if b < 6.70000000000000018Initial program 80.9%
Applied rewrites80.7%
Taylor expanded in a around 0
+-commutativeN/A
unpow2N/A
lower-fma.f64N/A
*-commutativeN/A
lower-*.f64N/A
*-commutativeN/A
lower-*.f6481.0
Applied rewrites81.0%
lift-/.f64N/A
clear-numN/A
associate-/r/N/A
lower-*.f64N/A
lift-*.f64N/A
associate-/r*N/A
lower-/.f64N/A
metadata-eval81.1
lift-+.f64N/A
Applied rewrites81.1%
if 6.70000000000000018 < b Initial program 49.0%
Taylor expanded in b around inf
lower-/.f64N/A
Applied rewrites87.9%
(FPCore (a b c) :precision binary64 (if (<= b 6.7) (* (/ 0.3333333333333333 a) (- (sqrt (fma b b (* -3.0 (* a c)))) b)) (/ (* c (fma -0.375 (* a (/ c (* b b))) -0.5)) b)))
double code(double a, double b, double c) {
double tmp;
if (b <= 6.7) {
tmp = (0.3333333333333333 / a) * (sqrt(fma(b, b, (-3.0 * (a * c)))) - b);
} else {
tmp = (c * fma(-0.375, (a * (c / (b * b))), -0.5)) / b;
}
return tmp;
}
function code(a, b, c) tmp = 0.0 if (b <= 6.7) tmp = Float64(Float64(0.3333333333333333 / a) * Float64(sqrt(fma(b, b, Float64(-3.0 * Float64(a * c)))) - b)); else tmp = Float64(Float64(c * fma(-0.375, Float64(a * Float64(c / Float64(b * b))), -0.5)) / b); end return tmp end
code[a_, b_, c_] := If[LessEqual[b, 6.7], N[(N[(0.3333333333333333 / a), $MachinePrecision] * N[(N[Sqrt[N[(b * b + N[(-3.0 * N[(a * c), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] - b), $MachinePrecision]), $MachinePrecision], N[(N[(c * N[(-0.375 * N[(a * N[(c / N[(b * b), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + -0.5), $MachinePrecision]), $MachinePrecision] / b), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;b \leq 6.7:\\
\;\;\;\;\frac{0.3333333333333333}{a} \cdot \left(\sqrt{\mathsf{fma}\left(b, b, -3 \cdot \left(a \cdot c\right)\right)} - b\right)\\
\mathbf{else}:\\
\;\;\;\;\frac{c \cdot \mathsf{fma}\left(-0.375, a \cdot \frac{c}{b \cdot b}, -0.5\right)}{b}\\
\end{array}
\end{array}
if b < 6.70000000000000018Initial program 80.9%
Applied rewrites80.7%
Taylor expanded in a around 0
+-commutativeN/A
unpow2N/A
lower-fma.f64N/A
*-commutativeN/A
lower-*.f64N/A
*-commutativeN/A
lower-*.f6481.0
Applied rewrites81.0%
lift-/.f64N/A
clear-numN/A
associate-/r/N/A
lower-*.f64N/A
lift-*.f64N/A
associate-/r*N/A
lower-/.f64N/A
metadata-eval81.1
lift-+.f64N/A
Applied rewrites81.1%
if 6.70000000000000018 < b Initial program 49.0%
Taylor expanded in b around inf
Applied rewrites94.7%
Taylor expanded in c around 0
Applied rewrites87.8%
(FPCore (a b c) :precision binary64 (if (<= b 6.7) (* (/ -0.3333333333333333 a) (- b (sqrt (fma a (* -3.0 c) (* b b))))) (/ (* c (fma -0.375 (* a (/ c (* b b))) -0.5)) b)))
double code(double a, double b, double c) {
double tmp;
if (b <= 6.7) {
tmp = (-0.3333333333333333 / a) * (b - sqrt(fma(a, (-3.0 * c), (b * b))));
} else {
tmp = (c * fma(-0.375, (a * (c / (b * b))), -0.5)) / b;
}
return tmp;
}
function code(a, b, c) tmp = 0.0 if (b <= 6.7) tmp = Float64(Float64(-0.3333333333333333 / a) * Float64(b - sqrt(fma(a, Float64(-3.0 * c), Float64(b * b))))); else tmp = Float64(Float64(c * fma(-0.375, Float64(a * Float64(c / Float64(b * b))), -0.5)) / b); end return tmp end
code[a_, b_, c_] := If[LessEqual[b, 6.7], N[(N[(-0.3333333333333333 / a), $MachinePrecision] * N[(b - N[Sqrt[N[(a * N[(-3.0 * c), $MachinePrecision] + N[(b * b), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(c * N[(-0.375 * N[(a * N[(c / N[(b * b), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + -0.5), $MachinePrecision]), $MachinePrecision] / b), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;b \leq 6.7:\\
\;\;\;\;\frac{-0.3333333333333333}{a} \cdot \left(b - \sqrt{\mathsf{fma}\left(a, -3 \cdot c, b \cdot b\right)}\right)\\
\mathbf{else}:\\
\;\;\;\;\frac{c \cdot \mathsf{fma}\left(-0.375, a \cdot \frac{c}{b \cdot b}, -0.5\right)}{b}\\
\end{array}
\end{array}
if b < 6.70000000000000018Initial program 80.9%
Applied rewrites80.9%
if 6.70000000000000018 < b Initial program 49.0%
Taylor expanded in b around inf
Applied rewrites94.7%
Taylor expanded in c around 0
Applied rewrites87.8%
(FPCore (a b c) :precision binary64 (/ (* c (fma -0.375 (* a (/ c (* b b))) -0.5)) b))
double code(double a, double b, double c) {
return (c * fma(-0.375, (a * (c / (b * b))), -0.5)) / b;
}
function code(a, b, c) return Float64(Float64(c * fma(-0.375, Float64(a * Float64(c / Float64(b * b))), -0.5)) / b) end
code[a_, b_, c_] := N[(N[(c * N[(-0.375 * N[(a * N[(c / N[(b * b), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + -0.5), $MachinePrecision]), $MachinePrecision] / b), $MachinePrecision]
\begin{array}{l}
\\
\frac{c \cdot \mathsf{fma}\left(-0.375, a \cdot \frac{c}{b \cdot b}, -0.5\right)}{b}
\end{array}
Initial program 55.9%
Taylor expanded in b around inf
Applied rewrites91.0%
Taylor expanded in c around 0
Applied rewrites81.8%
(FPCore (a b c) :precision binary64 (* -0.5 (/ c b)))
double code(double a, double b, double c) {
return -0.5 * (c / b);
}
real(8) function code(a, b, c)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: c
code = (-0.5d0) * (c / b)
end function
public static double code(double a, double b, double c) {
return -0.5 * (c / b);
}
def code(a, b, c): return -0.5 * (c / b)
function code(a, b, c) return Float64(-0.5 * Float64(c / b)) end
function tmp = code(a, b, c) tmp = -0.5 * (c / b); end
code[a_, b_, c_] := N[(-0.5 * N[(c / b), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
-0.5 \cdot \frac{c}{b}
\end{array}
Initial program 55.9%
Taylor expanded in b around inf
lower-*.f64N/A
lower-/.f6464.1
Applied rewrites64.1%
herbie shell --seed 2024235
(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)))