
(FPCore (g h a) :precision binary64 (let* ((t_0 (/ 1.0 (* 2.0 a))) (t_1 (sqrt (- (* g g) (* h h))))) (+ (cbrt (* t_0 (+ (- g) t_1))) (cbrt (* t_0 (- (- g) t_1))))))
double code(double g, double h, double a) {
double t_0 = 1.0 / (2.0 * a);
double t_1 = sqrt(((g * g) - (h * h)));
return cbrt((t_0 * (-g + t_1))) + cbrt((t_0 * (-g - t_1)));
}
public static double code(double g, double h, double a) {
double t_0 = 1.0 / (2.0 * a);
double t_1 = Math.sqrt(((g * g) - (h * h)));
return Math.cbrt((t_0 * (-g + t_1))) + Math.cbrt((t_0 * (-g - t_1)));
}
function code(g, h, a) t_0 = Float64(1.0 / Float64(2.0 * a)) t_1 = sqrt(Float64(Float64(g * g) - Float64(h * h))) return Float64(cbrt(Float64(t_0 * Float64(Float64(-g) + t_1))) + cbrt(Float64(t_0 * Float64(Float64(-g) - t_1)))) end
code[g_, h_, a_] := Block[{t$95$0 = N[(1.0 / N[(2.0 * a), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[Sqrt[N[(N[(g * g), $MachinePrecision] - N[(h * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]}, N[(N[Power[N[(t$95$0 * N[((-g) + t$95$1), $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision] + N[Power[N[(t$95$0 * N[((-g) - t$95$1), $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision]), $MachinePrecision]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \frac{1}{2 \cdot a}\\
t_1 := \sqrt{g \cdot g - h \cdot h}\\
\sqrt[3]{t\_0 \cdot \left(\left(-g\right) + t\_1\right)} + \sqrt[3]{t\_0 \cdot \left(\left(-g\right) - t\_1\right)}
\end{array}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 5 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (g h a) :precision binary64 (let* ((t_0 (/ 1.0 (* 2.0 a))) (t_1 (sqrt (- (* g g) (* h h))))) (+ (cbrt (* t_0 (+ (- g) t_1))) (cbrt (* t_0 (- (- g) t_1))))))
double code(double g, double h, double a) {
double t_0 = 1.0 / (2.0 * a);
double t_1 = sqrt(((g * g) - (h * h)));
return cbrt((t_0 * (-g + t_1))) + cbrt((t_0 * (-g - t_1)));
}
public static double code(double g, double h, double a) {
double t_0 = 1.0 / (2.0 * a);
double t_1 = Math.sqrt(((g * g) - (h * h)));
return Math.cbrt((t_0 * (-g + t_1))) + Math.cbrt((t_0 * (-g - t_1)));
}
function code(g, h, a) t_0 = Float64(1.0 / Float64(2.0 * a)) t_1 = sqrt(Float64(Float64(g * g) - Float64(h * h))) return Float64(cbrt(Float64(t_0 * Float64(Float64(-g) + t_1))) + cbrt(Float64(t_0 * Float64(Float64(-g) - t_1)))) end
code[g_, h_, a_] := Block[{t$95$0 = N[(1.0 / N[(2.0 * a), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[Sqrt[N[(N[(g * g), $MachinePrecision] - N[(h * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]}, N[(N[Power[N[(t$95$0 * N[((-g) + t$95$1), $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision] + N[Power[N[(t$95$0 * N[((-g) - t$95$1), $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision]), $MachinePrecision]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \frac{1}{2 \cdot a}\\
t_1 := \sqrt{g \cdot g - h \cdot h}\\
\sqrt[3]{t\_0 \cdot \left(\left(-g\right) + t\_1\right)} + \sqrt[3]{t\_0 \cdot \left(\left(-g\right) - t\_1\right)}
\end{array}
\end{array}
h_m = (fabs.f64 h)
(FPCore (g h_m a)
:precision binary64
(let* ((t_0 (pow (* 0.0 g) 2.0)))
(if (<= h_m 1.48e-170)
(fma
(* (/ (cbrt g) (cbrt a)) (cbrt -0.5))
(pow 2.0 0.3333333333333333)
(* (cbrt (* (/ h_m g) (/ h_m a))) (* (cbrt 0.5) (cbrt -0.5))))
(/
(+
(cbrt
(*
h_m
(fma
-0.5
(* h_m (/ (fma 0.25 (/ t_0 (* g g)) 1.0) g))
(* 0.5 (/ (* 0.0 g) g)))))
(cbrt
(-
(- g)
(fma
h_m
(/
(fma -0.5 (* h_m (fma (/ 0.25 g) (/ t_0 g) 1.0)) (* 0.5 (* 0.0 g)))
g)
g))))
(cbrt (* a 2.0))))))h_m = fabs(h);
double code(double g, double h_m, double a) {
double t_0 = pow((0.0 * g), 2.0);
double tmp;
if (h_m <= 1.48e-170) {
tmp = fma(((cbrt(g) / cbrt(a)) * cbrt(-0.5)), pow(2.0, 0.3333333333333333), (cbrt(((h_m / g) * (h_m / a))) * (cbrt(0.5) * cbrt(-0.5))));
} else {
tmp = (cbrt((h_m * fma(-0.5, (h_m * (fma(0.25, (t_0 / (g * g)), 1.0) / g)), (0.5 * ((0.0 * g) / g))))) + cbrt((-g - fma(h_m, (fma(-0.5, (h_m * fma((0.25 / g), (t_0 / g), 1.0)), (0.5 * (0.0 * g))) / g), g)))) / cbrt((a * 2.0));
}
return tmp;
}
h_m = abs(h) function code(g, h_m, a) t_0 = Float64(0.0 * g) ^ 2.0 tmp = 0.0 if (h_m <= 1.48e-170) tmp = fma(Float64(Float64(cbrt(g) / cbrt(a)) * cbrt(-0.5)), (2.0 ^ 0.3333333333333333), Float64(cbrt(Float64(Float64(h_m / g) * Float64(h_m / a))) * Float64(cbrt(0.5) * cbrt(-0.5)))); else tmp = Float64(Float64(cbrt(Float64(h_m * fma(-0.5, Float64(h_m * Float64(fma(0.25, Float64(t_0 / Float64(g * g)), 1.0) / g)), Float64(0.5 * Float64(Float64(0.0 * g) / g))))) + cbrt(Float64(Float64(-g) - fma(h_m, Float64(fma(-0.5, Float64(h_m * fma(Float64(0.25 / g), Float64(t_0 / g), 1.0)), Float64(0.5 * Float64(0.0 * g))) / g), g)))) / cbrt(Float64(a * 2.0))); end return tmp end
h_m = N[Abs[h], $MachinePrecision]
code[g_, h$95$m_, a_] := Block[{t$95$0 = N[Power[N[(0.0 * g), $MachinePrecision], 2.0], $MachinePrecision]}, If[LessEqual[h$95$m, 1.48e-170], N[(N[(N[(N[Power[g, 1/3], $MachinePrecision] / N[Power[a, 1/3], $MachinePrecision]), $MachinePrecision] * N[Power[-0.5, 1/3], $MachinePrecision]), $MachinePrecision] * N[Power[2.0, 0.3333333333333333], $MachinePrecision] + N[(N[Power[N[(N[(h$95$m / g), $MachinePrecision] * N[(h$95$m / a), $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision] * N[(N[Power[0.5, 1/3], $MachinePrecision] * N[Power[-0.5, 1/3], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(N[Power[N[(h$95$m * N[(-0.5 * N[(h$95$m * N[(N[(0.25 * N[(t$95$0 / N[(g * g), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision] / g), $MachinePrecision]), $MachinePrecision] + N[(0.5 * N[(N[(0.0 * g), $MachinePrecision] / g), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision] + N[Power[N[((-g) - N[(h$95$m * N[(N[(-0.5 * N[(h$95$m * N[(N[(0.25 / g), $MachinePrecision] * N[(t$95$0 / g), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision] + N[(0.5 * N[(0.0 * g), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / g), $MachinePrecision] + g), $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision]), $MachinePrecision] / N[Power[N[(a * 2.0), $MachinePrecision], 1/3], $MachinePrecision]), $MachinePrecision]]]
\begin{array}{l}
h_m = \left|h\right|
\\
\begin{array}{l}
t_0 := {\left(0 \cdot g\right)}^{2}\\
\mathbf{if}\;h\_m \leq 1.48 \cdot 10^{-170}:\\
\;\;\;\;\mathsf{fma}\left(\frac{\sqrt[3]{g}}{\sqrt[3]{a}} \cdot \sqrt[3]{-0.5}, {2}^{0.3333333333333333}, \sqrt[3]{\frac{h\_m}{g} \cdot \frac{h\_m}{a}} \cdot \left(\sqrt[3]{0.5} \cdot \sqrt[3]{-0.5}\right)\right)\\
\mathbf{else}:\\
\;\;\;\;\frac{\sqrt[3]{h\_m \cdot \mathsf{fma}\left(-0.5, h\_m \cdot \frac{\mathsf{fma}\left(0.25, \frac{t\_0}{g \cdot g}, 1\right)}{g}, 0.5 \cdot \frac{0 \cdot g}{g}\right)} + \sqrt[3]{\left(-g\right) - \mathsf{fma}\left(h\_m, \frac{\mathsf{fma}\left(-0.5, h\_m \cdot \mathsf{fma}\left(\frac{0.25}{g}, \frac{t\_0}{g}, 1\right), 0.5 \cdot \left(0 \cdot g\right)\right)}{g}, g\right)}}{\sqrt[3]{a \cdot 2}}\\
\end{array}
\end{array}
if h < 1.48000000000000005e-170Initial program 48.4%
Taylor expanded in h around 0
associate-*r*N/A
lower-fma.f64N/A
lower-*.f64N/A
lower-cbrt.f64N/A
lower-/.f64N/A
lower-cbrt.f64N/A
lower-cbrt.f64N/A
lower-*.f64N/A
lower-cbrt.f64N/A
unpow2N/A
*-commutativeN/A
times-fracN/A
lower-*.f64N/A
lower-/.f64N/A
lower-/.f64N/A
*-commutativeN/A
lower-*.f64N/A
lower-cbrt.f64N/A
lower-cbrt.f6475.4
Applied rewrites75.4%
Applied rewrites93.5%
Applied rewrites94.2%
if 1.48000000000000005e-170 < h Initial program 42.6%
Applied rewrites9.7%
Taylor expanded in h around 0
lower-*.f64N/A
lower-fma.f64N/A
Applied rewrites20.9%
Taylor expanded in h around 0
+-commutativeN/A
lower-fma.f64N/A
Applied rewrites98.8%
h_m = (fabs.f64 h)
(FPCore (g h_m a)
:precision binary64
(let* ((t_0 (pow (* 0.0 g) 2.0)))
(if (<= h_m 5e-169)
(/ (cbrt (* g -1.0)) (cbrt a))
(/
(+
(cbrt
(*
h_m
(fma
-0.5
(* h_m (/ (fma 0.25 (/ t_0 (* g g)) 1.0) g))
(* 0.5 (/ (* 0.0 g) g)))))
(cbrt
(-
(- g)
(fma
h_m
(/
(fma -0.5 (* h_m (fma (/ 0.25 g) (/ t_0 g) 1.0)) (* 0.5 (* 0.0 g)))
g)
g))))
(cbrt (* a 2.0))))))h_m = fabs(h);
double code(double g, double h_m, double a) {
double t_0 = pow((0.0 * g), 2.0);
double tmp;
if (h_m <= 5e-169) {
tmp = cbrt((g * -1.0)) / cbrt(a);
} else {
tmp = (cbrt((h_m * fma(-0.5, (h_m * (fma(0.25, (t_0 / (g * g)), 1.0) / g)), (0.5 * ((0.0 * g) / g))))) + cbrt((-g - fma(h_m, (fma(-0.5, (h_m * fma((0.25 / g), (t_0 / g), 1.0)), (0.5 * (0.0 * g))) / g), g)))) / cbrt((a * 2.0));
}
return tmp;
}
h_m = abs(h) function code(g, h_m, a) t_0 = Float64(0.0 * g) ^ 2.0 tmp = 0.0 if (h_m <= 5e-169) tmp = Float64(cbrt(Float64(g * -1.0)) / cbrt(a)); else tmp = Float64(Float64(cbrt(Float64(h_m * fma(-0.5, Float64(h_m * Float64(fma(0.25, Float64(t_0 / Float64(g * g)), 1.0) / g)), Float64(0.5 * Float64(Float64(0.0 * g) / g))))) + cbrt(Float64(Float64(-g) - fma(h_m, Float64(fma(-0.5, Float64(h_m * fma(Float64(0.25 / g), Float64(t_0 / g), 1.0)), Float64(0.5 * Float64(0.0 * g))) / g), g)))) / cbrt(Float64(a * 2.0))); end return tmp end
h_m = N[Abs[h], $MachinePrecision]
code[g_, h$95$m_, a_] := Block[{t$95$0 = N[Power[N[(0.0 * g), $MachinePrecision], 2.0], $MachinePrecision]}, If[LessEqual[h$95$m, 5e-169], N[(N[Power[N[(g * -1.0), $MachinePrecision], 1/3], $MachinePrecision] / N[Power[a, 1/3], $MachinePrecision]), $MachinePrecision], N[(N[(N[Power[N[(h$95$m * N[(-0.5 * N[(h$95$m * N[(N[(0.25 * N[(t$95$0 / N[(g * g), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision] / g), $MachinePrecision]), $MachinePrecision] + N[(0.5 * N[(N[(0.0 * g), $MachinePrecision] / g), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision] + N[Power[N[((-g) - N[(h$95$m * N[(N[(-0.5 * N[(h$95$m * N[(N[(0.25 / g), $MachinePrecision] * N[(t$95$0 / g), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision] + N[(0.5 * N[(0.0 * g), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / g), $MachinePrecision] + g), $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision]), $MachinePrecision] / N[Power[N[(a * 2.0), $MachinePrecision], 1/3], $MachinePrecision]), $MachinePrecision]]]
\begin{array}{l}
h_m = \left|h\right|
\\
\begin{array}{l}
t_0 := {\left(0 \cdot g\right)}^{2}\\
\mathbf{if}\;h\_m \leq 5 \cdot 10^{-169}:\\
\;\;\;\;\frac{\sqrt[3]{g \cdot -1}}{\sqrt[3]{a}}\\
\mathbf{else}:\\
\;\;\;\;\frac{\sqrt[3]{h\_m \cdot \mathsf{fma}\left(-0.5, h\_m \cdot \frac{\mathsf{fma}\left(0.25, \frac{t\_0}{g \cdot g}, 1\right)}{g}, 0.5 \cdot \frac{0 \cdot g}{g}\right)} + \sqrt[3]{\left(-g\right) - \mathsf{fma}\left(h\_m, \frac{\mathsf{fma}\left(-0.5, h\_m \cdot \mathsf{fma}\left(\frac{0.25}{g}, \frac{t\_0}{g}, 1\right), 0.5 \cdot \left(0 \cdot g\right)\right)}{g}, g\right)}}{\sqrt[3]{a \cdot 2}}\\
\end{array}
\end{array}
if h < 5.0000000000000002e-169Initial program 48.4%
Applied rewrites10.8%
Taylor expanded in g around inf
lower-*.f64N/A
lower-cbrt.f64N/A
lower-/.f64N/A
lower-cbrt.f6474.7
Applied rewrites74.7%
Applied rewrites96.0%
if 5.0000000000000002e-169 < h Initial program 42.6%
Applied rewrites9.7%
Taylor expanded in h around 0
lower-*.f64N/A
lower-fma.f64N/A
Applied rewrites20.9%
Taylor expanded in h around 0
+-commutativeN/A
lower-fma.f64N/A
Applied rewrites98.8%
h_m = (fabs.f64 h)
(FPCore (g h_m a)
:precision binary64
(if (<= h_m 5e-169)
(/ (cbrt (* g -1.0)) (cbrt a))
(/
(+
(cbrt
(*
h_m
(fma
-0.5
(* h_m (/ (fma 0.25 (/ (pow (* 0.0 g) 2.0) (* g g)) 1.0) g))
(* 0.5 (/ (* 0.0 g) g)))))
(cbrt (- (- g) (* (- g) -1.0))))
(cbrt (* a 2.0)))))h_m = fabs(h);
double code(double g, double h_m, double a) {
double tmp;
if (h_m <= 5e-169) {
tmp = cbrt((g * -1.0)) / cbrt(a);
} else {
tmp = (cbrt((h_m * fma(-0.5, (h_m * (fma(0.25, (pow((0.0 * g), 2.0) / (g * g)), 1.0) / g)), (0.5 * ((0.0 * g) / g))))) + cbrt((-g - (-g * -1.0)))) / cbrt((a * 2.0));
}
return tmp;
}
h_m = abs(h) function code(g, h_m, a) tmp = 0.0 if (h_m <= 5e-169) tmp = Float64(cbrt(Float64(g * -1.0)) / cbrt(a)); else tmp = Float64(Float64(cbrt(Float64(h_m * fma(-0.5, Float64(h_m * Float64(fma(0.25, Float64((Float64(0.0 * g) ^ 2.0) / Float64(g * g)), 1.0) / g)), Float64(0.5 * Float64(Float64(0.0 * g) / g))))) + cbrt(Float64(Float64(-g) - Float64(Float64(-g) * -1.0)))) / cbrt(Float64(a * 2.0))); end return tmp end
h_m = N[Abs[h], $MachinePrecision] code[g_, h$95$m_, a_] := If[LessEqual[h$95$m, 5e-169], N[(N[Power[N[(g * -1.0), $MachinePrecision], 1/3], $MachinePrecision] / N[Power[a, 1/3], $MachinePrecision]), $MachinePrecision], N[(N[(N[Power[N[(h$95$m * N[(-0.5 * N[(h$95$m * N[(N[(0.25 * N[(N[Power[N[(0.0 * g), $MachinePrecision], 2.0], $MachinePrecision] / N[(g * g), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision] / g), $MachinePrecision]), $MachinePrecision] + N[(0.5 * N[(N[(0.0 * g), $MachinePrecision] / g), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision] + N[Power[N[((-g) - N[((-g) * -1.0), $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision]), $MachinePrecision] / N[Power[N[(a * 2.0), $MachinePrecision], 1/3], $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
h_m = \left|h\right|
\\
\begin{array}{l}
\mathbf{if}\;h\_m \leq 5 \cdot 10^{-169}:\\
\;\;\;\;\frac{\sqrt[3]{g \cdot -1}}{\sqrt[3]{a}}\\
\mathbf{else}:\\
\;\;\;\;\frac{\sqrt[3]{h\_m \cdot \mathsf{fma}\left(-0.5, h\_m \cdot \frac{\mathsf{fma}\left(0.25, \frac{{\left(0 \cdot g\right)}^{2}}{g \cdot g}, 1\right)}{g}, 0.5 \cdot \frac{0 \cdot g}{g}\right)} + \sqrt[3]{\left(-g\right) - \left(-g\right) \cdot -1}}{\sqrt[3]{a \cdot 2}}\\
\end{array}
\end{array}
if h < 5.0000000000000002e-169Initial program 48.4%
Applied rewrites10.8%
Taylor expanded in g around inf
lower-*.f64N/A
lower-cbrt.f64N/A
lower-/.f64N/A
lower-cbrt.f6474.7
Applied rewrites74.7%
Applied rewrites96.0%
if 5.0000000000000002e-169 < h Initial program 42.6%
Applied rewrites9.7%
Taylor expanded in h around 0
lower-*.f64N/A
lower-fma.f64N/A
Applied rewrites20.9%
Taylor expanded in g around -inf
associate-*r*N/A
unpow2N/A
rem-square-sqrtN/A
lower-*.f64N/A
mul-1-negN/A
lower-neg.f6498.8
Applied rewrites98.8%
h_m = (fabs.f64 h) (FPCore (g h_m a) :precision binary64 (/ (cbrt (* g -1.0)) (cbrt a)))
h_m = fabs(h);
double code(double g, double h_m, double a) {
return cbrt((g * -1.0)) / cbrt(a);
}
h_m = Math.abs(h);
public static double code(double g, double h_m, double a) {
return Math.cbrt((g * -1.0)) / Math.cbrt(a);
}
h_m = abs(h) function code(g, h_m, a) return Float64(cbrt(Float64(g * -1.0)) / cbrt(a)) end
h_m = N[Abs[h], $MachinePrecision] code[g_, h$95$m_, a_] := N[(N[Power[N[(g * -1.0), $MachinePrecision], 1/3], $MachinePrecision] / N[Power[a, 1/3], $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
h_m = \left|h\right|
\\
\frac{\sqrt[3]{g \cdot -1}}{\sqrt[3]{a}}
\end{array}
Initial program 46.8%
Applied rewrites10.5%
Taylor expanded in g around inf
lower-*.f64N/A
lower-cbrt.f64N/A
lower-/.f64N/A
lower-cbrt.f6472.7
Applied rewrites72.7%
Applied rewrites96.0%
h_m = (fabs.f64 h) (FPCore (g h_m a) :precision binary64 (cbrt (* (/ g a) -1.0)))
h_m = fabs(h);
double code(double g, double h_m, double a) {
return cbrt(((g / a) * -1.0));
}
h_m = Math.abs(h);
public static double code(double g, double h_m, double a) {
return Math.cbrt(((g / a) * -1.0));
}
h_m = abs(h) function code(g, h_m, a) return cbrt(Float64(Float64(g / a) * -1.0)) end
h_m = N[Abs[h], $MachinePrecision] code[g_, h$95$m_, a_] := N[Power[N[(N[(g / a), $MachinePrecision] * -1.0), $MachinePrecision], 1/3], $MachinePrecision]
\begin{array}{l}
h_m = \left|h\right|
\\
\sqrt[3]{\frac{g}{a} \cdot -1}
\end{array}
Initial program 46.8%
Applied rewrites10.5%
Taylor expanded in g around inf
lower-*.f64N/A
lower-cbrt.f64N/A
lower-/.f64N/A
lower-cbrt.f6472.7
Applied rewrites72.7%
Applied rewrites72.7%
herbie shell --seed 2025003
(FPCore (g h a)
:name "2-ancestry mixing, positive discriminant"
:precision binary64
(+ (cbrt (* (/ 1.0 (* 2.0 a)) (+ (- g) (sqrt (- (* g g) (* h h)))))) (cbrt (* (/ 1.0 (* 2.0 a)) (- (- g) (sqrt (- (* g g) (* h h))))))))