\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}
\mathbf{if}\;g \le -2.5831109636145956 \cdot 10^{-186}:\\
\;\;\;\;\frac{\sqrt[3]{1 \cdot \left(\left(-g\right) + \sqrt{g \cdot g - h \cdot h}\right)}}{\sqrt[3]{2 \cdot a}} + \sqrt[3]{\frac{1}{2 \cdot a}} \cdot \sqrt[3]{\left(-g\right) - -1 \cdot g}\\
\mathbf{else}:\\
\;\;\;\;\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 \sqrt[3]{\left(-g\right) - g}\\
\end{array}double f(double g, double h, double a) {
double r154908 = 1.0;
double r154909 = 2.0;
double r154910 = a;
double r154911 = r154909 * r154910;
double r154912 = r154908 / r154911;
double r154913 = g;
double r154914 = -r154913;
double r154915 = r154913 * r154913;
double r154916 = h;
double r154917 = r154916 * r154916;
double r154918 = r154915 - r154917;
double r154919 = sqrt(r154918);
double r154920 = r154914 + r154919;
double r154921 = r154912 * r154920;
double r154922 = cbrt(r154921);
double r154923 = r154914 - r154919;
double r154924 = r154912 * r154923;
double r154925 = cbrt(r154924);
double r154926 = r154922 + r154925;
return r154926;
}
double f(double g, double h, double a) {
double r154927 = g;
double r154928 = -2.5831109636145956e-186;
bool r154929 = r154927 <= r154928;
double r154930 = 1.0;
double r154931 = -r154927;
double r154932 = r154927 * r154927;
double r154933 = h;
double r154934 = r154933 * r154933;
double r154935 = r154932 - r154934;
double r154936 = sqrt(r154935);
double r154937 = r154931 + r154936;
double r154938 = r154930 * r154937;
double r154939 = cbrt(r154938);
double r154940 = 2.0;
double r154941 = a;
double r154942 = r154940 * r154941;
double r154943 = cbrt(r154942);
double r154944 = r154939 / r154943;
double r154945 = r154930 / r154942;
double r154946 = cbrt(r154945);
double r154947 = -1.0;
double r154948 = r154947 * r154927;
double r154949 = r154931 - r154948;
double r154950 = cbrt(r154949);
double r154951 = r154946 * r154950;
double r154952 = r154944 + r154951;
double r154953 = r154945 * r154937;
double r154954 = cbrt(r154953);
double r154955 = r154931 - r154927;
double r154956 = cbrt(r154955);
double r154957 = r154946 * r154956;
double r154958 = r154954 + r154957;
double r154959 = r154929 ? r154952 : r154958;
return r154959;
}



Bits error versus g



Bits error versus h



Bits error versus a
Results
if g < -2.5831109636145956e-186Initial program 34.9
rmApplied cbrt-prod34.8
rmApplied associate-*l/34.8
Applied cbrt-div31.4
Taylor expanded around -inf 31.7
if -2.5831109636145956e-186 < g Initial program 36.4
rmApplied cbrt-prod33.0
Taylor expanded around inf 31.9
Final simplification31.8
herbie shell --seed 2020046 +o rules:numerics
(FPCore (g h a)
:name "2-ancestry mixing, positive discriminant"
:precision binary64
(+ (cbrt (* (/ 1 (* 2 a)) (+ (- g) (sqrt (- (* g g) (* h h)))))) (cbrt (* (/ 1 (* 2 a)) (- (- g) (sqrt (- (* g g) (* h h))))))))