
(FPCore (a b c) :precision binary64 (/ (+ (- b) (sqrt (- (* b b) (* (* 4.0 a) c)))) (* 2.0 a)))
double code(double a, double b, double c) {
return (-b + sqrt(((b * b) - ((4.0 * a) * c)))) / (2.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) - ((4.0d0 * a) * c)))) / (2.0d0 * a)
end function
public static double code(double a, double b, double c) {
return (-b + Math.sqrt(((b * b) - ((4.0 * a) * c)))) / (2.0 * a);
}
def code(a, b, c): return (-b + math.sqrt(((b * b) - ((4.0 * a) * c)))) / (2.0 * a)
function code(a, b, c) return Float64(Float64(Float64(-b) + sqrt(Float64(Float64(b * b) - Float64(Float64(4.0 * a) * c)))) / Float64(2.0 * a)) end
function tmp = code(a, b, c) tmp = (-b + sqrt(((b * b) - ((4.0 * a) * c)))) / (2.0 * a); end
code[a_, b_, c_] := N[(N[((-b) + N[Sqrt[N[(N[(b * b), $MachinePrecision] - N[(N[(4.0 * a), $MachinePrecision] * c), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] / N[(2.0 * a), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{\left(-b\right) + \sqrt{b \cdot b - \left(4 \cdot a\right) \cdot c}}{2 \cdot a}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 6 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (a b c) :precision binary64 (/ (+ (- b) (sqrt (- (* b b) (* (* 4.0 a) c)))) (* 2.0 a)))
double code(double a, double b, double c) {
return (-b + sqrt(((b * b) - ((4.0 * a) * c)))) / (2.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) - ((4.0d0 * a) * c)))) / (2.0d0 * a)
end function
public static double code(double a, double b, double c) {
return (-b + Math.sqrt(((b * b) - ((4.0 * a) * c)))) / (2.0 * a);
}
def code(a, b, c): return (-b + math.sqrt(((b * b) - ((4.0 * a) * c)))) / (2.0 * a)
function code(a, b, c) return Float64(Float64(Float64(-b) + sqrt(Float64(Float64(b * b) - Float64(Float64(4.0 * a) * c)))) / Float64(2.0 * a)) end
function tmp = code(a, b, c) tmp = (-b + sqrt(((b * b) - ((4.0 * a) * c)))) / (2.0 * a); end
code[a_, b_, c_] := N[(N[((-b) + N[Sqrt[N[(N[(b * b), $MachinePrecision] - N[(N[(4.0 * a), $MachinePrecision] * c), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] / N[(2.0 * a), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{\left(-b\right) + \sqrt{b \cdot b - \left(4 \cdot a\right) \cdot c}}{2 \cdot a}
\end{array}
(FPCore (a b c)
:precision binary64
(let* ((t_0 (* c (* a 4.0))))
(/
(* (/ (* a c) a) 2.0)
(-
(- b)
(sqrt
(/
(fma -64.0 (pow (* a c) 3.0) (pow b 6.0))
(fma t_0 (fma b b t_0) (pow b 4.0))))))))
double code(double a, double b, double c) {
double t_0 = c * (a * 4.0);
return (((a * c) / a) * 2.0) / (-b - sqrt((fma(-64.0, pow((a * c), 3.0), pow(b, 6.0)) / fma(t_0, fma(b, b, t_0), pow(b, 4.0)))));
}
function code(a, b, c) t_0 = Float64(c * Float64(a * 4.0)) return Float64(Float64(Float64(Float64(a * c) / a) * 2.0) / Float64(Float64(-b) - sqrt(Float64(fma(-64.0, (Float64(a * c) ^ 3.0), (b ^ 6.0)) / fma(t_0, fma(b, b, t_0), (b ^ 4.0)))))) end
code[a_, b_, c_] := Block[{t$95$0 = N[(c * N[(a * 4.0), $MachinePrecision]), $MachinePrecision]}, N[(N[(N[(N[(a * c), $MachinePrecision] / a), $MachinePrecision] * 2.0), $MachinePrecision] / N[((-b) - N[Sqrt[N[(N[(-64.0 * N[Power[N[(a * c), $MachinePrecision], 3.0], $MachinePrecision] + N[Power[b, 6.0], $MachinePrecision]), $MachinePrecision] / N[(t$95$0 * N[(b * b + t$95$0), $MachinePrecision] + N[Power[b, 4.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := c \cdot \left(a \cdot 4\right)\\
\frac{\frac{a \cdot c}{a} \cdot 2}{\left(-b\right) - \sqrt{\frac{\mathsf{fma}\left(-64, {\left(a \cdot c\right)}^{3}, {b}^{6}\right)}{\mathsf{fma}\left(t_0, \mathsf{fma}\left(b, b, t_0\right), {b}^{4}\right)}}}
\end{array}
\end{array}
Initial program 28.8%
*-commutative28.8%
Simplified28.8%
flip3--28.9%
div-inv28.9%
pow228.9%
pow-pow28.9%
metadata-eval28.9%
associate-*l*28.9%
unpow-prod-down28.9%
metadata-eval28.9%
pow228.9%
pow228.9%
pow-prod-up28.8%
metadata-eval28.8%
distribute-rgt-out28.8%
Applied egg-rr28.8%
flip-+28.8%
Applied egg-rr29.7%
Taylor expanded in b around 0 99.4%
associate-*r*99.4%
*-commutative99.4%
*-commutative99.4%
Simplified99.4%
Applied egg-rr37.2%
expm1-def86.1%
expm1-log1p99.3%
associate-/r*99.6%
associate-*r*99.6%
*-commutative99.6%
times-frac99.6%
metadata-eval99.6%
*-commutative99.6%
Simplified99.6%
Final simplification99.6%
(FPCore (a b c) :precision binary64 (let* ((t_0 (/ (- (sqrt (- (* b b) (* c (* a 4.0)))) b) (* a 2.0)))) (if (<= t_0 -2e-10) t_0 (/ (- c) b))))
double code(double a, double b, double c) {
double t_0 = (sqrt(((b * b) - (c * (a * 4.0)))) - b) / (a * 2.0);
double tmp;
if (t_0 <= -2e-10) {
tmp = t_0;
} else {
tmp = -c / b;
}
return tmp;
}
real(8) function code(a, b, c)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: c
real(8) :: t_0
real(8) :: tmp
t_0 = (sqrt(((b * b) - (c * (a * 4.0d0)))) - b) / (a * 2.0d0)
if (t_0 <= (-2d-10)) then
tmp = t_0
else
tmp = -c / b
end if
code = tmp
end function
public static double code(double a, double b, double c) {
double t_0 = (Math.sqrt(((b * b) - (c * (a * 4.0)))) - b) / (a * 2.0);
double tmp;
if (t_0 <= -2e-10) {
tmp = t_0;
} else {
tmp = -c / b;
}
return tmp;
}
def code(a, b, c): t_0 = (math.sqrt(((b * b) - (c * (a * 4.0)))) - b) / (a * 2.0) tmp = 0 if t_0 <= -2e-10: tmp = t_0 else: tmp = -c / b return tmp
function code(a, b, c) t_0 = Float64(Float64(sqrt(Float64(Float64(b * b) - Float64(c * Float64(a * 4.0)))) - b) / Float64(a * 2.0)) tmp = 0.0 if (t_0 <= -2e-10) tmp = t_0; else tmp = Float64(Float64(-c) / b); end return tmp end
function tmp_2 = code(a, b, c) t_0 = (sqrt(((b * b) - (c * (a * 4.0)))) - b) / (a * 2.0); tmp = 0.0; if (t_0 <= -2e-10) tmp = t_0; else tmp = -c / b; end tmp_2 = tmp; end
code[a_, b_, c_] := Block[{t$95$0 = N[(N[(N[Sqrt[N[(N[(b * b), $MachinePrecision] - N[(c * N[(a * 4.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] - b), $MachinePrecision] / N[(a * 2.0), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[t$95$0, -2e-10], t$95$0, N[((-c) / b), $MachinePrecision]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \frac{\sqrt{b \cdot b - c \cdot \left(a \cdot 4\right)} - b}{a \cdot 2}\\
\mathbf{if}\;t_0 \leq -2 \cdot 10^{-10}:\\
\;\;\;\;t_0\\
\mathbf{else}:\\
\;\;\;\;\frac{-c}{b}\\
\end{array}
\end{array}
if (/.f64 (+.f64 (neg.f64 b) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 4 a) c)))) (*.f64 2 a)) < -2.00000000000000007e-10Initial program 67.2%
if -2.00000000000000007e-10 < (/.f64 (+.f64 (neg.f64 b) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 4 a) c)))) (*.f64 2 a)) Initial program 9.4%
*-commutative9.4%
Simplified9.4%
Taylor expanded in b around inf 96.6%
mul-1-neg96.6%
distribute-neg-frac96.6%
Simplified96.6%
Final simplification86.7%
(FPCore (a b c) :precision binary64 (/ (/ (* c (* a 4.0)) (- (- b) (sqrt (+ (* (* a c) -4.0) (pow b 2.0))))) (* a 2.0)))
double code(double a, double b, double c) {
return ((c * (a * 4.0)) / (-b - sqrt((((a * c) * -4.0) + pow(b, 2.0))))) / (a * 2.0);
}
real(8) function code(a, b, c)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: c
code = ((c * (a * 4.0d0)) / (-b - sqrt((((a * c) * (-4.0d0)) + (b ** 2.0d0))))) / (a * 2.0d0)
end function
public static double code(double a, double b, double c) {
return ((c * (a * 4.0)) / (-b - Math.sqrt((((a * c) * -4.0) + Math.pow(b, 2.0))))) / (a * 2.0);
}
def code(a, b, c): return ((c * (a * 4.0)) / (-b - math.sqrt((((a * c) * -4.0) + math.pow(b, 2.0))))) / (a * 2.0)
function code(a, b, c) return Float64(Float64(Float64(c * Float64(a * 4.0)) / Float64(Float64(-b) - sqrt(Float64(Float64(Float64(a * c) * -4.0) + (b ^ 2.0))))) / Float64(a * 2.0)) end
function tmp = code(a, b, c) tmp = ((c * (a * 4.0)) / (-b - sqrt((((a * c) * -4.0) + (b ^ 2.0))))) / (a * 2.0); end
code[a_, b_, c_] := N[(N[(N[(c * N[(a * 4.0), $MachinePrecision]), $MachinePrecision] / N[((-b) - N[Sqrt[N[(N[(N[(a * c), $MachinePrecision] * -4.0), $MachinePrecision] + N[Power[b, 2.0], $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(a * 2.0), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{\frac{c \cdot \left(a \cdot 4\right)}{\left(-b\right) - \sqrt{\left(a \cdot c\right) \cdot -4 + {b}^{2}}}}{a \cdot 2}
\end{array}
Initial program 28.8%
*-commutative28.8%
Simplified28.8%
flip3--28.9%
div-inv28.9%
pow228.9%
pow-pow28.9%
metadata-eval28.9%
associate-*l*28.9%
unpow-prod-down28.9%
metadata-eval28.9%
pow228.9%
pow228.9%
pow-prod-up28.8%
metadata-eval28.8%
distribute-rgt-out28.8%
Applied egg-rr28.8%
flip-+28.8%
Applied egg-rr29.7%
Taylor expanded in b around 0 99.4%
associate-*r*99.4%
*-commutative99.4%
*-commutative99.4%
Simplified99.4%
Taylor expanded in b around 0 99.5%
Final simplification99.5%
(FPCore (a b c) :precision binary64 (if (<= b 0.000225) (/ (- (sqrt (fma b b (* c (* a -4.0)))) b) (* a 2.0)) (- (/ (- c) b) (/ a (/ (pow b 3.0) (pow c 2.0))))))
double code(double a, double b, double c) {
double tmp;
if (b <= 0.000225) {
tmp = (sqrt(fma(b, b, (c * (a * -4.0)))) - b) / (a * 2.0);
} else {
tmp = (-c / b) - (a / (pow(b, 3.0) / pow(c, 2.0)));
}
return tmp;
}
function code(a, b, c) tmp = 0.0 if (b <= 0.000225) tmp = Float64(Float64(sqrt(fma(b, b, Float64(c * Float64(a * -4.0)))) - b) / Float64(a * 2.0)); else tmp = Float64(Float64(Float64(-c) / b) - Float64(a / Float64((b ^ 3.0) / (c ^ 2.0)))); end return tmp end
code[a_, b_, c_] := If[LessEqual[b, 0.000225], N[(N[(N[Sqrt[N[(b * b + N[(c * N[(a * -4.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] - b), $MachinePrecision] / N[(a * 2.0), $MachinePrecision]), $MachinePrecision], N[(N[((-c) / b), $MachinePrecision] - N[(a / N[(N[Power[b, 3.0], $MachinePrecision] / N[Power[c, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;b \leq 0.000225:\\
\;\;\;\;\frac{\sqrt{\mathsf{fma}\left(b, b, c \cdot \left(a \cdot -4\right)\right)} - b}{a \cdot 2}\\
\mathbf{else}:\\
\;\;\;\;\frac{-c}{b} - \frac{a}{\frac{{b}^{3}}{{c}^{2}}}\\
\end{array}
\end{array}
if b < 2.2499999999999999e-4Initial program 81.2%
Simplified81.3%
if 2.2499999999999999e-4 < b Initial program 24.4%
*-commutative24.4%
Simplified24.4%
Taylor expanded in b around inf 94.5%
mul-1-neg94.5%
unsub-neg94.5%
mul-1-neg94.5%
distribute-neg-frac94.5%
associate-/l*94.5%
Simplified94.5%
Final simplification93.5%
(FPCore (a b c) :precision binary64 (if (<= b 0.000225) (/ (- (sqrt (- (* b b) (* c (* a 4.0)))) b) (* a 2.0)) (- (/ (- c) b) (/ a (/ (pow b 3.0) (pow c 2.0))))))
double code(double a, double b, double c) {
double tmp;
if (b <= 0.000225) {
tmp = (sqrt(((b * b) - (c * (a * 4.0)))) - b) / (a * 2.0);
} else {
tmp = (-c / b) - (a / (pow(b, 3.0) / pow(c, 2.0)));
}
return tmp;
}
real(8) function code(a, b, c)
real(8), intent (in) :: a
real(8), intent (in) :: b
real(8), intent (in) :: c
real(8) :: tmp
if (b <= 0.000225d0) then
tmp = (sqrt(((b * b) - (c * (a * 4.0d0)))) - b) / (a * 2.0d0)
else
tmp = (-c / b) - (a / ((b ** 3.0d0) / (c ** 2.0d0)))
end if
code = tmp
end function
public static double code(double a, double b, double c) {
double tmp;
if (b <= 0.000225) {
tmp = (Math.sqrt(((b * b) - (c * (a * 4.0)))) - b) / (a * 2.0);
} else {
tmp = (-c / b) - (a / (Math.pow(b, 3.0) / Math.pow(c, 2.0)));
}
return tmp;
}
def code(a, b, c): tmp = 0 if b <= 0.000225: tmp = (math.sqrt(((b * b) - (c * (a * 4.0)))) - b) / (a * 2.0) else: tmp = (-c / b) - (a / (math.pow(b, 3.0) / math.pow(c, 2.0))) return tmp
function code(a, b, c) tmp = 0.0 if (b <= 0.000225) tmp = Float64(Float64(sqrt(Float64(Float64(b * b) - Float64(c * Float64(a * 4.0)))) - b) / Float64(a * 2.0)); else tmp = Float64(Float64(Float64(-c) / b) - Float64(a / Float64((b ^ 3.0) / (c ^ 2.0)))); end return tmp end
function tmp_2 = code(a, b, c) tmp = 0.0; if (b <= 0.000225) tmp = (sqrt(((b * b) - (c * (a * 4.0)))) - b) / (a * 2.0); else tmp = (-c / b) - (a / ((b ^ 3.0) / (c ^ 2.0))); end tmp_2 = tmp; end
code[a_, b_, c_] := If[LessEqual[b, 0.000225], N[(N[(N[Sqrt[N[(N[(b * b), $MachinePrecision] - N[(c * N[(a * 4.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] - b), $MachinePrecision] / N[(a * 2.0), $MachinePrecision]), $MachinePrecision], N[(N[((-c) / b), $MachinePrecision] - N[(a / N[(N[Power[b, 3.0], $MachinePrecision] / N[Power[c, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;b \leq 0.000225:\\
\;\;\;\;\frac{\sqrt{b \cdot b - c \cdot \left(a \cdot 4\right)} - b}{a \cdot 2}\\
\mathbf{else}:\\
\;\;\;\;\frac{-c}{b} - \frac{a}{\frac{{b}^{3}}{{c}^{2}}}\\
\end{array}
\end{array}
if b < 2.2499999999999999e-4Initial program 81.2%
if 2.2499999999999999e-4 < b Initial program 24.4%
*-commutative24.4%
Simplified24.4%
Taylor expanded in b around inf 94.5%
mul-1-neg94.5%
unsub-neg94.5%
mul-1-neg94.5%
distribute-neg-frac94.5%
associate-/l*94.5%
Simplified94.5%
Final simplification93.5%
(FPCore (a b c) :precision binary64 (/ (- c) b))
double code(double a, double b, double c) {
return -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 = -c / b
end function
public static double code(double a, double b, double c) {
return -c / b;
}
def code(a, b, c): return -c / b
function code(a, b, c) return Float64(Float64(-c) / b) end
function tmp = code(a, b, c) tmp = -c / b; end
code[a_, b_, c_] := N[((-c) / b), $MachinePrecision]
\begin{array}{l}
\\
\frac{-c}{b}
\end{array}
Initial program 28.8%
*-commutative28.8%
Simplified28.8%
Taylor expanded in b around inf 82.6%
mul-1-neg82.6%
distribute-neg-frac82.6%
Simplified82.6%
Final simplification82.6%
herbie shell --seed 2023309
(FPCore (a b c)
:name "Quadratic roots, medium range"
:precision binary64
:pre (and (and (and (< 1.1102230246251565e-16 a) (< a 9007199254740992.0)) (and (< 1.1102230246251565e-16 b) (< b 9007199254740992.0))) (and (< 1.1102230246251565e-16 c) (< c 9007199254740992.0)))
(/ (+ (- b) (sqrt (- (* b b) (* (* 4.0 a) c)))) (* 2.0 a)))