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 \mathsf{log1p}\left(\mathsf{expm1}\left(\sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)\right)\right) \cdot \sin \left(\frac{\lambda_1 - \lambda_2}{2}\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 \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 r80811 = R;
double r80812 = 2.0;
double r80813 = phi1;
double r80814 = phi2;
double r80815 = r80813 - r80814;
double r80816 = r80815 / r80812;
double r80817 = sin(r80816);
double r80818 = pow(r80817, r80812);
double r80819 = cos(r80813);
double r80820 = cos(r80814);
double r80821 = r80819 * r80820;
double r80822 = lambda1;
double r80823 = lambda2;
double r80824 = r80822 - r80823;
double r80825 = r80824 / r80812;
double r80826 = sin(r80825);
double r80827 = r80821 * r80826;
double r80828 = r80827 * r80826;
double r80829 = r80818 + r80828;
double r80830 = sqrt(r80829);
double r80831 = 1.0;
double r80832 = r80831 - r80829;
double r80833 = sqrt(r80832);
double r80834 = atan2(r80830, r80833);
double r80835 = r80812 * r80834;
double r80836 = r80811 * r80835;
return r80836;
}
double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r80837 = R;
double r80838 = 2.0;
double r80839 = phi1;
double r80840 = r80839 / r80838;
double r80841 = sin(r80840);
double r80842 = phi2;
double r80843 = r80842 / r80838;
double r80844 = cos(r80843);
double r80845 = r80841 * r80844;
double r80846 = cos(r80840);
double r80847 = sin(r80843);
double r80848 = r80846 * r80847;
double r80849 = r80845 - r80848;
double r80850 = pow(r80849, r80838);
double r80851 = cos(r80839);
double r80852 = cos(r80842);
double r80853 = r80851 * r80852;
double r80854 = lambda1;
double r80855 = lambda2;
double r80856 = r80854 - r80855;
double r80857 = r80856 / r80838;
double r80858 = sin(r80857);
double r80859 = expm1(r80858);
double r80860 = log1p(r80859);
double r80861 = r80853 * r80860;
double r80862 = r80861 * r80858;
double r80863 = r80850 + r80862;
double r80864 = sqrt(r80863);
double r80865 = 1.0;
double r80866 = r80853 * r80858;
double r80867 = r80866 * r80860;
double r80868 = r80850 + r80867;
double r80869 = r80865 - r80868;
double r80870 = sqrt(r80869);
double r80871 = atan2(r80864, r80870);
double r80872 = r80838 * r80871;
double r80873 = r80837 * r80872;
return r80873;
}



Bits error versus R



Bits error versus lambda1



Bits error versus lambda2



Bits error versus phi1



Bits error versus phi2
Results
Initial program 24.3
rmApplied div-sub24.3
Applied sin-diff23.7
rmApplied div-sub23.7
Applied sin-diff13.8
rmApplied log1p-expm1-u13.8
rmApplied log1p-expm1-u13.9
Final simplification13.9
herbie shell --seed 2020036 +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))))))))))