
(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 (if (<= (pow q_m 2.0) 1e-73) (fma (+ (fabs r) (+ (fabs p) p)) 0.5 (* -0.5 r)) (if (<= (pow q_m 2.0) 1e+141) (/ (* 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) <= 1e-73) {
tmp = fma((fabs(r) + (fabs(p) + p)), 0.5, (-0.5 * r));
} else if (pow(q_m, 2.0) <= 1e+141) {
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) <= 1e-73) tmp = fma(Float64(abs(r) + Float64(abs(p) + p)), 0.5, Float64(-0.5 * r)); elseif ((q_m ^ 2.0) <= 1e+141) tmp = Float64(Float64(q_m * q_m) / 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_] := If[LessEqual[N[Power[q$95$m, 2.0], $MachinePrecision], 1e-73], N[(N[(N[Abs[r], $MachinePrecision] + N[(N[Abs[p], $MachinePrecision] + p), $MachinePrecision]), $MachinePrecision] * 0.5 + N[(-0.5 * r), $MachinePrecision]), $MachinePrecision], If[LessEqual[N[Power[q$95$m, 2.0], $MachinePrecision], 1e+141], 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 10^{-73}:\\
\;\;\;\;\mathsf{fma}\left(\left|r\right| + \left(\left|p\right| + p\right), 0.5, -0.5 \cdot r\right)\\
\mathbf{elif}\;{q\_m}^{2} \leq 10^{+141}:\\
\;\;\;\;\frac{q\_m \cdot q\_m}{p}\\
\mathbf{else}:\\
\;\;\;\;-q\_m\\
\end{array}
\end{array}
if (pow.f64 q #s(literal 2 binary64)) < 9.99999999999999997e-74Initial program 28.3%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f648.6
Applied rewrites8.6%
Taylor expanded in r around inf
*-commutativeN/A
lower-*.f64N/A
Applied rewrites36.8%
Taylor expanded in r around 0
Applied rewrites16.5%
Applied rewrites36.9%
if 9.99999999999999997e-74 < (pow.f64 q #s(literal 2 binary64)) < 1.00000000000000002e141Initial program 23.8%
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
Applied rewrites3.6%
Taylor expanded in p around 0
Applied rewrites5.2%
Taylor expanded in p around inf
Applied rewrites11.1%
if 1.00000000000000002e141 < (pow.f64 q #s(literal 2 binary64)) Initial program 22.4%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6432.3
Applied rewrites32.3%
Final simplification31.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 (if (<= (pow q_m 2.0) 1e+141) (/ (* 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) <= 1e+141) {
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) <= 1d+141) 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) <= 1e+141) {
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) <= 1e+141: 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) <= 1e+141) 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) <= 1e+141)
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], 1e+141], 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 10^{+141}:\\
\;\;\;\;\frac{q\_m \cdot q\_m}{p}\\
\mathbf{else}:\\
\;\;\;\;-q\_m\\
\end{array}
\end{array}
if (pow.f64 q #s(literal 2 binary64)) < 1.00000000000000002e141Initial program 27.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
Applied rewrites14.6%
Taylor expanded in p around 0
Applied rewrites29.9%
Taylor expanded in p around inf
Applied rewrites33.5%
if 1.00000000000000002e141 < (pow.f64 q #s(literal 2 binary64)) Initial program 22.4%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6432.3
Applied rewrites32.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 (<= q_m 5.5e-173) (* (- (+ (fabs r) (fabs p)) r) 0.5) (if (<= q_m 7e-60) (* (+ (+ (fabs r) p) (fabs 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 <= 5.5e-173) {
tmp = ((fabs(r) + fabs(p)) - r) * 0.5;
} else if (q_m <= 7e-60) {
tmp = ((fabs(r) + p) + fabs(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 <= 5.5d-173) then
tmp = ((abs(r) + abs(p)) - r) * 0.5d0
else if (q_m <= 7d-60) then
tmp = ((abs(r) + p) + abs(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 <= 5.5e-173) {
tmp = ((Math.abs(r) + Math.abs(p)) - r) * 0.5;
} else if (q_m <= 7e-60) {
tmp = ((Math.abs(r) + p) + Math.abs(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 <= 5.5e-173: tmp = ((math.fabs(r) + math.fabs(p)) - r) * 0.5 elif q_m <= 7e-60: tmp = ((math.fabs(r) + p) + math.fabs(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 <= 5.5e-173) tmp = Float64(Float64(Float64(abs(r) + abs(p)) - r) * 0.5); elseif (q_m <= 7e-60) tmp = Float64(Float64(Float64(abs(r) + p) + abs(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 <= 5.5e-173)
tmp = ((abs(r) + abs(p)) - r) * 0.5;
elseif (q_m <= 7e-60)
tmp = ((abs(r) + p) + abs(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, 5.5e-173], N[(N[(N[(N[Abs[r], $MachinePrecision] + N[Abs[p], $MachinePrecision]), $MachinePrecision] - r), $MachinePrecision] * 0.5), $MachinePrecision], If[LessEqual[q$95$m, 7e-60], N[(N[(N[(N[Abs[r], $MachinePrecision] + p), $MachinePrecision] + N[Abs[p], $MachinePrecision]), $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 5.5 \cdot 10^{-173}:\\
\;\;\;\;\left(\left(\left|r\right| + \left|p\right|\right) - r\right) \cdot 0.5\\
\mathbf{elif}\;q\_m \leq 7 \cdot 10^{-60}:\\
\;\;\;\;\left(\left(\left|r\right| + p\right) + \left|p\right|\right) \cdot 0.5\\
\mathbf{else}:\\
\;\;\;\;-q\_m\\
\end{array}
\end{array}
if q < 5.50000000000000022e-173Initial program 27.3%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f643.5
Applied rewrites3.5%
Taylor expanded in r around inf
*-commutativeN/A
lower-*.f64N/A
Applied rewrites25.0%
Taylor expanded in r around 0
Applied rewrites11.3%
Taylor expanded in p around 0
Applied rewrites9.1%
if 5.50000000000000022e-173 < q < 6.99999999999999952e-60Initial program 28.2%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6415.7
Applied rewrites15.7%
Taylor expanded in r around inf
*-commutativeN/A
lower-*.f64N/A
Applied rewrites30.4%
Taylor expanded in r around 0
Applied rewrites10.2%
if 6.99999999999999952e-60 < q Initial program 20.1%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6449.6
Applied rewrites49.6%
Final simplification20.9%
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 6.8e-69) (* (- (+ (fabs r) (fabs p)) r) 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 <= 6.8e-69) {
tmp = ((fabs(r) + fabs(p)) - r) * 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 <= 6.8d-69) then
tmp = ((abs(r) + abs(p)) - r) * 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 <= 6.8e-69) {
tmp = ((Math.abs(r) + Math.abs(p)) - r) * 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 <= 6.8e-69: tmp = ((math.fabs(r) + math.fabs(p)) - r) * 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 <= 6.8e-69) tmp = Float64(Float64(Float64(abs(r) + abs(p)) - r) * 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 <= 6.8e-69)
tmp = ((abs(r) + abs(p)) - r) * 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, 6.8e-69], N[(N[(N[(N[Abs[r], $MachinePrecision] + N[Abs[p], $MachinePrecision]), $MachinePrecision] - r), $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 6.8 \cdot 10^{-69}:\\
\;\;\;\;\left(\left(\left|r\right| + \left|p\right|\right) - r\right) \cdot 0.5\\
\mathbf{else}:\\
\;\;\;\;-q\_m\\
\end{array}
\end{array}
if q < 6.80000000000000016e-69Initial program 27.6%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f645.5
Applied rewrites5.5%
Taylor expanded in r around inf
*-commutativeN/A
lower-*.f64N/A
Applied rewrites26.0%
Taylor expanded in r around 0
Applied rewrites11.9%
Taylor expanded in p around 0
Applied rewrites9.8%
if 6.80000000000000016e-69 < q Initial program 19.9%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6449.0
Applied rewrites49.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 (- 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.f6418.2
Applied rewrites18.2%
herbie shell --seed 2024284
(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)))))))