
(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}
NOTE: p, r, and q should be sorted in increasing order before calling this function. (FPCore (p r q) :precision binary64 (if (<= (pow q 2.0) 1e-43) (* (+ (- (fabs r) r) (+ (fabs p) p)) 0.5) (* (- (+ (fabs r) (fabs p)) (hypot (* -2.0 q) (- p r))) 0.5)))
assert(p < r && r < q);
double code(double p, double r, double q) {
double tmp;
if (pow(q, 2.0) <= 1e-43) {
tmp = ((fabs(r) - r) + (fabs(p) + p)) * 0.5;
} else {
tmp = ((fabs(r) + fabs(p)) - hypot((-2.0 * q), (p - r))) * 0.5;
}
return tmp;
}
assert p < r && r < q;
public static double code(double p, double r, double q) {
double tmp;
if (Math.pow(q, 2.0) <= 1e-43) {
tmp = ((Math.abs(r) - r) + (Math.abs(p) + p)) * 0.5;
} else {
tmp = ((Math.abs(r) + Math.abs(p)) - Math.hypot((-2.0 * q), (p - r))) * 0.5;
}
return tmp;
}
[p, r, q] = sort([p, r, q]) def code(p, r, q): tmp = 0 if math.pow(q, 2.0) <= 1e-43: tmp = ((math.fabs(r) - r) + (math.fabs(p) + p)) * 0.5 else: tmp = ((math.fabs(r) + math.fabs(p)) - math.hypot((-2.0 * q), (p - r))) * 0.5 return tmp
p, r, q = sort([p, r, q]) function code(p, r, q) tmp = 0.0 if ((q ^ 2.0) <= 1e-43) tmp = Float64(Float64(Float64(abs(r) - r) + Float64(abs(p) + p)) * 0.5); else tmp = Float64(Float64(Float64(abs(r) + abs(p)) - hypot(Float64(-2.0 * q), Float64(p - r))) * 0.5); end return tmp end
p, r, q = num2cell(sort([p, r, q])){:}
function tmp_2 = code(p, r, q)
tmp = 0.0;
if ((q ^ 2.0) <= 1e-43)
tmp = ((abs(r) - r) + (abs(p) + p)) * 0.5;
else
tmp = ((abs(r) + abs(p)) - hypot((-2.0 * q), (p - r))) * 0.5;
end
tmp_2 = tmp;
end
NOTE: p, r, and q should be sorted in increasing order before calling this function. code[p_, r_, q_] := If[LessEqual[N[Power[q, 2.0], $MachinePrecision], 1e-43], N[(N[(N[(N[Abs[r], $MachinePrecision] - r), $MachinePrecision] + N[(N[Abs[p], $MachinePrecision] + p), $MachinePrecision]), $MachinePrecision] * 0.5), $MachinePrecision], N[(N[(N[(N[Abs[r], $MachinePrecision] + N[Abs[p], $MachinePrecision]), $MachinePrecision] - N[Sqrt[N[(-2.0 * q), $MachinePrecision] ^ 2 + N[(p - r), $MachinePrecision] ^ 2], $MachinePrecision]), $MachinePrecision] * 0.5), $MachinePrecision]]
\begin{array}{l}
[p, r, q] = \mathsf{sort}([p, r, q])\\
\\
\begin{array}{l}
\mathbf{if}\;{q}^{2} \leq 10^{-43}:\\
\;\;\;\;\left(\left(\left|r\right| - r\right) + \left(\left|p\right| + p\right)\right) \cdot 0.5\\
\mathbf{else}:\\
\;\;\;\;\left(\left(\left|r\right| + \left|p\right|\right) - \mathsf{hypot}\left(-2 \cdot q, p - r\right)\right) \cdot 0.5\\
\end{array}
\end{array}
if (pow.f64 q #s(literal 2 binary64)) < 1.00000000000000008e-43Initial program 27.7%
Taylor expanded in p around 0
distribute-lft-outN/A
*-commutativeN/A
lower-*.f64N/A
Applied rewrites11.1%
Taylor expanded in r around 0
Applied rewrites8.1%
Applied rewrites8.1%
Taylor expanded in q around 0
Applied rewrites41.2%
if 1.00000000000000008e-43 < (pow.f64 q #s(literal 2 binary64)) Initial program 23.4%
lift--.f64N/A
sub-negN/A
+-commutativeN/A
lift-+.f64N/A
flip-+N/A
lift-fabs.f64N/A
lift-fabs.f64N/A
sqr-absN/A
lift-fabs.f64N/A
lift-fabs.f64N/A
sqr-absN/A
div-subN/A
associate-+r-N/A
Applied rewrites21.8%
Applied rewrites59.1%
Final simplification51.3%
NOTE: p, r, and q should be sorted in increasing order before calling this function. (FPCore (p r q) :precision binary64 (if (<= (pow q 2.0) 1e-43) (* (+ (- (fabs r) r) (+ (fabs p) p)) 0.5) (- q)))
assert(p < r && r < q);
double code(double p, double r, double q) {
double tmp;
if (pow(q, 2.0) <= 1e-43) {
tmp = ((fabs(r) - r) + (fabs(p) + p)) * 0.5;
} else {
tmp = -q;
}
return tmp;
}
NOTE: p, r, and q should be sorted in increasing order before calling this function.
real(8) function code(p, r, q)
real(8), intent (in) :: p
real(8), intent (in) :: r
real(8), intent (in) :: q
real(8) :: tmp
if ((q ** 2.0d0) <= 1d-43) then
tmp = ((abs(r) - r) + (abs(p) + p)) * 0.5d0
else
tmp = -q
end if
code = tmp
end function
assert p < r && r < q;
public static double code(double p, double r, double q) {
double tmp;
if (Math.pow(q, 2.0) <= 1e-43) {
tmp = ((Math.abs(r) - r) + (Math.abs(p) + p)) * 0.5;
} else {
tmp = -q;
}
return tmp;
}
[p, r, q] = sort([p, r, q]) def code(p, r, q): tmp = 0 if math.pow(q, 2.0) <= 1e-43: tmp = ((math.fabs(r) - r) + (math.fabs(p) + p)) * 0.5 else: tmp = -q return tmp
p, r, q = sort([p, r, q]) function code(p, r, q) tmp = 0.0 if ((q ^ 2.0) <= 1e-43) tmp = Float64(Float64(Float64(abs(r) - r) + Float64(abs(p) + p)) * 0.5); else tmp = Float64(-q); end return tmp end
p, r, q = num2cell(sort([p, r, q])){:}
function tmp_2 = code(p, r, q)
tmp = 0.0;
if ((q ^ 2.0) <= 1e-43)
tmp = ((abs(r) - r) + (abs(p) + p)) * 0.5;
else
tmp = -q;
end
tmp_2 = tmp;
end
NOTE: p, r, and q should be sorted in increasing order before calling this function. code[p_, r_, q_] := If[LessEqual[N[Power[q, 2.0], $MachinePrecision], 1e-43], N[(N[(N[(N[Abs[r], $MachinePrecision] - r), $MachinePrecision] + N[(N[Abs[p], $MachinePrecision] + p), $MachinePrecision]), $MachinePrecision] * 0.5), $MachinePrecision], (-q)]
\begin{array}{l}
[p, r, q] = \mathsf{sort}([p, r, q])\\
\\
\begin{array}{l}
\mathbf{if}\;{q}^{2} \leq 10^{-43}:\\
\;\;\;\;\left(\left(\left|r\right| - r\right) + \left(\left|p\right| + p\right)\right) \cdot 0.5\\
\mathbf{else}:\\
\;\;\;\;-q\\
\end{array}
\end{array}
if (pow.f64 q #s(literal 2 binary64)) < 1.00000000000000008e-43Initial program 27.7%
Taylor expanded in p around 0
distribute-lft-outN/A
*-commutativeN/A
lower-*.f64N/A
Applied rewrites11.1%
Taylor expanded in r around 0
Applied rewrites8.1%
Applied rewrites8.1%
Taylor expanded in q around 0
Applied rewrites41.2%
if 1.00000000000000008e-43 < (pow.f64 q #s(literal 2 binary64)) Initial program 23.4%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6429.8
Applied rewrites29.8%
Final simplification34.8%
NOTE: p, r, and q should be sorted in increasing order before calling this function. (FPCore (p r q) :precision binary64 (if (<= (pow q 2.0) 20000000.0) (/ (* (- q) q) r) (- q)))
assert(p < r && r < q);
double code(double p, double r, double q) {
double tmp;
if (pow(q, 2.0) <= 20000000.0) {
tmp = (-q * q) / r;
} else {
tmp = -q;
}
return tmp;
}
NOTE: p, r, and q should be sorted in increasing order before calling this function.
real(8) function code(p, r, q)
real(8), intent (in) :: p
real(8), intent (in) :: r
real(8), intent (in) :: q
real(8) :: tmp
if ((q ** 2.0d0) <= 20000000.0d0) then
tmp = (-q * q) / r
else
tmp = -q
end if
code = tmp
end function
assert p < r && r < q;
public static double code(double p, double r, double q) {
double tmp;
if (Math.pow(q, 2.0) <= 20000000.0) {
tmp = (-q * q) / r;
} else {
tmp = -q;
}
return tmp;
}
[p, r, q] = sort([p, r, q]) def code(p, r, q): tmp = 0 if math.pow(q, 2.0) <= 20000000.0: tmp = (-q * q) / r else: tmp = -q return tmp
p, r, q = sort([p, r, q]) function code(p, r, q) tmp = 0.0 if ((q ^ 2.0) <= 20000000.0) tmp = Float64(Float64(Float64(-q) * q) / r); else tmp = Float64(-q); end return tmp end
p, r, q = num2cell(sort([p, r, q])){:}
function tmp_2 = code(p, r, q)
tmp = 0.0;
if ((q ^ 2.0) <= 20000000.0)
tmp = (-q * q) / r;
else
tmp = -q;
end
tmp_2 = tmp;
end
NOTE: p, r, and q should be sorted in increasing order before calling this function. code[p_, r_, q_] := If[LessEqual[N[Power[q, 2.0], $MachinePrecision], 20000000.0], N[(N[((-q) * q), $MachinePrecision] / r), $MachinePrecision], (-q)]
\begin{array}{l}
[p, r, q] = \mathsf{sort}([p, r, q])\\
\\
\begin{array}{l}
\mathbf{if}\;{q}^{2} \leq 20000000:\\
\;\;\;\;\frac{\left(-q\right) \cdot q}{r}\\
\mathbf{else}:\\
\;\;\;\;-q\\
\end{array}
\end{array}
if (pow.f64 q #s(literal 2 binary64)) < 2e7Initial program 28.6%
Applied rewrites7.5%
Taylor expanded in r around inf
distribute-lft-outN/A
lower-*.f64N/A
lower-+.f64N/A
lower-/.f64N/A
lower--.f64N/A
lower-fabs.f64N/A
*-commutativeN/A
lower-*.f64N/A
*-commutativeN/A
lower-fma.f64N/A
unpow2N/A
lower-*.f64N/A
unpow2N/A
lower-*.f6421.8
Applied rewrites21.8%
Taylor expanded in q around inf
Applied rewrites38.6%
if 2e7 < (pow.f64 q #s(literal 2 binary64)) Initial program 22.2%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6430.6
Applied rewrites30.6%
NOTE: p, r, and q should be sorted in increasing order before calling this function. (FPCore (p r q) :precision binary64 (- q))
assert(p < r && r < q);
double code(double p, double r, double q) {
return -q;
}
NOTE: p, r, and q should be sorted in increasing order before calling this function.
real(8) function code(p, r, q)
real(8), intent (in) :: p
real(8), intent (in) :: r
real(8), intent (in) :: q
code = -q
end function
assert p < r && r < q;
public static double code(double p, double r, double q) {
return -q;
}
[p, r, q] = sort([p, r, q]) def code(p, r, q): return -q
p, r, q = sort([p, r, q]) function code(p, r, q) return Float64(-q) end
p, r, q = num2cell(sort([p, r, q])){:}
function tmp = code(p, r, q)
tmp = -q;
end
NOTE: p, r, and q should be sorted in increasing order before calling this function. code[p_, r_, q_] := (-q)
\begin{array}{l}
[p, r, q] = \mathsf{sort}([p, r, q])\\
\\
-q
\end{array}
Initial program 25.3%
Taylor expanded in q around inf
mul-1-negN/A
lower-neg.f6420.8
Applied rewrites20.8%
herbie shell --seed 2024285
(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)))))))