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 \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}{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 \sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)}}\right)double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r85944 = R;
double r85945 = 2.0;
double r85946 = phi1;
double r85947 = phi2;
double r85948 = r85946 - r85947;
double r85949 = r85948 / r85945;
double r85950 = sin(r85949);
double r85951 = pow(r85950, r85945);
double r85952 = cos(r85946);
double r85953 = cos(r85947);
double r85954 = r85952 * r85953;
double r85955 = lambda1;
double r85956 = lambda2;
double r85957 = r85955 - r85956;
double r85958 = r85957 / r85945;
double r85959 = sin(r85958);
double r85960 = r85954 * r85959;
double r85961 = r85960 * r85959;
double r85962 = r85951 + r85961;
double r85963 = sqrt(r85962);
double r85964 = 1.0;
double r85965 = r85964 - r85962;
double r85966 = sqrt(r85965);
double r85967 = atan2(r85963, r85966);
double r85968 = r85945 * r85967;
double r85969 = r85944 * r85968;
return r85969;
}
double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r85970 = R;
double r85971 = 2.0;
double r85972 = phi1;
double r85973 = r85972 / r85971;
double r85974 = sin(r85973);
double r85975 = phi2;
double r85976 = r85975 / r85971;
double r85977 = cos(r85976);
double r85978 = r85974 * r85977;
double r85979 = cos(r85973);
double r85980 = sin(r85976);
double r85981 = r85979 * r85980;
double r85982 = r85978 - r85981;
double r85983 = pow(r85982, r85971);
double r85984 = cos(r85972);
double r85985 = cos(r85975);
double r85986 = r85984 * r85985;
double r85987 = lambda1;
double r85988 = lambda2;
double r85989 = r85987 - r85988;
double r85990 = r85989 / r85971;
double r85991 = sin(r85990);
double r85992 = r85986 * r85991;
double r85993 = log1p(r85991);
double r85994 = expm1(r85993);
double r85995 = r85992 * r85994;
double r85996 = r85983 + r85995;
double r85997 = sqrt(r85996);
double r85998 = 1.0;
double r85999 = exp(r85991);
double r86000 = log(r85999);
double r86001 = r85986 * r86000;
double r86002 = r86001 * r85991;
double r86003 = r85983 + r86002;
double r86004 = r85998 - r86003;
double r86005 = sqrt(r86004);
double r86006 = atan2(r85997, r86005);
double r86007 = r85971 * r86006;
double r86008 = r85970 * r86007;
return r86008;
}



Bits error versus R



Bits error versus lambda1



Bits error versus lambda2



Bits error versus phi1



Bits error versus phi2
Results
Initial program 24.0
rmApplied div-sub24.0
Applied sin-diff23.5
rmApplied div-sub23.5
Applied sin-diff13.9
rmApplied add-log-exp13.9
rmApplied expm1-log1p-u13.9
Final simplification13.9
herbie shell --seed 2020024 +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))))))))))