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{log1p}\left(\mathsf{expm1}\left(\mathsf{expm1}\left(\mathsf{log1p}\left(\sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)\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 \sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)}}\right)double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r72903 = R;
double r72904 = 2.0;
double r72905 = phi1;
double r72906 = phi2;
double r72907 = r72905 - r72906;
double r72908 = r72907 / r72904;
double r72909 = sin(r72908);
double r72910 = pow(r72909, r72904);
double r72911 = cos(r72905);
double r72912 = cos(r72906);
double r72913 = r72911 * r72912;
double r72914 = lambda1;
double r72915 = lambda2;
double r72916 = r72914 - r72915;
double r72917 = r72916 / r72904;
double r72918 = sin(r72917);
double r72919 = r72913 * r72918;
double r72920 = r72919 * r72918;
double r72921 = r72910 + r72920;
double r72922 = sqrt(r72921);
double r72923 = 1.0;
double r72924 = r72923 - r72921;
double r72925 = sqrt(r72924);
double r72926 = atan2(r72922, r72925);
double r72927 = r72904 * r72926;
double r72928 = r72903 * r72927;
return r72928;
}
double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r72929 = R;
double r72930 = 2.0;
double r72931 = phi1;
double r72932 = phi2;
double r72933 = r72931 - r72932;
double r72934 = r72933 / r72930;
double r72935 = sin(r72934);
double r72936 = pow(r72935, r72930);
double r72937 = cos(r72931);
double r72938 = cos(r72932);
double r72939 = r72937 * r72938;
double r72940 = lambda1;
double r72941 = lambda2;
double r72942 = r72940 - r72941;
double r72943 = r72942 / r72930;
double r72944 = sin(r72943);
double r72945 = log1p(r72944);
double r72946 = expm1(r72945);
double r72947 = expm1(r72946);
double r72948 = log1p(r72947);
double r72949 = r72939 * r72948;
double r72950 = r72949 * r72946;
double r72951 = r72936 + r72950;
double r72952 = sqrt(r72951);
double r72953 = 1.0;
double r72954 = r72939 * r72944;
double r72955 = r72954 * r72944;
double r72956 = r72936 + r72955;
double r72957 = r72953 - r72956;
double r72958 = sqrt(r72957);
double r72959 = atan2(r72952, r72958);
double r72960 = r72930 * r72959;
double r72961 = r72929 * r72960;
return r72961;
}



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 log1p-expm1-u24.7
rmApplied expm1-log1p-u24.7
rmApplied expm1-log1p-u24.7
Final simplification24.7
herbie shell --seed 2020089 +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))))))))))