\sqrt[3]{\frac{1}{2 \cdot a} \cdot \left(\left(-g\right) + \sqrt{g \cdot g - h \cdot h}\right)} + \sqrt[3]{\frac{1}{2 \cdot a} \cdot \left(\left(-g\right) - \sqrt{g \cdot g - h \cdot h}\right)}
\begin{array}{l}
t_0 := \sqrt{{g}^{2} - {h}^{2}}\\
t_1 := \log \left(g + t_0\right)\\
t_2 := \sqrt{g \cdot g - h \cdot h}\\
t_3 := \sqrt[3]{\frac{1}{2 \cdot a} \cdot \left(t_2 - g\right)}\\
\mathbf{if}\;2 \cdot a \leq -3.991902884307353 \cdot 10^{-306}:\\
\;\;\;\;t_3 + \sqrt[3]{0.5} \cdot e^{0.3333333333333333 \cdot \left(t_1 + \log \left(\frac{-1}{a}\right)\right)}\\
\mathbf{elif}\;2 \cdot a \leq 2.8974824116342774 \cdot 10^{-38}:\\
\;\;\;\;t_3 + \sqrt[3]{-0.5} \cdot e^{0.3333333333333333 \cdot \left(t_1 - \log a\right)}\\
\mathbf{else}:\\
\;\;\;\;\sqrt[3]{0.5} \cdot e^{0.3333333333333333 \cdot \left(\log \left(t_0 - g\right) - \log a\right)} + \sqrt[3]{\left(g + t_2\right) \cdot \frac{-1}{2 \cdot a}}\\
\end{array}
(FPCore (g h a) :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))))))))
(FPCore (g h a)
:precision binary64
(let* ((t_0 (sqrt (- (pow g 2.0) (pow h 2.0))))
(t_1 (log (+ g t_0)))
(t_2 (sqrt (- (* g g) (* h h))))
(t_3 (cbrt (* (/ 1.0 (* 2.0 a)) (- t_2 g)))))
(if (<= (* 2.0 a) -3.991902884307353e-306)
(+
t_3
(* (cbrt 0.5) (exp (* 0.3333333333333333 (+ t_1 (log (/ -1.0 a)))))))
(if (<= (* 2.0 a) 2.8974824116342774e-38)
(+ t_3 (* (cbrt -0.5) (exp (* 0.3333333333333333 (- t_1 (log a))))))
(+
(* (cbrt 0.5) (exp (* 0.3333333333333333 (- (log (- t_0 g)) (log a)))))
(cbrt (* (+ g t_2) (/ -1.0 (* 2.0 a)))))))))double code(double g, double h, double a) {
return cbrt(((1.0 / (2.0 * a)) * (-g + sqrt(((g * g) - (h * h)))))) + cbrt(((1.0 / (2.0 * a)) * (-g - sqrt(((g * g) - (h * h))))));
}
double code(double g, double h, double a) {
double t_0 = sqrt((pow(g, 2.0) - pow(h, 2.0)));
double t_1 = log((g + t_0));
double t_2 = sqrt(((g * g) - (h * h)));
double t_3 = cbrt(((1.0 / (2.0 * a)) * (t_2 - g)));
double tmp;
if ((2.0 * a) <= -3.991902884307353e-306) {
tmp = t_3 + (cbrt(0.5) * exp((0.3333333333333333 * (t_1 + log((-1.0 / a))))));
} else if ((2.0 * a) <= 2.8974824116342774e-38) {
tmp = t_3 + (cbrt(-0.5) * exp((0.3333333333333333 * (t_1 - log(a)))));
} else {
tmp = (cbrt(0.5) * exp((0.3333333333333333 * (log((t_0 - g)) - log(a))))) + cbrt(((g + t_2) * (-1.0 / (2.0 * a))));
}
return tmp;
}



Bits error versus g



Bits error versus h



Bits error versus a
Results
if (*.f64 2 a) < -3.99190288430735316e-306Initial program 35.6
Taylor expanded in a around -inf 35.4
if -3.99190288430735316e-306 < (*.f64 2 a) < 2.8974824116342774e-38Initial program 37.4
Taylor expanded in a around 0 36.9
if 2.8974824116342774e-38 < (*.f64 2 a) Initial program 34.5
Taylor expanded in a around 0 35.0
Final simplification35.6
herbie shell --seed 2022129
(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))))))))