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 \sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right) \cdot \sqrt[3]{{\left(\mathsf{log1p}\left(\mathsf{expm1}\left(\sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)\right)\right)}^{3}}}}{\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 \mathsf{expm1}\left(\mathsf{log1p}\left(\sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)\right)\right) \cdot \mathsf{log1p}\left(\mathsf{expm1}\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 r111659 = R;
double r111660 = 2.0;
double r111661 = phi1;
double r111662 = phi2;
double r111663 = r111661 - r111662;
double r111664 = r111663 / r111660;
double r111665 = sin(r111664);
double r111666 = pow(r111665, r111660);
double r111667 = cos(r111661);
double r111668 = cos(r111662);
double r111669 = r111667 * r111668;
double r111670 = lambda1;
double r111671 = lambda2;
double r111672 = r111670 - r111671;
double r111673 = r111672 / r111660;
double r111674 = sin(r111673);
double r111675 = r111669 * r111674;
double r111676 = r111675 * r111674;
double r111677 = r111666 + r111676;
double r111678 = sqrt(r111677);
double r111679 = 1.0;
double r111680 = r111679 - r111677;
double r111681 = sqrt(r111680);
double r111682 = atan2(r111678, r111681);
double r111683 = r111660 * r111682;
double r111684 = r111659 * r111683;
return r111684;
}
double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r111685 = R;
double r111686 = 2.0;
double r111687 = phi1;
double r111688 = phi2;
double r111689 = r111687 - r111688;
double r111690 = r111689 / r111686;
double r111691 = sin(r111690);
double r111692 = pow(r111691, r111686);
double r111693 = cos(r111687);
double r111694 = cos(r111688);
double r111695 = r111693 * r111694;
double r111696 = lambda1;
double r111697 = lambda2;
double r111698 = r111696 - r111697;
double r111699 = r111698 / r111686;
double r111700 = sin(r111699);
double r111701 = r111695 * r111700;
double r111702 = expm1(r111700);
double r111703 = log1p(r111702);
double r111704 = 3.0;
double r111705 = pow(r111703, r111704);
double r111706 = cbrt(r111705);
double r111707 = r111701 * r111706;
double r111708 = r111692 + r111707;
double r111709 = sqrt(r111708);
double r111710 = 1.0;
double r111711 = log1p(r111700);
double r111712 = expm1(r111711);
double r111713 = r111695 * r111712;
double r111714 = r111713 * r111703;
double r111715 = r111692 + r111714;
double r111716 = r111710 - r111715;
double r111717 = sqrt(r111716);
double r111718 = atan2(r111709, r111717);
double r111719 = r111686 * r111718;
double r111720 = r111685 * r111719;
return r111720;
}



Bits error versus R



Bits error versus lambda1



Bits error versus lambda2



Bits error versus phi1



Bits error versus phi2
Results
Initial program 24.1
rmApplied add-cbrt-cube24.3
Simplified24.3
rmApplied expm1-log1p-u24.3
rmApplied log1p-expm1-u24.3
rmApplied log1p-expm1-u24.3
Final simplification24.3
herbie shell --seed 2020057 +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))))))))))