
(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 5 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
(let* ((t_0 (fma (/ (+ (fabs p) (- (fabs r) r)) p) 0.5 0.5)))
(if (<= (pow q_m 2.0) 100.0)
(* t_0 p)
(if (<= (pow q_m 2.0) 5e+48)
(* (* q_m q_m) (fma (/ t_0 q_m) (/ p q_m) (/ 1.0 p)))
(- q_m)))))q_m = fabs(q);
assert(p < r && r < q_m);
double code(double p, double r, double q_m) {
double t_0 = fma(((fabs(p) + (fabs(r) - r)) / p), 0.5, 0.5);
double tmp;
if (pow(q_m, 2.0) <= 100.0) {
tmp = t_0 * p;
} else if (pow(q_m, 2.0) <= 5e+48) {
tmp = (q_m * q_m) * fma((t_0 / q_m), (p / q_m), (1.0 / p));
} 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) t_0 = fma(Float64(Float64(abs(p) + Float64(abs(r) - r)) / p), 0.5, 0.5) tmp = 0.0 if ((q_m ^ 2.0) <= 100.0) tmp = Float64(t_0 * p); elseif ((q_m ^ 2.0) <= 5e+48) tmp = Float64(Float64(q_m * q_m) * fma(Float64(t_0 / q_m), Float64(p / q_m), Float64(1.0 / p))); else tmp = 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_] := Block[{t$95$0 = N[(N[(N[(N[Abs[p], $MachinePrecision] + N[(N[Abs[r], $MachinePrecision] - r), $MachinePrecision]), $MachinePrecision] / p), $MachinePrecision] * 0.5 + 0.5), $MachinePrecision]}, If[LessEqual[N[Power[q$95$m, 2.0], $MachinePrecision], 100.0], N[(t$95$0 * p), $MachinePrecision], If[LessEqual[N[Power[q$95$m, 2.0], $MachinePrecision], 5e+48], N[(N[(q$95$m * q$95$m), $MachinePrecision] * N[(N[(t$95$0 / q$95$m), $MachinePrecision] * N[(p / q$95$m), $MachinePrecision] + N[(1.0 / p), $MachinePrecision]), $MachinePrecision]), $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}
t_0 := \mathsf{fma}\left(\frac{\left|p\right| + \left(\left|r\right| - r\right)}{p}, 0.5, 0.5\right)\\
\mathbf{if}\;{q\_m}^{2} \leq 100:\\
\;\;\;\;t\_0 \cdot p\\
\mathbf{elif}\;{q\_m}^{2} \leq 5 \cdot 10^{+48}:\\
\;\;\;\;\left(q\_m \cdot q\_m\right) \cdot \mathsf{fma}\left(\frac{t\_0}{q\_m}, \frac{p}{q\_m}, \frac{1}{p}\right)\\
\mathbf{else}:\\
\;\;\;\;-q\_m\\
\end{array}
\end{array}
if (pow.f64 q #s(literal 2 binary64)) < 100Initial program 24.0%
Taylor expanded in p around -inf
associate-*r*N/A
mul-1-negN/A
lower-*.f64N/A
lower-neg.f64N/A
sub-negN/A
metadata-evalN/A
+-commutativeN/A
mul-1-negN/A
unsub-negN/A
lower--.f64N/A
lower-/.f64N/A
Applied rewrites12.4%
Taylor expanded in p around inf
Applied rewrites30.0%
if 100 < (pow.f64 q #s(literal 2 binary64)) < 4.99999999999999973e48Initial program 12.1%
Taylor expanded in p around -inf
associate-*r*N/A
mul-1-negN/A
lower-*.f64N/A
lower-neg.f64N/A
sub-negN/A
metadata-evalN/A
+-commutativeN/A
mul-1-negN/A
unsub-negN/A
lower--.f64N/A
lower-/.f64N/A
Applied rewrites3.0%
Taylor expanded in q around inf
Applied rewrites21.9%
if 4.99999999999999973e48 < (pow.f64 q #s(literal 2 binary64)) Initial program 26.8%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6432.1
Applied rewrites32.1%
Final simplification30.3%
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-58) (* (fma (/ (+ (fabs p) (- (fabs r) r)) p) 0.5 0.5) p) (if (<= (pow q_m 2.0) 8e+31) (/ (* (- q_m) q_m) 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 (pow(q_m, 2.0) <= 2e-58) {
tmp = fma(((fabs(p) + (fabs(r) - r)) / p), 0.5, 0.5) * p;
} else if (pow(q_m, 2.0) <= 8e+31) {
tmp = (-q_m * q_m) / 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 ((q_m ^ 2.0) <= 2e-58) tmp = Float64(fma(Float64(Float64(abs(p) + Float64(abs(r) - r)) / p), 0.5, 0.5) * p); elseif ((q_m ^ 2.0) <= 8e+31) tmp = Float64(Float64(Float64(-q_m) * q_m) / r); else tmp = 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-58], N[(N[(N[(N[(N[Abs[p], $MachinePrecision] + N[(N[Abs[r], $MachinePrecision] - r), $MachinePrecision]), $MachinePrecision] / p), $MachinePrecision] * 0.5 + 0.5), $MachinePrecision] * p), $MachinePrecision], If[LessEqual[N[Power[q$95$m, 2.0], $MachinePrecision], 8e+31], N[(N[((-q$95$m) * q$95$m), $MachinePrecision] / 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}\;{q\_m}^{2} \leq 2 \cdot 10^{-58}:\\
\;\;\;\;\mathsf{fma}\left(\frac{\left|p\right| + \left(\left|r\right| - r\right)}{p}, 0.5, 0.5\right) \cdot p\\
\mathbf{elif}\;{q\_m}^{2} \leq 8 \cdot 10^{+31}:\\
\;\;\;\;\frac{\left(-q\_m\right) \cdot q\_m}{r}\\
\mathbf{else}:\\
\;\;\;\;-q\_m\\
\end{array}
\end{array}
if (pow.f64 q #s(literal 2 binary64)) < 2.0000000000000001e-58Initial program 24.2%
Taylor expanded in p around -inf
associate-*r*N/A
mul-1-negN/A
lower-*.f64N/A
lower-neg.f64N/A
sub-negN/A
metadata-evalN/A
+-commutativeN/A
mul-1-negN/A
unsub-negN/A
lower--.f64N/A
lower-/.f64N/A
Applied rewrites13.2%
Taylor expanded in p around inf
Applied rewrites31.4%
if 2.0000000000000001e-58 < (pow.f64 q #s(literal 2 binary64)) < 7.9999999999999997e31Initial program 10.8%
Taylor expanded in r around inf
Applied rewrites2.5%
Taylor expanded in r around 0
Applied rewrites3.3%
Taylor expanded in p around 0
Applied rewrites7.5%
if 7.9999999999999997e31 < (pow.f64 q #s(literal 2 binary64)) Initial program 27.0%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6432.0
Applied rewrites32.0%
Final simplification29.6%
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) 8e+31) (/ (* (- q_m) q_m) 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 (pow(q_m, 2.0) <= 8e+31) {
tmp = (-q_m * q_m) / 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 ((q_m ** 2.0d0) <= 8d+31) then
tmp = (-q_m * q_m) / 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 (Math.pow(q_m, 2.0) <= 8e+31) {
tmp = (-q_m * q_m) / 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 math.pow(q_m, 2.0) <= 8e+31: tmp = (-q_m * q_m) / 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 ((q_m ^ 2.0) <= 8e+31) tmp = Float64(Float64(Float64(-q_m) * q_m) / 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 ((q_m ^ 2.0) <= 8e+31)
tmp = (-q_m * q_m) / 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[Power[q$95$m, 2.0], $MachinePrecision], 8e+31], N[(N[((-q$95$m) * q$95$m), $MachinePrecision] / 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}\;{q\_m}^{2} \leq 8 \cdot 10^{+31}:\\
\;\;\;\;\frac{\left(-q\_m\right) \cdot q\_m}{r}\\
\mathbf{else}:\\
\;\;\;\;-q\_m\\
\end{array}
\end{array}
if (pow.f64 q #s(literal 2 binary64)) < 7.9999999999999997e31Initial program 21.9%
Taylor expanded in r around inf
Applied rewrites14.9%
Taylor expanded in r around 0
Applied rewrites24.7%
Taylor expanded in p around 0
Applied rewrites31.4%
if 7.9999999999999997e31 < (pow.f64 q #s(literal 2 binary64)) Initial program 27.0%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6432.0
Applied rewrites32.0%
Final simplification31.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) 5e+48) (/ (* q_m q_m) 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) <= 5e+48) {
tmp = (q_m * q_m) / p;
} 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 ** 2.0d0) <= 5d+48) then
tmp = (q_m * q_m) / p
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 (Math.pow(q_m, 2.0) <= 5e+48) {
tmp = (q_m * q_m) / p;
} 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 math.pow(q_m, 2.0) <= 5e+48: tmp = (q_m * q_m) / p 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 ^ 2.0) <= 5e+48) tmp = Float64(Float64(q_m * q_m) / p); 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 ^ 2.0) <= 5e+48)
tmp = (q_m * q_m) / p;
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[Power[q$95$m, 2.0], $MachinePrecision], 5e+48], N[(N[(q$95$m * q$95$m), $MachinePrecision] / p), $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}^{2} \leq 5 \cdot 10^{+48}:\\
\;\;\;\;\frac{q\_m \cdot q\_m}{p}\\
\mathbf{else}:\\
\;\;\;\;-q\_m\\
\end{array}
\end{array}
if (pow.f64 q #s(literal 2 binary64)) < 4.99999999999999973e48Initial program 22.3%
Taylor expanded in p around -inf
associate-*r*N/A
mul-1-negN/A
lower-*.f64N/A
lower-neg.f64N/A
sub-negN/A
metadata-evalN/A
+-commutativeN/A
mul-1-negN/A
unsub-negN/A
lower--.f64N/A
lower-/.f64N/A
Applied rewrites11.1%
Taylor expanded in p around 0
Applied rewrites34.2%
if 4.99999999999999973e48 < (pow.f64 q #s(literal 2 binary64)) Initial program 26.8%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6432.1
Applied rewrites32.1%
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 24.3%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6419.5
Applied rewrites19.5%
herbie shell --seed 2024332
(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)))))))