
(FPCore (x y z) :precision binary64 (/ (- (+ (* x x) (* y y)) (* z z)) (* y 2.0)))
double code(double x, double y, double z) {
return (((x * x) + (y * y)) - (z * z)) / (y * 2.0);
}
real(8) function code(x, y, z)
real(8), intent (in) :: x
real(8), intent (in) :: y
real(8), intent (in) :: z
code = (((x * x) + (y * y)) - (z * z)) / (y * 2.0d0)
end function
public static double code(double x, double y, double z) {
return (((x * x) + (y * y)) - (z * z)) / (y * 2.0);
}
def code(x, y, z): return (((x * x) + (y * y)) - (z * z)) / (y * 2.0)
function code(x, y, z) return Float64(Float64(Float64(Float64(x * x) + Float64(y * y)) - Float64(z * z)) / Float64(y * 2.0)) end
function tmp = code(x, y, z) tmp = (((x * x) + (y * y)) - (z * z)) / (y * 2.0); end
code[x_, y_, z_] := N[(N[(N[(N[(x * x), $MachinePrecision] + N[(y * y), $MachinePrecision]), $MachinePrecision] - N[(z * z), $MachinePrecision]), $MachinePrecision] / N[(y * 2.0), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 4 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (x y z) :precision binary64 (/ (- (+ (* x x) (* y y)) (* z z)) (* y 2.0)))
double code(double x, double y, double z) {
return (((x * x) + (y * y)) - (z * z)) / (y * 2.0);
}
real(8) function code(x, y, z)
real(8), intent (in) :: x
real(8), intent (in) :: y
real(8), intent (in) :: z
code = (((x * x) + (y * y)) - (z * z)) / (y * 2.0d0)
end function
public static double code(double x, double y, double z) {
return (((x * x) + (y * y)) - (z * z)) / (y * 2.0);
}
def code(x, y, z): return (((x * x) + (y * y)) - (z * z)) / (y * 2.0)
function code(x, y, z) return Float64(Float64(Float64(Float64(x * x) + Float64(y * y)) - Float64(z * z)) / Float64(y * 2.0)) end
function tmp = code(x, y, z) tmp = (((x * x) + (y * y)) - (z * z)) / (y * 2.0); end
code[x_, y_, z_] := N[(N[(N[(N[(x * x), $MachinePrecision] + N[(y * y), $MachinePrecision]), $MachinePrecision] - N[(z * z), $MachinePrecision]), $MachinePrecision] / N[(y * 2.0), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\frac{\left(x \cdot x + y \cdot y\right) - z \cdot z}{y \cdot 2}
\end{array}
y\_m = (fabs.f64 y)
y\_s = (copysign.f64 #s(literal 1 binary64) y)
(FPCore (y_s x y_m z)
:precision binary64
(*
y_s
(if (<= (/ (- (+ (* x x) (* y_m y_m)) (* z z)) (* y_m 2.0)) -5e-98)
(/ 1.0 (* (/ 1.0 z) (* (/ 1.0 z) (* y_m -2.0))))
(* 0.5 (fma x (/ x y_m) y_m)))))y\_m = fabs(y);
y\_s = copysign(1.0, y);
double code(double y_s, double x, double y_m, double z) {
double tmp;
if (((((x * x) + (y_m * y_m)) - (z * z)) / (y_m * 2.0)) <= -5e-98) {
tmp = 1.0 / ((1.0 / z) * ((1.0 / z) * (y_m * -2.0)));
} else {
tmp = 0.5 * fma(x, (x / y_m), y_m);
}
return y_s * tmp;
}
y\_m = abs(y) y\_s = copysign(1.0, y) function code(y_s, x, y_m, z) tmp = 0.0 if (Float64(Float64(Float64(Float64(x * x) + Float64(y_m * y_m)) - Float64(z * z)) / Float64(y_m * 2.0)) <= -5e-98) tmp = Float64(1.0 / Float64(Float64(1.0 / z) * Float64(Float64(1.0 / z) * Float64(y_m * -2.0)))); else tmp = Float64(0.5 * fma(x, Float64(x / y_m), y_m)); end return Float64(y_s * tmp) end
y\_m = N[Abs[y], $MachinePrecision]
y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
code[y$95$s_, x_, y$95$m_, z_] := N[(y$95$s * If[LessEqual[N[(N[(N[(N[(x * x), $MachinePrecision] + N[(y$95$m * y$95$m), $MachinePrecision]), $MachinePrecision] - N[(z * z), $MachinePrecision]), $MachinePrecision] / N[(y$95$m * 2.0), $MachinePrecision]), $MachinePrecision], -5e-98], N[(1.0 / N[(N[(1.0 / z), $MachinePrecision] * N[(N[(1.0 / z), $MachinePrecision] * N[(y$95$m * -2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(0.5 * N[(x * N[(x / y$95$m), $MachinePrecision] + y$95$m), $MachinePrecision]), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
y\_m = \left|y\right|
\\
y\_s = \mathsf{copysign}\left(1, y\right)
\\
y\_s \cdot \begin{array}{l}
\mathbf{if}\;\frac{\left(x \cdot x + y\_m \cdot y\_m\right) - z \cdot z}{y\_m \cdot 2} \leq -5 \cdot 10^{-98}:\\
\;\;\;\;\frac{1}{\frac{1}{z} \cdot \left(\frac{1}{z} \cdot \left(y\_m \cdot -2\right)\right)}\\
\mathbf{else}:\\
\;\;\;\;0.5 \cdot \mathsf{fma}\left(x, \frac{x}{y\_m}, y\_m\right)\\
\end{array}
\end{array}
if (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) < -5.00000000000000018e-98Initial program 75.5%
Taylor expanded in z around inf 27.7%
associate-*r/27.7%
metadata-eval27.7%
distribute-lft-neg-in27.7%
*-commutative27.7%
distribute-neg-frac27.7%
associate-*r/27.7%
distribute-rgt-neg-in27.7%
distribute-neg-frac27.7%
metadata-eval27.7%
Simplified27.7%
associate-*r/27.7%
clear-num27.7%
Applied egg-rr27.7%
add-cube-cbrt27.5%
pow327.5%
*-un-lft-identity27.5%
times-frac27.6%
pow-flip27.6%
metadata-eval27.6%
Applied egg-rr27.6%
rem-cube-cbrt27.7%
sqr-pow27.7%
associate-*l*28.5%
metadata-eval28.5%
unpow-128.5%
metadata-eval28.5%
unpow-128.5%
div-inv28.5%
metadata-eval28.5%
Applied egg-rr28.5%
if -5.00000000000000018e-98 < (/.f64 (-.f64 (+.f64 (*.f64 x x) (*.f64 y y)) (*.f64 z z)) (*.f64 y #s(literal 2 binary64))) Initial program 59.2%
remove-double-neg59.2%
distribute-lft-neg-out59.2%
distribute-frac-neg259.2%
distribute-frac-neg59.2%
neg-mul-159.2%
distribute-lft-neg-out59.2%
*-commutative59.2%
distribute-lft-neg-in59.2%
times-frac59.2%
metadata-eval59.2%
metadata-eval59.2%
associate--l+59.2%
fma-define64.0%
Simplified64.0%
Taylor expanded in z around inf 49.7%
associate--l+49.7%
unpow249.7%
times-frac59.1%
fma-neg59.1%
distribute-neg-frac59.1%
metadata-eval59.1%
Simplified59.1%
*-un-lft-identity59.1%
unpow259.1%
times-frac59.5%
Applied egg-rr59.5%
Taylor expanded in z around 0 58.4%
+-commutative58.4%
unpow258.4%
associate-/l*65.3%
fma-define65.3%
Simplified65.3%
Final simplification49.5%
y\_m = (fabs.f64 y)
y\_s = (copysign.f64 #s(literal 1 binary64) y)
(FPCore (y_s x y_m z)
:precision binary64
(*
y_s
(if (or (<= z 1.85e+49) (and (not (<= z 1.5e+71)) (<= z 2.55e+124)))
(* y_m 0.5)
(/ 1.0 (* (/ 1.0 z) (* (/ 1.0 z) (* y_m -2.0)))))))y\_m = fabs(y);
y\_s = copysign(1.0, y);
double code(double y_s, double x, double y_m, double z) {
double tmp;
if ((z <= 1.85e+49) || (!(z <= 1.5e+71) && (z <= 2.55e+124))) {
tmp = y_m * 0.5;
} else {
tmp = 1.0 / ((1.0 / z) * ((1.0 / z) * (y_m * -2.0)));
}
return y_s * tmp;
}
y\_m = abs(y)
y\_s = copysign(1.0d0, y)
real(8) function code(y_s, x, y_m, z)
real(8), intent (in) :: y_s
real(8), intent (in) :: x
real(8), intent (in) :: y_m
real(8), intent (in) :: z
real(8) :: tmp
if ((z <= 1.85d+49) .or. (.not. (z <= 1.5d+71)) .and. (z <= 2.55d+124)) then
tmp = y_m * 0.5d0
else
tmp = 1.0d0 / ((1.0d0 / z) * ((1.0d0 / z) * (y_m * (-2.0d0))))
end if
code = y_s * tmp
end function
y\_m = Math.abs(y);
y\_s = Math.copySign(1.0, y);
public static double code(double y_s, double x, double y_m, double z) {
double tmp;
if ((z <= 1.85e+49) || (!(z <= 1.5e+71) && (z <= 2.55e+124))) {
tmp = y_m * 0.5;
} else {
tmp = 1.0 / ((1.0 / z) * ((1.0 / z) * (y_m * -2.0)));
}
return y_s * tmp;
}
y\_m = math.fabs(y) y\_s = math.copysign(1.0, y) def code(y_s, x, y_m, z): tmp = 0 if (z <= 1.85e+49) or (not (z <= 1.5e+71) and (z <= 2.55e+124)): tmp = y_m * 0.5 else: tmp = 1.0 / ((1.0 / z) * ((1.0 / z) * (y_m * -2.0))) return y_s * tmp
y\_m = abs(y) y\_s = copysign(1.0, y) function code(y_s, x, y_m, z) tmp = 0.0 if ((z <= 1.85e+49) || (!(z <= 1.5e+71) && (z <= 2.55e+124))) tmp = Float64(y_m * 0.5); else tmp = Float64(1.0 / Float64(Float64(1.0 / z) * Float64(Float64(1.0 / z) * Float64(y_m * -2.0)))); end return Float64(y_s * tmp) end
y\_m = abs(y); y\_s = sign(y) * abs(1.0); function tmp_2 = code(y_s, x, y_m, z) tmp = 0.0; if ((z <= 1.85e+49) || (~((z <= 1.5e+71)) && (z <= 2.55e+124))) tmp = y_m * 0.5; else tmp = 1.0 / ((1.0 / z) * ((1.0 / z) * (y_m * -2.0))); end tmp_2 = y_s * tmp; end
y\_m = N[Abs[y], $MachinePrecision]
y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
code[y$95$s_, x_, y$95$m_, z_] := N[(y$95$s * If[Or[LessEqual[z, 1.85e+49], And[N[Not[LessEqual[z, 1.5e+71]], $MachinePrecision], LessEqual[z, 2.55e+124]]], N[(y$95$m * 0.5), $MachinePrecision], N[(1.0 / N[(N[(1.0 / z), $MachinePrecision] * N[(N[(1.0 / z), $MachinePrecision] * N[(y$95$m * -2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
y\_m = \left|y\right|
\\
y\_s = \mathsf{copysign}\left(1, y\right)
\\
y\_s \cdot \begin{array}{l}
\mathbf{if}\;z \leq 1.85 \cdot 10^{+49} \lor \neg \left(z \leq 1.5 \cdot 10^{+71}\right) \land z \leq 2.55 \cdot 10^{+124}:\\
\;\;\;\;y\_m \cdot 0.5\\
\mathbf{else}:\\
\;\;\;\;\frac{1}{\frac{1}{z} \cdot \left(\frac{1}{z} \cdot \left(y\_m \cdot -2\right)\right)}\\
\end{array}
\end{array}
if z < 1.85000000000000009e49 or 1.50000000000000006e71 < z < 2.5499999999999999e124Initial program 67.1%
Taylor expanded in y around inf 44.3%
*-commutative44.3%
Simplified44.3%
if 1.85000000000000009e49 < z < 1.50000000000000006e71 or 2.5499999999999999e124 < z Initial program 62.4%
Taylor expanded in z around inf 70.0%
associate-*r/70.0%
metadata-eval70.0%
distribute-lft-neg-in70.0%
*-commutative70.0%
distribute-neg-frac70.0%
associate-*r/70.0%
distribute-rgt-neg-in70.0%
distribute-neg-frac70.0%
metadata-eval70.0%
Simplified70.0%
associate-*r/70.0%
clear-num70.0%
Applied egg-rr70.0%
add-cube-cbrt69.9%
pow369.9%
*-un-lft-identity69.9%
times-frac69.9%
pow-flip70.7%
metadata-eval70.7%
Applied egg-rr70.7%
rem-cube-cbrt70.8%
sqr-pow70.8%
associate-*l*75.5%
metadata-eval75.5%
unpow-175.5%
metadata-eval75.5%
unpow-175.5%
div-inv75.5%
metadata-eval75.5%
Applied egg-rr75.5%
Final simplification50.3%
y\_m = (fabs.f64 y)
y\_s = (copysign.f64 #s(literal 1 binary64) y)
(FPCore (y_s x y_m z)
:precision binary64
(*
y_s
(if (<= y_m 1.05e+160)
(/ (- (+ (* x x) (* y_m y_m)) (* z z)) (* y_m 2.0))
(* y_m 0.5))))y\_m = fabs(y);
y\_s = copysign(1.0, y);
double code(double y_s, double x, double y_m, double z) {
double tmp;
if (y_m <= 1.05e+160) {
tmp = (((x * x) + (y_m * y_m)) - (z * z)) / (y_m * 2.0);
} else {
tmp = y_m * 0.5;
}
return y_s * tmp;
}
y\_m = abs(y)
y\_s = copysign(1.0d0, y)
real(8) function code(y_s, x, y_m, z)
real(8), intent (in) :: y_s
real(8), intent (in) :: x
real(8), intent (in) :: y_m
real(8), intent (in) :: z
real(8) :: tmp
if (y_m <= 1.05d+160) then
tmp = (((x * x) + (y_m * y_m)) - (z * z)) / (y_m * 2.0d0)
else
tmp = y_m * 0.5d0
end if
code = y_s * tmp
end function
y\_m = Math.abs(y);
y\_s = Math.copySign(1.0, y);
public static double code(double y_s, double x, double y_m, double z) {
double tmp;
if (y_m <= 1.05e+160) {
tmp = (((x * x) + (y_m * y_m)) - (z * z)) / (y_m * 2.0);
} else {
tmp = y_m * 0.5;
}
return y_s * tmp;
}
y\_m = math.fabs(y) y\_s = math.copysign(1.0, y) def code(y_s, x, y_m, z): tmp = 0 if y_m <= 1.05e+160: tmp = (((x * x) + (y_m * y_m)) - (z * z)) / (y_m * 2.0) else: tmp = y_m * 0.5 return y_s * tmp
y\_m = abs(y) y\_s = copysign(1.0, y) function code(y_s, x, y_m, z) tmp = 0.0 if (y_m <= 1.05e+160) tmp = Float64(Float64(Float64(Float64(x * x) + Float64(y_m * y_m)) - Float64(z * z)) / Float64(y_m * 2.0)); else tmp = Float64(y_m * 0.5); end return Float64(y_s * tmp) end
y\_m = abs(y); y\_s = sign(y) * abs(1.0); function tmp_2 = code(y_s, x, y_m, z) tmp = 0.0; if (y_m <= 1.05e+160) tmp = (((x * x) + (y_m * y_m)) - (z * z)) / (y_m * 2.0); else tmp = y_m * 0.5; end tmp_2 = y_s * tmp; end
y\_m = N[Abs[y], $MachinePrecision]
y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
code[y$95$s_, x_, y$95$m_, z_] := N[(y$95$s * If[LessEqual[y$95$m, 1.05e+160], N[(N[(N[(N[(x * x), $MachinePrecision] + N[(y$95$m * y$95$m), $MachinePrecision]), $MachinePrecision] - N[(z * z), $MachinePrecision]), $MachinePrecision] / N[(y$95$m * 2.0), $MachinePrecision]), $MachinePrecision], N[(y$95$m * 0.5), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
y\_m = \left|y\right|
\\
y\_s = \mathsf{copysign}\left(1, y\right)
\\
y\_s \cdot \begin{array}{l}
\mathbf{if}\;y\_m \leq 1.05 \cdot 10^{+160}:\\
\;\;\;\;\frac{\left(x \cdot x + y\_m \cdot y\_m\right) - z \cdot z}{y\_m \cdot 2}\\
\mathbf{else}:\\
\;\;\;\;y\_m \cdot 0.5\\
\end{array}
\end{array}
if y < 1.04999999999999998e160Initial program 74.4%
if 1.04999999999999998e160 < y Initial program 8.8%
Taylor expanded in y around inf 85.1%
*-commutative85.1%
Simplified85.1%
Final simplification75.7%
y\_m = (fabs.f64 y) y\_s = (copysign.f64 #s(literal 1 binary64) y) (FPCore (y_s x y_m z) :precision binary64 (* y_s (* y_m 0.5)))
y\_m = fabs(y);
y\_s = copysign(1.0, y);
double code(double y_s, double x, double y_m, double z) {
return y_s * (y_m * 0.5);
}
y\_m = abs(y)
y\_s = copysign(1.0d0, y)
real(8) function code(y_s, x, y_m, z)
real(8), intent (in) :: y_s
real(8), intent (in) :: x
real(8), intent (in) :: y_m
real(8), intent (in) :: z
code = y_s * (y_m * 0.5d0)
end function
y\_m = Math.abs(y);
y\_s = Math.copySign(1.0, y);
public static double code(double y_s, double x, double y_m, double z) {
return y_s * (y_m * 0.5);
}
y\_m = math.fabs(y) y\_s = math.copysign(1.0, y) def code(y_s, x, y_m, z): return y_s * (y_m * 0.5)
y\_m = abs(y) y\_s = copysign(1.0, y) function code(y_s, x, y_m, z) return Float64(y_s * Float64(y_m * 0.5)) end
y\_m = abs(y); y\_s = sign(y) * abs(1.0); function tmp = code(y_s, x, y_m, z) tmp = y_s * (y_m * 0.5); end
y\_m = N[Abs[y], $MachinePrecision]
y\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[y]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
code[y$95$s_, x_, y$95$m_, z_] := N[(y$95$s * N[(y$95$m * 0.5), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
y\_m = \left|y\right|
\\
y\_s = \mathsf{copysign}\left(1, y\right)
\\
y\_s \cdot \left(y\_m \cdot 0.5\right)
\end{array}
Initial program 66.2%
Taylor expanded in y around inf 37.6%
*-commutative37.6%
Simplified37.6%
Final simplification37.6%
(FPCore (x y z) :precision binary64 (- (* y 0.5) (* (* (/ 0.5 y) (+ z x)) (- z x))))
double code(double x, double y, double z) {
return (y * 0.5) - (((0.5 / y) * (z + x)) * (z - x));
}
real(8) function code(x, y, z)
real(8), intent (in) :: x
real(8), intent (in) :: y
real(8), intent (in) :: z
code = (y * 0.5d0) - (((0.5d0 / y) * (z + x)) * (z - x))
end function
public static double code(double x, double y, double z) {
return (y * 0.5) - (((0.5 / y) * (z + x)) * (z - x));
}
def code(x, y, z): return (y * 0.5) - (((0.5 / y) * (z + x)) * (z - x))
function code(x, y, z) return Float64(Float64(y * 0.5) - Float64(Float64(Float64(0.5 / y) * Float64(z + x)) * Float64(z - x))) end
function tmp = code(x, y, z) tmp = (y * 0.5) - (((0.5 / y) * (z + x)) * (z - x)); end
code[x_, y_, z_] := N[(N[(y * 0.5), $MachinePrecision] - N[(N[(N[(0.5 / y), $MachinePrecision] * N[(z + x), $MachinePrecision]), $MachinePrecision] * N[(z - x), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
y \cdot 0.5 - \left(\frac{0.5}{y} \cdot \left(z + x\right)\right) \cdot \left(z - x\right)
\end{array}
herbie shell --seed 2024076
(FPCore (x y z)
:name "Diagrams.TwoD.Apollonian:initialConfig from diagrams-contrib-1.3.0.5, A"
:precision binary64
:alt
(- (* y 0.5) (* (* (/ 0.5 y) (+ z x)) (- z x)))
(/ (- (+ (* x x) (* y y)) (* z z)) (* y 2.0)))