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 \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}{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 \log \left(e^{\sin \left(\frac{\lambda_1 - \lambda_2}{2}\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 r74878 = R;
double r74879 = 2.0;
double r74880 = phi1;
double r74881 = phi2;
double r74882 = r74880 - r74881;
double r74883 = r74882 / r74879;
double r74884 = sin(r74883);
double r74885 = pow(r74884, r74879);
double r74886 = cos(r74880);
double r74887 = cos(r74881);
double r74888 = r74886 * r74887;
double r74889 = lambda1;
double r74890 = lambda2;
double r74891 = r74889 - r74890;
double r74892 = r74891 / r74879;
double r74893 = sin(r74892);
double r74894 = r74888 * r74893;
double r74895 = r74894 * r74893;
double r74896 = r74885 + r74895;
double r74897 = sqrt(r74896);
double r74898 = 1.0;
double r74899 = r74898 - r74896;
double r74900 = sqrt(r74899);
double r74901 = atan2(r74897, r74900);
double r74902 = r74879 * r74901;
double r74903 = r74878 * r74902;
return r74903;
}
double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r74904 = R;
double r74905 = 2.0;
double r74906 = phi1;
double r74907 = r74906 / r74905;
double r74908 = sin(r74907);
double r74909 = phi2;
double r74910 = r74909 / r74905;
double r74911 = cos(r74910);
double r74912 = r74908 * r74911;
double r74913 = cos(r74907);
double r74914 = sin(r74910);
double r74915 = r74913 * r74914;
double r74916 = r74912 - r74915;
double r74917 = pow(r74916, r74905);
double r74918 = cos(r74906);
double r74919 = cos(r74909);
double r74920 = r74918 * r74919;
double r74921 = lambda1;
double r74922 = lambda2;
double r74923 = r74921 - r74922;
double r74924 = r74923 / r74905;
double r74925 = sin(r74924);
double r74926 = r74920 * r74925;
double r74927 = r74926 * r74925;
double r74928 = r74917 + r74927;
double r74929 = sqrt(r74928);
double r74930 = 1.0;
double r74931 = exp(r74925);
double r74932 = log(r74931);
double r74933 = r74920 * r74932;
double r74934 = expm1(r74925);
double r74935 = log1p(r74934);
double r74936 = r74933 * r74935;
double r74937 = r74917 + r74936;
double r74938 = r74930 - r74937;
double r74939 = sqrt(r74938);
double r74940 = atan2(r74929, r74939);
double r74941 = r74905 * r74940;
double r74942 = r74904 * r74941;
return r74942;
}



Bits error versus R



Bits error versus lambda1



Bits error versus lambda2



Bits error versus phi1



Bits error versus phi2
Results
Initial program 24.9
rmApplied div-sub24.9
Applied sin-diff24.2
rmApplied div-sub24.2
Applied sin-diff14.3
rmApplied add-log-exp14.3
rmApplied log1p-expm1-u14.3
Final simplification14.3
herbie shell --seed 2020025 +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))))))))))