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 - \phi_2}{2}\right)\right)}^{2} + \left(\left(\cos \phi_1 \cdot \cos \phi_2\right) \cdot \mathsf{expm1}\left(\mathsf{log1p}\left(\sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)\right)\right) \cdot \mathsf{expm1}\left(\mathsf{log1p}\left(\sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)\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 \mathsf{expm1}\left(\mathsf{log1p}\left(\sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)\right)\right)}}\right)double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r81584 = R;
double r81585 = 2.0;
double r81586 = phi1;
double r81587 = phi2;
double r81588 = r81586 - r81587;
double r81589 = r81588 / r81585;
double r81590 = sin(r81589);
double r81591 = pow(r81590, r81585);
double r81592 = cos(r81586);
double r81593 = cos(r81587);
double r81594 = r81592 * r81593;
double r81595 = lambda1;
double r81596 = lambda2;
double r81597 = r81595 - r81596;
double r81598 = r81597 / r81585;
double r81599 = sin(r81598);
double r81600 = r81594 * r81599;
double r81601 = r81600 * r81599;
double r81602 = r81591 + r81601;
double r81603 = sqrt(r81602);
double r81604 = 1.0;
double r81605 = r81604 - r81602;
double r81606 = sqrt(r81605);
double r81607 = atan2(r81603, r81606);
double r81608 = r81585 * r81607;
double r81609 = r81584 * r81608;
return r81609;
}
double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r81610 = R;
double r81611 = 2.0;
double r81612 = phi1;
double r81613 = phi2;
double r81614 = r81612 - r81613;
double r81615 = r81614 / r81611;
double r81616 = sin(r81615);
double r81617 = pow(r81616, r81611);
double r81618 = cos(r81612);
double r81619 = cos(r81613);
double r81620 = r81618 * r81619;
double r81621 = lambda1;
double r81622 = lambda2;
double r81623 = r81621 - r81622;
double r81624 = r81623 / r81611;
double r81625 = sin(r81624);
double r81626 = log1p(r81625);
double r81627 = expm1(r81626);
double r81628 = r81620 * r81627;
double r81629 = r81628 * r81627;
double r81630 = r81617 + r81629;
double r81631 = sqrt(r81630);
double r81632 = 1.0;
double r81633 = r81620 * r81625;
double r81634 = r81633 * r81627;
double r81635 = r81617 + r81634;
double r81636 = r81632 - r81635;
double r81637 = sqrt(r81636);
double r81638 = atan2(r81631, r81637);
double r81639 = r81611 * r81638;
double r81640 = r81610 * r81639;
return r81640;
}



Bits error versus R



Bits error versus lambda1



Bits error versus lambda2



Bits error versus phi1



Bits error versus phi2
Results
Initial program 24.7
rmApplied expm1-log1p-u24.7
rmApplied expm1-log1p-u24.7
rmApplied expm1-log1p-u24.7
Final simplification24.7
herbie shell --seed 2020018 +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))))))))))