
(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 (<= r -2.3e-160)
(* (- (* (+ (+ (fabs p) p) (fabs r)) (/ 0.5 r)) 0.5) r)
(if (<= r 5.5e-62)
(- q_m)
(*
(pow 2.0 -1.0)
(fma (/ (+ p r) r) (* (* -2.0 q_m) (/ q_m r)) (+ p (fabs p)))))))q_m = fabs(q);
assert(p < r && r < q_m);
double code(double p, double r, double q_m) {
double tmp;
if (r <= -2.3e-160) {
tmp = ((((fabs(p) + p) + fabs(r)) * (0.5 / r)) - 0.5) * r;
} else if (r <= 5.5e-62) {
tmp = -q_m;
} else {
tmp = pow(2.0, -1.0) * fma(((p + r) / r), ((-2.0 * q_m) * (q_m / r)), (p + fabs(p)));
}
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 (r <= -2.3e-160) tmp = Float64(Float64(Float64(Float64(Float64(abs(p) + p) + abs(r)) * Float64(0.5 / r)) - 0.5) * r); elseif (r <= 5.5e-62) tmp = Float64(-q_m); else tmp = Float64((2.0 ^ -1.0) * fma(Float64(Float64(p + r) / r), Float64(Float64(-2.0 * q_m) * Float64(q_m / r)), Float64(p + abs(p)))); 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[r, -2.3e-160], N[(N[(N[(N[(N[(N[Abs[p], $MachinePrecision] + p), $MachinePrecision] + N[Abs[r], $MachinePrecision]), $MachinePrecision] * N[(0.5 / r), $MachinePrecision]), $MachinePrecision] - 0.5), $MachinePrecision] * r), $MachinePrecision], If[LessEqual[r, 5.5e-62], (-q$95$m), N[(N[Power[2.0, -1.0], $MachinePrecision] * N[(N[(N[(p + r), $MachinePrecision] / r), $MachinePrecision] * N[(N[(-2.0 * q$95$m), $MachinePrecision] * N[(q$95$m / r), $MachinePrecision]), $MachinePrecision] + N[(p + N[Abs[p], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]
\begin{array}{l}
q_m = \left|q\right|
\\
[p, r, q_m] = \mathsf{sort}([p, r, q_m])\\
\\
\begin{array}{l}
\mathbf{if}\;r \leq -2.3 \cdot 10^{-160}:\\
\;\;\;\;\left(\left(\left(\left|p\right| + p\right) + \left|r\right|\right) \cdot \frac{0.5}{r} - 0.5\right) \cdot r\\
\mathbf{elif}\;r \leq 5.5 \cdot 10^{-62}:\\
\;\;\;\;-q\_m\\
\mathbf{else}:\\
\;\;\;\;{2}^{-1} \cdot \mathsf{fma}\left(\frac{p + r}{r}, \left(-2 \cdot q\_m\right) \cdot \frac{q\_m}{r}, p + \left|p\right|\right)\\
\end{array}
\end{array}
if r < -2.29999999999999985e-160Initial program 15.7%
Taylor expanded in r around inf
*-commutativeN/A
lower-*.f64N/A
Applied rewrites4.2%
Applied rewrites13.6%
if -2.29999999999999985e-160 < r < 5.50000000000000022e-62Initial program 34.7%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6420.4
Applied rewrites20.4%
if 5.50000000000000022e-62 < r Initial program 15.7%
Taylor expanded in r around inf
Applied rewrites16.9%
Taylor expanded in r around 0
Applied rewrites26.7%
Applied rewrites31.7%
Applied rewrites59.0%
Final simplification31.0%
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 (<= (* 4.0 (pow q_m 2.0)) 2e-108) (* 0.0 r) (- q_m)))
q_m = fabs(q);
assert(p < r && r < q_m);
double code(double p, double r, double q_m) {
double tmp;
if ((4.0 * pow(q_m, 2.0)) <= 2e-108) {
tmp = 0.0 * r;
} else {
tmp = -q_m;
}
return tmp;
}
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
real(8) :: tmp
if ((4.0d0 * (q_m ** 2.0d0)) <= 2d-108) then
tmp = 0.0d0 * r
else
tmp = -q_m
end if
code = tmp
end function
q_m = Math.abs(q);
assert p < r && r < q_m;
public static double code(double p, double r, double q_m) {
double tmp;
if ((4.0 * Math.pow(q_m, 2.0)) <= 2e-108) {
tmp = 0.0 * r;
} else {
tmp = -q_m;
}
return tmp;
}
q_m = math.fabs(q) [p, r, q_m] = sort([p, r, q_m]) def code(p, r, q_m): tmp = 0 if (4.0 * math.pow(q_m, 2.0)) <= 2e-108: tmp = 0.0 * r else: tmp = -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 (Float64(4.0 * (q_m ^ 2.0)) <= 2e-108) tmp = Float64(0.0 * r); else tmp = Float64(-q_m); end return tmp end
q_m = abs(q);
p, r, q_m = num2cell(sort([p, r, q_m])){:}
function tmp_2 = code(p, r, q_m)
tmp = 0.0;
if ((4.0 * (q_m ^ 2.0)) <= 2e-108)
tmp = 0.0 * r;
else
tmp = -q_m;
end
tmp_2 = 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[(4.0 * N[Power[q$95$m, 2.0], $MachinePrecision]), $MachinePrecision], 2e-108], N[(0.0 * r), $MachinePrecision], (-q$95$m)]
\begin{array}{l}
q_m = \left|q\right|
\\
[p, r, q_m] = \mathsf{sort}([p, r, q_m])\\
\\
\begin{array}{l}
\mathbf{if}\;4 \cdot {q\_m}^{2} \leq 2 \cdot 10^{-108}:\\
\;\;\;\;0 \cdot r\\
\mathbf{else}:\\
\;\;\;\;-q\_m\\
\end{array}
\end{array}
if (*.f64 #s(literal 4 binary64) (pow.f64 q #s(literal 2 binary64))) < 2.00000000000000008e-108Initial program 20.7%
Taylor expanded in r around inf
*-commutativeN/A
lower-*.f64N/A
Applied rewrites11.5%
Applied rewrites12.8%
Taylor expanded in r around -inf
Applied rewrites37.4%
if 2.00000000000000008e-108 < (*.f64 #s(literal 4 binary64) (pow.f64 q #s(literal 2 binary64))) Initial program 22.6%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6423.4
Applied rewrites23.4%
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 (<= q_m 1.4e-35) (* (+ (fabs p) p) 0.5) (- q_m)))
q_m = fabs(q);
assert(p < r && r < q_m);
double code(double p, double r, double q_m) {
double tmp;
if (q_m <= 1.4e-35) {
tmp = (fabs(p) + p) * 0.5;
} else {
tmp = -q_m;
}
return tmp;
}
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
real(8) :: tmp
if (q_m <= 1.4d-35) then
tmp = (abs(p) + p) * 0.5d0
else
tmp = -q_m
end if
code = tmp
end function
q_m = Math.abs(q);
assert p < r && r < q_m;
public static double code(double p, double r, double q_m) {
double tmp;
if (q_m <= 1.4e-35) {
tmp = (Math.abs(p) + p) * 0.5;
} else {
tmp = -q_m;
}
return tmp;
}
q_m = math.fabs(q) [p, r, q_m] = sort([p, r, q_m]) def code(p, r, q_m): tmp = 0 if q_m <= 1.4e-35: tmp = (math.fabs(p) + p) * 0.5 else: tmp = -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 <= 1.4e-35) tmp = Float64(Float64(abs(p) + p) * 0.5); else tmp = Float64(-q_m); end return tmp end
q_m = abs(q);
p, r, q_m = num2cell(sort([p, r, q_m])){:}
function tmp_2 = code(p, r, q_m)
tmp = 0.0;
if (q_m <= 1.4e-35)
tmp = (abs(p) + p) * 0.5;
else
tmp = -q_m;
end
tmp_2 = 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[q$95$m, 1.4e-35], N[(N[(N[Abs[p], $MachinePrecision] + p), $MachinePrecision] * 0.5), $MachinePrecision], (-q$95$m)]
\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 \leq 1.4 \cdot 10^{-35}:\\
\;\;\;\;\left(\left|p\right| + p\right) \cdot 0.5\\
\mathbf{else}:\\
\;\;\;\;-q\_m\\
\end{array}
\end{array}
if q < 1.4e-35Initial program 22.3%
Taylor expanded in r around inf
*-commutativeN/A
lower-*.f64N/A
Applied rewrites7.6%
Applied rewrites7.7%
Taylor expanded in r around 0
Applied rewrites21.0%
if 1.4e-35 < q Initial program 20.4%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6452.8
Applied rewrites52.8%
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 21.8%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6416.4
Applied rewrites16.4%
herbie shell --seed 2024337
(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)))))))