
(FPCore (a b c) :precision binary64 (let* ((t_0 (sqrt (- (* b b) (* (* 4.0 a) c))))) (if (>= b 0.0) (/ (* 2.0 c) (- (- b) t_0)) (/ (+ (- b) t_0) (* 2.0 a)))))
double code(double a, double b, double c) {
double t_0 = sqrt(((b * b) - ((4.0 * a) * c)));
double tmp;
if (b >= 0.0) {
tmp = (2.0 * c) / (-b - t_0);
} else {
tmp = (-b + t_0) / (2.0 * a);
}
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) - ((4.0d0 * a) * c)))
if (b >= 0.0d0) then
tmp = (2.0d0 * c) / (-b - t_0)
else
tmp = (-b + t_0) / (2.0d0 * a)
end if
code = tmp
end function
public static double code(double a, double b, double c) {
double t_0 = Math.sqrt(((b * b) - ((4.0 * a) * c)));
double tmp;
if (b >= 0.0) {
tmp = (2.0 * c) / (-b - t_0);
} else {
tmp = (-b + t_0) / (2.0 * a);
}
return tmp;
}
def code(a, b, c): t_0 = math.sqrt(((b * b) - ((4.0 * a) * c))) tmp = 0 if b >= 0.0: tmp = (2.0 * c) / (-b - t_0) else: tmp = (-b + t_0) / (2.0 * a) return tmp
function code(a, b, c) t_0 = sqrt(Float64(Float64(b * b) - Float64(Float64(4.0 * a) * c))) tmp = 0.0 if (b >= 0.0) tmp = Float64(Float64(2.0 * c) / Float64(Float64(-b) - t_0)); else tmp = Float64(Float64(Float64(-b) + t_0) / Float64(2.0 * a)); end return tmp end
function tmp_2 = code(a, b, c) t_0 = sqrt(((b * b) - ((4.0 * a) * c))); tmp = 0.0; if (b >= 0.0) tmp = (2.0 * c) / (-b - t_0); else tmp = (-b + t_0) / (2.0 * a); end tmp_2 = tmp; end
code[a_, b_, c_] := Block[{t$95$0 = N[Sqrt[N[(N[(b * b), $MachinePrecision] - N[(N[(4.0 * a), $MachinePrecision] * c), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]}, If[GreaterEqual[b, 0.0], N[(N[(2.0 * c), $MachinePrecision] / N[((-b) - t$95$0), $MachinePrecision]), $MachinePrecision], N[(N[((-b) + t$95$0), $MachinePrecision] / N[(2.0 * a), $MachinePrecision]), $MachinePrecision]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \sqrt{b \cdot b - \left(4 \cdot a\right) \cdot c}\\
\mathbf{if}\;b \geq 0:\\
\;\;\;\;\frac{2 \cdot c}{\left(-b\right) - t\_0}\\
\mathbf{else}:\\
\;\;\;\;\frac{\left(-b\right) + t\_0}{2 \cdot a}\\
\end{array}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 6 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (a b c) :precision binary64 (let* ((t_0 (sqrt (- (* b b) (* (* 4.0 a) c))))) (if (>= b 0.0) (/ (* 2.0 c) (- (- b) t_0)) (/ (+ (- b) t_0) (* 2.0 a)))))
double code(double a, double b, double c) {
double t_0 = sqrt(((b * b) - ((4.0 * a) * c)));
double tmp;
if (b >= 0.0) {
tmp = (2.0 * c) / (-b - t_0);
} else {
tmp = (-b + t_0) / (2.0 * a);
}
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) - ((4.0d0 * a) * c)))
if (b >= 0.0d0) then
tmp = (2.0d0 * c) / (-b - t_0)
else
tmp = (-b + t_0) / (2.0d0 * a)
end if
code = tmp
end function
public static double code(double a, double b, double c) {
double t_0 = Math.sqrt(((b * b) - ((4.0 * a) * c)));
double tmp;
if (b >= 0.0) {
tmp = (2.0 * c) / (-b - t_0);
} else {
tmp = (-b + t_0) / (2.0 * a);
}
return tmp;
}
def code(a, b, c): t_0 = math.sqrt(((b * b) - ((4.0 * a) * c))) tmp = 0 if b >= 0.0: tmp = (2.0 * c) / (-b - t_0) else: tmp = (-b + t_0) / (2.0 * a) return tmp
function code(a, b, c) t_0 = sqrt(Float64(Float64(b * b) - Float64(Float64(4.0 * a) * c))) tmp = 0.0 if (b >= 0.0) tmp = Float64(Float64(2.0 * c) / Float64(Float64(-b) - t_0)); else tmp = Float64(Float64(Float64(-b) + t_0) / Float64(2.0 * a)); end return tmp end
function tmp_2 = code(a, b, c) t_0 = sqrt(((b * b) - ((4.0 * a) * c))); tmp = 0.0; if (b >= 0.0) tmp = (2.0 * c) / (-b - t_0); else tmp = (-b + t_0) / (2.0 * a); end tmp_2 = tmp; end
code[a_, b_, c_] := Block[{t$95$0 = N[Sqrt[N[(N[(b * b), $MachinePrecision] - N[(N[(4.0 * a), $MachinePrecision] * c), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]}, If[GreaterEqual[b, 0.0], N[(N[(2.0 * c), $MachinePrecision] / N[((-b) - t$95$0), $MachinePrecision]), $MachinePrecision], N[(N[((-b) + t$95$0), $MachinePrecision] / N[(2.0 * a), $MachinePrecision]), $MachinePrecision]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \sqrt{b \cdot b - \left(4 \cdot a\right) \cdot c}\\
\mathbf{if}\;b \geq 0:\\
\;\;\;\;\frac{2 \cdot c}{\left(-b\right) - t\_0}\\
\mathbf{else}:\\
\;\;\;\;\frac{\left(-b\right) + t\_0}{2 \cdot a}\\
\end{array}
\end{array}
(FPCore (a b c)
:precision binary64
(let* ((t_0 (sqrt (fma (* c a) -4.0 (* b b)))) (t_1 (- (/ b a))))
(if (<= b -9.8e+141)
(if (>= b 0.0) t_1 t_1)
(if (<= b 1.16e+89)
(if (>= b 0.0) (/ (* c -2.0) (+ b t_0)) (/ (* 0.5 (- t_0 b)) a))
(if (>= b 0.0) (/ c (- b)) t_1)))))
double code(double a, double b, double c) {
double t_0 = sqrt(fma((c * a), -4.0, (b * b)));
double t_1 = -(b / a);
double tmp_1;
if (b <= -9.8e+141) {
double tmp_2;
if (b >= 0.0) {
tmp_2 = t_1;
} else {
tmp_2 = t_1;
}
tmp_1 = tmp_2;
} else if (b <= 1.16e+89) {
double tmp_3;
if (b >= 0.0) {
tmp_3 = (c * -2.0) / (b + t_0);
} else {
tmp_3 = (0.5 * (t_0 - b)) / a;
}
tmp_1 = tmp_3;
} else if (b >= 0.0) {
tmp_1 = c / -b;
} else {
tmp_1 = t_1;
}
return tmp_1;
}
function code(a, b, c) t_0 = sqrt(fma(Float64(c * a), -4.0, Float64(b * b))) t_1 = Float64(-Float64(b / a)) tmp_1 = 0.0 if (b <= -9.8e+141) tmp_2 = 0.0 if (b >= 0.0) tmp_2 = t_1; else tmp_2 = t_1; end tmp_1 = tmp_2; elseif (b <= 1.16e+89) tmp_3 = 0.0 if (b >= 0.0) tmp_3 = Float64(Float64(c * -2.0) / Float64(b + t_0)); else tmp_3 = Float64(Float64(0.5 * Float64(t_0 - b)) / a); end tmp_1 = tmp_3; elseif (b >= 0.0) tmp_1 = Float64(c / Float64(-b)); else tmp_1 = t_1; end return tmp_1 end
code[a_, b_, c_] := Block[{t$95$0 = N[Sqrt[N[(N[(c * a), $MachinePrecision] * -4.0 + N[(b * b), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]}, Block[{t$95$1 = (-N[(b / a), $MachinePrecision])}, If[LessEqual[b, -9.8e+141], If[GreaterEqual[b, 0.0], t$95$1, t$95$1], If[LessEqual[b, 1.16e+89], If[GreaterEqual[b, 0.0], N[(N[(c * -2.0), $MachinePrecision] / N[(b + t$95$0), $MachinePrecision]), $MachinePrecision], N[(N[(0.5 * N[(t$95$0 - b), $MachinePrecision]), $MachinePrecision] / a), $MachinePrecision]], If[GreaterEqual[b, 0.0], N[(c / (-b)), $MachinePrecision], t$95$1]]]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \sqrt{\mathsf{fma}\left(c \cdot a, -4, b \cdot b\right)}\\
t_1 := -\frac{b}{a}\\
\mathbf{if}\;b \leq -9.8 \cdot 10^{+141}:\\
\;\;\;\;\begin{array}{l}
\mathbf{if}\;b \geq 0:\\
\;\;\;\;t\_1\\
\mathbf{else}:\\
\;\;\;\;t\_1\\
\end{array}\\
\mathbf{elif}\;b \leq 1.16 \cdot 10^{+89}:\\
\;\;\;\;\begin{array}{l}
\mathbf{if}\;b \geq 0:\\
\;\;\;\;\frac{c \cdot -2}{b + t\_0}\\
\mathbf{else}:\\
\;\;\;\;\frac{0.5 \cdot \left(t\_0 - b\right)}{a}\\
\end{array}\\
\mathbf{elif}\;b \geq 0:\\
\;\;\;\;\frac{c}{-b}\\
\mathbf{else}:\\
\;\;\;\;t\_1\\
\end{array}
\end{array}
if b < -9.8000000000000002e141Initial program 45.0%
Taylor expanded in c around 0
mul-1-negN/A
distribute-neg-frac2N/A
lower-/.f64N/A
lower-neg.f6445.0
Applied rewrites45.0%
Taylor expanded in b around -inf
mul-1-negN/A
distribute-neg-frac2N/A
mul-1-negN/A
lower-/.f64N/A
mul-1-negN/A
lower-neg.f6493.9
Applied rewrites93.9%
Taylor expanded in b around -inf
mul-1-negN/A
distribute-neg-frac2N/A
mul-1-negN/A
lower-/.f64N/A
mul-1-negN/A
lower-neg.f6493.9
Applied rewrites93.9%
if -9.8000000000000002e141 < b < 1.16e89Initial program 85.5%
Taylor expanded in c around 0
mul-1-negN/A
distribute-neg-frac2N/A
lower-/.f64N/A
lower-neg.f6469.6
Applied rewrites69.6%
Taylor expanded in b around 0
lower->=.f64N/A
associate-*r/N/A
lower-/.f64N/A
lower-*.f64N/A
lower-+.f64N/A
lower-sqrt.f64N/A
cancel-sign-sub-invN/A
metadata-evalN/A
+-commutativeN/A
*-commutativeN/A
lower-fma.f64N/A
lower-*.f64N/A
unpow2N/A
lower-*.f64N/A
Applied rewrites85.5%
if 1.16e89 < b Initial program 52.6%
Taylor expanded in c around 0
mul-1-negN/A
distribute-neg-frac2N/A
lower-/.f64N/A
lower-neg.f6497.2
Applied rewrites97.2%
Taylor expanded in b around -inf
mul-1-negN/A
distribute-neg-frac2N/A
mul-1-negN/A
lower-/.f64N/A
mul-1-negN/A
lower-neg.f6497.2
Applied rewrites97.2%
Final simplification90.0%
(FPCore (a b c)
:precision binary64
(let* ((t_0 (- (/ b a))))
(if (<= b -1.4e-77)
(if (>= b 0.0) t_0 t_0)
(if (>= b 0.0) (/ c (- b)) (/ (- (sqrt (* (* c a) -4.0)) b) (* 2.0 a))))))
double code(double a, double b, double c) {
double t_0 = -(b / a);
double tmp_1;
if (b <= -1.4e-77) {
double tmp_2;
if (b >= 0.0) {
tmp_2 = t_0;
} else {
tmp_2 = t_0;
}
tmp_1 = tmp_2;
} else if (b >= 0.0) {
tmp_1 = c / -b;
} else {
tmp_1 = (sqrt(((c * a) * -4.0)) - b) / (2.0 * a);
}
return tmp_1;
}
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
real(8) :: tmp_1
real(8) :: tmp_2
t_0 = -(b / a)
if (b <= (-1.4d-77)) then
if (b >= 0.0d0) then
tmp_2 = t_0
else
tmp_2 = t_0
end if
tmp_1 = tmp_2
else if (b >= 0.0d0) then
tmp_1 = c / -b
else
tmp_1 = (sqrt(((c * a) * (-4.0d0))) - b) / (2.0d0 * a)
end if
code = tmp_1
end function
public static double code(double a, double b, double c) {
double t_0 = -(b / a);
double tmp_1;
if (b <= -1.4e-77) {
double tmp_2;
if (b >= 0.0) {
tmp_2 = t_0;
} else {
tmp_2 = t_0;
}
tmp_1 = tmp_2;
} else if (b >= 0.0) {
tmp_1 = c / -b;
} else {
tmp_1 = (Math.sqrt(((c * a) * -4.0)) - b) / (2.0 * a);
}
return tmp_1;
}
def code(a, b, c): t_0 = -(b / a) tmp_1 = 0 if b <= -1.4e-77: tmp_2 = 0 if b >= 0.0: tmp_2 = t_0 else: tmp_2 = t_0 tmp_1 = tmp_2 elif b >= 0.0: tmp_1 = c / -b else: tmp_1 = (math.sqrt(((c * a) * -4.0)) - b) / (2.0 * a) return tmp_1
function code(a, b, c) t_0 = Float64(-Float64(b / a)) tmp_1 = 0.0 if (b <= -1.4e-77) tmp_2 = 0.0 if (b >= 0.0) tmp_2 = t_0; else tmp_2 = t_0; end tmp_1 = tmp_2; elseif (b >= 0.0) tmp_1 = Float64(c / Float64(-b)); else tmp_1 = Float64(Float64(sqrt(Float64(Float64(c * a) * -4.0)) - b) / Float64(2.0 * a)); end return tmp_1 end
function tmp_4 = code(a, b, c) t_0 = -(b / a); tmp_2 = 0.0; if (b <= -1.4e-77) tmp_3 = 0.0; if (b >= 0.0) tmp_3 = t_0; else tmp_3 = t_0; end tmp_2 = tmp_3; elseif (b >= 0.0) tmp_2 = c / -b; else tmp_2 = (sqrt(((c * a) * -4.0)) - b) / (2.0 * a); end tmp_4 = tmp_2; end
code[a_, b_, c_] := Block[{t$95$0 = (-N[(b / a), $MachinePrecision])}, If[LessEqual[b, -1.4e-77], If[GreaterEqual[b, 0.0], t$95$0, t$95$0], If[GreaterEqual[b, 0.0], N[(c / (-b)), $MachinePrecision], N[(N[(N[Sqrt[N[(N[(c * a), $MachinePrecision] * -4.0), $MachinePrecision]], $MachinePrecision] - b), $MachinePrecision] / N[(2.0 * a), $MachinePrecision]), $MachinePrecision]]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := -\frac{b}{a}\\
\mathbf{if}\;b \leq -1.4 \cdot 10^{-77}:\\
\;\;\;\;\begin{array}{l}
\mathbf{if}\;b \geq 0:\\
\;\;\;\;t\_0\\
\mathbf{else}:\\
\;\;\;\;t\_0\\
\end{array}\\
\mathbf{elif}\;b \geq 0:\\
\;\;\;\;\frac{c}{-b}\\
\mathbf{else}:\\
\;\;\;\;\frac{\sqrt{\left(c \cdot a\right) \cdot -4} - b}{2 \cdot a}\\
\end{array}
\end{array}
if b < -1.4e-77Initial program 68.5%
Taylor expanded in c around 0
mul-1-negN/A
distribute-neg-frac2N/A
lower-/.f64N/A
lower-neg.f6468.5
Applied rewrites68.5%
Taylor expanded in b around -inf
mul-1-negN/A
distribute-neg-frac2N/A
mul-1-negN/A
lower-/.f64N/A
mul-1-negN/A
lower-neg.f6485.5
Applied rewrites85.5%
Taylor expanded in b around -inf
mul-1-negN/A
distribute-neg-frac2N/A
mul-1-negN/A
lower-/.f64N/A
mul-1-negN/A
lower-neg.f6485.5
Applied rewrites85.5%
if -1.4e-77 < b Initial program 74.7%
Taylor expanded in c around 0
mul-1-negN/A
distribute-neg-frac2N/A
lower-/.f64N/A
lower-neg.f6469.8
Applied rewrites69.8%
Taylor expanded in b around 0
*-commutativeN/A
lower-*.f64N/A
lower-*.f6467.5
Applied rewrites67.5%
Final simplification74.2%
herbie shell --seed 2024229
(FPCore (a b c)
:name "jeff quadratic root 2"
:precision binary64
(if (>= b 0.0) (/ (* 2.0 c) (- (- b) (sqrt (- (* b b) (* (* 4.0 a) c))))) (/ (+ (- b) (sqrt (- (* b b) (* (* 4.0 a) c)))) (* 2.0 a))))