R \cdot \left(2 \cdot \tan^{-1}_* \frac{\sqrt{{\left(\sin \left(\frac{\phi_1 - \phi_2}{2}\right)\right)}^{2} + \left(\left(\cos \phi_1 \cdot \cos \phi_2\right) \cdot \sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right) \cdot \sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)}}{\sqrt{1 - \left({\left(\sin \left(\frac{\phi_1 - \phi_2}{2}\right)\right)}^{2} + \left(\left(\cos \phi_1 \cdot \cos \phi_2\right) \cdot \sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right) \cdot \sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)}}\right)R \cdot \left(2 \cdot \tan^{-1}_* \frac{\sqrt{{\left(\sin \left(\frac{\phi_1}{2}\right) \cdot \cos \left(\frac{\phi_2}{2}\right) - \cos \left(\frac{\phi_1}{2}\right) \cdot \sin \left(\frac{\phi_2}{2}\right)\right)}^{2} + \left(\left(\cos \phi_1 \cdot \cos \phi_2\right) \cdot \sqrt[3]{{\left(\sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)}^{3}}\right) \cdot \mathsf{log1p}\left(\mathsf{expm1}\left(\sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)\right)}}{\sqrt{1 - \left({\left(\sin \left(\frac{\phi_1}{2}\right) \cdot \cos \left(\frac{\phi_2}{2}\right) - \cos \left(\frac{\phi_1}{2}\right) \cdot \sin \left(\frac{\phi_2}{2}\right)\right)}^{2} + \left(\left(\cos \phi_1 \cdot \cos \phi_2\right) \cdot \sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right) \cdot \sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)}}\right)double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r113596 = R;
double r113597 = 2.0;
double r113598 = phi1;
double r113599 = phi2;
double r113600 = r113598 - r113599;
double r113601 = r113600 / r113597;
double r113602 = sin(r113601);
double r113603 = pow(r113602, r113597);
double r113604 = cos(r113598);
double r113605 = cos(r113599);
double r113606 = r113604 * r113605;
double r113607 = lambda1;
double r113608 = lambda2;
double r113609 = r113607 - r113608;
double r113610 = r113609 / r113597;
double r113611 = sin(r113610);
double r113612 = r113606 * r113611;
double r113613 = r113612 * r113611;
double r113614 = r113603 + r113613;
double r113615 = sqrt(r113614);
double r113616 = 1.0;
double r113617 = r113616 - r113614;
double r113618 = sqrt(r113617);
double r113619 = atan2(r113615, r113618);
double r113620 = r113597 * r113619;
double r113621 = r113596 * r113620;
return r113621;
}
double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r113622 = R;
double r113623 = 2.0;
double r113624 = phi1;
double r113625 = r113624 / r113623;
double r113626 = sin(r113625);
double r113627 = phi2;
double r113628 = r113627 / r113623;
double r113629 = cos(r113628);
double r113630 = r113626 * r113629;
double r113631 = cos(r113625);
double r113632 = sin(r113628);
double r113633 = r113631 * r113632;
double r113634 = r113630 - r113633;
double r113635 = pow(r113634, r113623);
double r113636 = cos(r113624);
double r113637 = cos(r113627);
double r113638 = r113636 * r113637;
double r113639 = lambda1;
double r113640 = lambda2;
double r113641 = r113639 - r113640;
double r113642 = r113641 / r113623;
double r113643 = sin(r113642);
double r113644 = 3.0;
double r113645 = pow(r113643, r113644);
double r113646 = cbrt(r113645);
double r113647 = r113638 * r113646;
double r113648 = expm1(r113643);
double r113649 = log1p(r113648);
double r113650 = r113647 * r113649;
double r113651 = r113635 + r113650;
double r113652 = sqrt(r113651);
double r113653 = 1.0;
double r113654 = r113638 * r113643;
double r113655 = r113654 * r113643;
double r113656 = r113635 + r113655;
double r113657 = r113653 - r113656;
double r113658 = sqrt(r113657);
double r113659 = atan2(r113652, r113658);
double r113660 = r113623 * r113659;
double r113661 = r113622 * r113660;
return r113661;
}



Bits error versus R



Bits error versus lambda1



Bits error versus lambda2



Bits error versus phi1



Bits error versus phi2
Results
Initial program 23.8
rmApplied div-sub23.8
Applied sin-diff23.3
rmApplied div-sub23.3
Applied sin-diff13.9
rmApplied log1p-expm1-u13.9
rmApplied add-cbrt-cube14.1
Simplified14.1
Final simplification14.1
herbie shell --seed 2020003 +o rules:numerics
(FPCore (R lambda1 lambda2 phi1 phi2)
:name "Distance on a great circle"
:precision binary64
(* R (* 2 (atan2 (sqrt (+ (pow (sin (/ (- phi1 phi2) 2)) 2) (* (* (* (cos phi1) (cos phi2)) (sin (/ (- lambda1 lambda2) 2))) (sin (/ (- lambda1 lambda2) 2))))) (sqrt (- 1 (+ (pow (sin (/ (- phi1 phi2) 2)) 2) (* (* (* (cos phi1) (cos phi2)) (sin (/ (- lambda1 lambda2) 2))) (sin (/ (- lambda1 lambda2) 2))))))))))