
(FPCore (p r q) :precision binary64 (* (/ 1.0 2.0) (- (+ (fabs p) (fabs r)) (sqrt (+ (pow (- p r) 2.0) (* 4.0 (pow q 2.0)))))))
double code(double p, double r, double q) {
return (1.0 / 2.0) * ((fabs(p) + fabs(r)) - sqrt((pow((p - r), 2.0) + (4.0 * pow(q, 2.0)))));
}
real(8) function code(p, r, q)
real(8), intent (in) :: p
real(8), intent (in) :: r
real(8), intent (in) :: q
code = (1.0d0 / 2.0d0) * ((abs(p) + abs(r)) - sqrt((((p - r) ** 2.0d0) + (4.0d0 * (q ** 2.0d0)))))
end function
public static double code(double p, double r, double q) {
return (1.0 / 2.0) * ((Math.abs(p) + Math.abs(r)) - Math.sqrt((Math.pow((p - r), 2.0) + (4.0 * Math.pow(q, 2.0)))));
}
def code(p, r, q): return (1.0 / 2.0) * ((math.fabs(p) + math.fabs(r)) - math.sqrt((math.pow((p - r), 2.0) + (4.0 * math.pow(q, 2.0)))))
function code(p, r, q) return Float64(Float64(1.0 / 2.0) * Float64(Float64(abs(p) + abs(r)) - sqrt(Float64((Float64(p - r) ^ 2.0) + Float64(4.0 * (q ^ 2.0)))))) end
function tmp = code(p, r, q) tmp = (1.0 / 2.0) * ((abs(p) + abs(r)) - sqrt((((p - r) ^ 2.0) + (4.0 * (q ^ 2.0))))); end
code[p_, r_, q_] := N[(N[(1.0 / 2.0), $MachinePrecision] * N[(N[(N[Abs[p], $MachinePrecision] + N[Abs[r], $MachinePrecision]), $MachinePrecision] - N[Sqrt[N[(N[Power[N[(p - r), $MachinePrecision], 2.0], $MachinePrecision] + N[(4.0 * N[Power[q, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{1}{2} \cdot \left(\left(\left|p\right| + \left|r\right|\right) - \sqrt{{\left(p - r\right)}^{2} + 4 \cdot {q}^{2}}\right)
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 4 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (p r q) :precision binary64 (* (/ 1.0 2.0) (- (+ (fabs p) (fabs r)) (sqrt (+ (pow (- p r) 2.0) (* 4.0 (pow q 2.0)))))))
double code(double p, double r, double q) {
return (1.0 / 2.0) * ((fabs(p) + fabs(r)) - sqrt((pow((p - r), 2.0) + (4.0 * pow(q, 2.0)))));
}
real(8) function code(p, r, q)
real(8), intent (in) :: p
real(8), intent (in) :: r
real(8), intent (in) :: q
code = (1.0d0 / 2.0d0) * ((abs(p) + abs(r)) - sqrt((((p - r) ** 2.0d0) + (4.0d0 * (q ** 2.0d0)))))
end function
public static double code(double p, double r, double q) {
return (1.0 / 2.0) * ((Math.abs(p) + Math.abs(r)) - Math.sqrt((Math.pow((p - r), 2.0) + (4.0 * Math.pow(q, 2.0)))));
}
def code(p, r, q): return (1.0 / 2.0) * ((math.fabs(p) + math.fabs(r)) - math.sqrt((math.pow((p - r), 2.0) + (4.0 * math.pow(q, 2.0)))))
function code(p, r, q) return Float64(Float64(1.0 / 2.0) * Float64(Float64(abs(p) + abs(r)) - sqrt(Float64((Float64(p - r) ^ 2.0) + Float64(4.0 * (q ^ 2.0)))))) end
function tmp = code(p, r, q) tmp = (1.0 / 2.0) * ((abs(p) + abs(r)) - sqrt((((p - r) ^ 2.0) + (4.0 * (q ^ 2.0))))); end
code[p_, r_, q_] := N[(N[(1.0 / 2.0), $MachinePrecision] * N[(N[(N[Abs[p], $MachinePrecision] + N[Abs[r], $MachinePrecision]), $MachinePrecision] - N[Sqrt[N[(N[Power[N[(p - r), $MachinePrecision], 2.0], $MachinePrecision] + N[(4.0 * N[Power[q, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{1}{2} \cdot \left(\left(\left|p\right| + \left|r\right|\right) - \sqrt{{\left(p - r\right)}^{2} + 4 \cdot {q}^{2}}\right)
\end{array}
q_m = (fabs.f64 q) NOTE: p, r, and q_m should be sorted in increasing order before calling this function. (FPCore (p r q_m) :precision binary64 (if (<= (pow q_m 2.0) 20000000000000.0) (* (- p) (fma (/ (+ (- (fabs r) r) (fabs p)) p) -0.5 -0.5)) (fma 0.5 (+ (fabs r) (fabs p)) (- q_m))))
q_m = fabs(q);
assert(p < r && r < q_m);
double code(double p, double r, double q_m) {
double tmp;
if (pow(q_m, 2.0) <= 20000000000000.0) {
tmp = -p * fma((((fabs(r) - r) + fabs(p)) / p), -0.5, -0.5);
} else {
tmp = fma(0.5, (fabs(r) + fabs(p)), -q_m);
}
return tmp;
}
q_m = abs(q) p, r, q_m = sort([p, r, q_m]) function code(p, r, q_m) tmp = 0.0 if ((q_m ^ 2.0) <= 20000000000000.0) tmp = Float64(Float64(-p) * fma(Float64(Float64(Float64(abs(r) - r) + abs(p)) / p), -0.5, -0.5)); else tmp = fma(0.5, Float64(abs(r) + abs(p)), Float64(-q_m)); end return tmp end
q_m = N[Abs[q], $MachinePrecision] NOTE: p, r, and q_m should be sorted in increasing order before calling this function. code[p_, r_, q$95$m_] := If[LessEqual[N[Power[q$95$m, 2.0], $MachinePrecision], 20000000000000.0], N[((-p) * N[(N[(N[(N[(N[Abs[r], $MachinePrecision] - r), $MachinePrecision] + N[Abs[p], $MachinePrecision]), $MachinePrecision] / p), $MachinePrecision] * -0.5 + -0.5), $MachinePrecision]), $MachinePrecision], N[(0.5 * N[(N[Abs[r], $MachinePrecision] + N[Abs[p], $MachinePrecision]), $MachinePrecision] + (-q$95$m)), $MachinePrecision]]
\begin{array}{l}
q_m = \left|q\right|
\\
[p, r, q_m] = \mathsf{sort}([p, r, q_m])\\
\\
\begin{array}{l}
\mathbf{if}\;{q\_m}^{2} \leq 20000000000000:\\
\;\;\;\;\left(-p\right) \cdot \mathsf{fma}\left(\frac{\left(\left|r\right| - r\right) + \left|p\right|}{p}, -0.5, -0.5\right)\\
\mathbf{else}:\\
\;\;\;\;\mathsf{fma}\left(0.5, \left|r\right| + \left|p\right|, -q\_m\right)\\
\end{array}
\end{array}
if (pow.f64 q #s(literal 2 binary64)) < 2e13Initial program 23.2%
Taylor expanded in q around inf
*-commutativeN/A
lower-*.f64N/A
sub-negN/A
*-commutativeN/A
metadata-evalN/A
lower-fma.f64N/A
lower-/.f64N/A
+-commutativeN/A
lower-+.f64N/A
lower-fabs.f64N/A
lower-fabs.f646.9
Applied rewrites6.9%
Taylor expanded in p around -inf
associate-*r*N/A
lower-*.f64N/A
mul-1-negN/A
lower-neg.f64N/A
sub-negN/A
*-commutativeN/A
metadata-evalN/A
lower-fma.f64N/A
lower-/.f64N/A
associate--l+N/A
+-commutativeN/A
lower-+.f64N/A
lower--.f64N/A
lower-fabs.f64N/A
lower-fabs.f6428.8
Applied rewrites28.8%
if 2e13 < (pow.f64 q #s(literal 2 binary64)) Initial program 27.7%
Taylor expanded in q around inf
*-commutativeN/A
lower-*.f64N/A
sub-negN/A
*-commutativeN/A
metadata-evalN/A
lower-fma.f64N/A
lower-/.f64N/A
+-commutativeN/A
lower-+.f64N/A
lower-fabs.f64N/A
lower-fabs.f6441.7
Applied rewrites41.7%
Taylor expanded in q around 0
Applied rewrites41.7%
q_m = (fabs.f64 q) NOTE: p, r, and q_m should be sorted in increasing order before calling this function. (FPCore (p r q_m) :precision binary64 (if (<= (pow q_m 2.0) 20000000000000.0) (/ (* (- q_m) q_m) r) (fma 0.5 (+ (fabs r) (fabs p)) (- q_m))))
q_m = fabs(q);
assert(p < r && r < q_m);
double code(double p, double r, double q_m) {
double tmp;
if (pow(q_m, 2.0) <= 20000000000000.0) {
tmp = (-q_m * q_m) / r;
} else {
tmp = fma(0.5, (fabs(r) + fabs(p)), -q_m);
}
return tmp;
}
q_m = abs(q) p, r, q_m = sort([p, r, q_m]) function code(p, r, q_m) tmp = 0.0 if ((q_m ^ 2.0) <= 20000000000000.0) tmp = Float64(Float64(Float64(-q_m) * q_m) / r); else tmp = fma(0.5, Float64(abs(r) + abs(p)), Float64(-q_m)); end return tmp end
q_m = N[Abs[q], $MachinePrecision] NOTE: p, r, and q_m should be sorted in increasing order before calling this function. code[p_, r_, q$95$m_] := If[LessEqual[N[Power[q$95$m, 2.0], $MachinePrecision], 20000000000000.0], N[(N[((-q$95$m) * q$95$m), $MachinePrecision] / r), $MachinePrecision], N[(0.5 * N[(N[Abs[r], $MachinePrecision] + N[Abs[p], $MachinePrecision]), $MachinePrecision] + (-q$95$m)), $MachinePrecision]]
\begin{array}{l}
q_m = \left|q\right|
\\
[p, r, q_m] = \mathsf{sort}([p, r, q_m])\\
\\
\begin{array}{l}
\mathbf{if}\;{q\_m}^{2} \leq 20000000000000:\\
\;\;\;\;\frac{\left(-q\_m\right) \cdot q\_m}{r}\\
\mathbf{else}:\\
\;\;\;\;\mathsf{fma}\left(0.5, \left|r\right| + \left|p\right|, -q\_m\right)\\
\end{array}
\end{array}
if (pow.f64 q #s(literal 2 binary64)) < 2e13Initial program 23.2%
Taylor expanded in r around inf
*-commutativeN/A
lower-*.f64N/A
Applied rewrites14.9%
Taylor expanded in r around 0
Applied rewrites37.1%
if 2e13 < (pow.f64 q #s(literal 2 binary64)) Initial program 27.7%
Taylor expanded in q around inf
*-commutativeN/A
lower-*.f64N/A
sub-negN/A
*-commutativeN/A
metadata-evalN/A
lower-fma.f64N/A
lower-/.f64N/A
+-commutativeN/A
lower-+.f64N/A
lower-fabs.f64N/A
lower-fabs.f6441.7
Applied rewrites41.7%
Taylor expanded in q around 0
Applied rewrites41.7%
q_m = (fabs.f64 q) NOTE: p, r, and q_m should be sorted in increasing order before calling this function. (FPCore (p r q_m) :precision binary64 (if (<= (pow q_m 2.0) 2e-232) (* (+ (+ p (fabs r)) (fabs p)) 0.5) (fma 0.5 (+ (fabs r) (fabs p)) (- q_m))))
q_m = fabs(q);
assert(p < r && r < q_m);
double code(double p, double r, double q_m) {
double tmp;
if (pow(q_m, 2.0) <= 2e-232) {
tmp = ((p + fabs(r)) + fabs(p)) * 0.5;
} else {
tmp = fma(0.5, (fabs(r) + fabs(p)), -q_m);
}
return tmp;
}
q_m = abs(q) p, r, q_m = sort([p, r, q_m]) function code(p, r, q_m) tmp = 0.0 if ((q_m ^ 2.0) <= 2e-232) tmp = Float64(Float64(Float64(p + abs(r)) + abs(p)) * 0.5); else tmp = fma(0.5, Float64(abs(r) + abs(p)), Float64(-q_m)); end return tmp end
q_m = N[Abs[q], $MachinePrecision] NOTE: p, r, and q_m should be sorted in increasing order before calling this function. code[p_, r_, q$95$m_] := If[LessEqual[N[Power[q$95$m, 2.0], $MachinePrecision], 2e-232], N[(N[(N[(p + N[Abs[r], $MachinePrecision]), $MachinePrecision] + N[Abs[p], $MachinePrecision]), $MachinePrecision] * 0.5), $MachinePrecision], N[(0.5 * N[(N[Abs[r], $MachinePrecision] + N[Abs[p], $MachinePrecision]), $MachinePrecision] + (-q$95$m)), $MachinePrecision]]
\begin{array}{l}
q_m = \left|q\right|
\\
[p, r, q_m] = \mathsf{sort}([p, r, q_m])\\
\\
\begin{array}{l}
\mathbf{if}\;{q\_m}^{2} \leq 2 \cdot 10^{-232}:\\
\;\;\;\;\left(\left(p + \left|r\right|\right) + \left|p\right|\right) \cdot 0.5\\
\mathbf{else}:\\
\;\;\;\;\mathsf{fma}\left(0.5, \left|r\right| + \left|p\right|, -q\_m\right)\\
\end{array}
\end{array}
if (pow.f64 q #s(literal 2 binary64)) < 2.00000000000000005e-232Initial program 28.9%
Taylor expanded in r around inf
*-commutativeN/A
lower-*.f64N/A
Applied rewrites17.6%
Taylor expanded in r around 0
Applied rewrites11.7%
if 2.00000000000000005e-232 < (pow.f64 q #s(literal 2 binary64)) Initial program 23.9%
Taylor expanded in q around inf
*-commutativeN/A
lower-*.f64N/A
sub-negN/A
*-commutativeN/A
metadata-evalN/A
lower-fma.f64N/A
lower-/.f64N/A
+-commutativeN/A
lower-+.f64N/A
lower-fabs.f64N/A
lower-fabs.f6433.2
Applied rewrites33.2%
Taylor expanded in q around 0
Applied rewrites33.2%
q_m = (fabs.f64 q) NOTE: p, r, and q_m should be sorted in increasing order before calling this function. (FPCore (p r q_m) :precision binary64 (- q_m))
q_m = fabs(q);
assert(p < r && r < q_m);
double code(double p, double r, double q_m) {
return -q_m;
}
q_m = abs(q)
NOTE: p, r, and q_m should be sorted in increasing order before calling this function.
real(8) function code(p, r, q_m)
real(8), intent (in) :: p
real(8), intent (in) :: r
real(8), intent (in) :: q_m
code = -q_m
end function
q_m = Math.abs(q);
assert p < r && r < q_m;
public static double code(double p, double r, double q_m) {
return -q_m;
}
q_m = math.fabs(q) [p, r, q_m] = sort([p, r, q_m]) def code(p, r, q_m): return -q_m
q_m = abs(q) p, r, q_m = sort([p, r, q_m]) function code(p, r, q_m) return Float64(-q_m) end
q_m = abs(q);
p, r, q_m = num2cell(sort([p, r, q_m])){:}
function tmp = code(p, r, q_m)
tmp = -q_m;
end
q_m = N[Abs[q], $MachinePrecision] NOTE: p, r, and q_m should be sorted in increasing order before calling this function. code[p_, r_, q$95$m_] := (-q$95$m)
\begin{array}{l}
q_m = \left|q\right|
\\
[p, r, q_m] = \mathsf{sort}([p, r, q_m])\\
\\
-q\_m
\end{array}
Initial program 25.4%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6424.5
Applied rewrites24.5%
herbie shell --seed 2024319
(FPCore (p r q)
:name "1/2(abs(p)+abs(r) - sqrt((p-r)^2 + 4q^2))"
:precision binary64
(* (/ 1.0 2.0) (- (+ (fabs p) (fabs r)) (sqrt (+ (pow (- p r) 2.0) (* 4.0 (pow q 2.0)))))))