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 \sqrt[3]{{\left(\sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)}^{3}}\right)}}\right)double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r108948 = R;
double r108949 = 2.0;
double r108950 = phi1;
double r108951 = phi2;
double r108952 = r108950 - r108951;
double r108953 = r108952 / r108949;
double r108954 = sin(r108953);
double r108955 = pow(r108954, r108949);
double r108956 = cos(r108950);
double r108957 = cos(r108951);
double r108958 = r108956 * r108957;
double r108959 = lambda1;
double r108960 = lambda2;
double r108961 = r108959 - r108960;
double r108962 = r108961 / r108949;
double r108963 = sin(r108962);
double r108964 = r108958 * r108963;
double r108965 = r108964 * r108963;
double r108966 = r108955 + r108965;
double r108967 = sqrt(r108966);
double r108968 = 1.0;
double r108969 = r108968 - r108966;
double r108970 = sqrt(r108969);
double r108971 = atan2(r108967, r108970);
double r108972 = r108949 * r108971;
double r108973 = r108948 * r108972;
return r108973;
}
double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r108974 = R;
double r108975 = 2.0;
double r108976 = phi1;
double r108977 = r108976 / r108975;
double r108978 = sin(r108977);
double r108979 = phi2;
double r108980 = r108979 / r108975;
double r108981 = cos(r108980);
double r108982 = r108978 * r108981;
double r108983 = cos(r108977);
double r108984 = sin(r108980);
double r108985 = r108983 * r108984;
double r108986 = r108982 - r108985;
double r108987 = pow(r108986, r108975);
double r108988 = cos(r108976);
double r108989 = cos(r108979);
double r108990 = r108988 * r108989;
double r108991 = lambda1;
double r108992 = lambda2;
double r108993 = r108991 - r108992;
double r108994 = r108993 / r108975;
double r108995 = sin(r108994);
double r108996 = expm1(r108995);
double r108997 = log1p(r108996);
double r108998 = r108990 * r108997;
double r108999 = r108998 * r108995;
double r109000 = r108987 + r108999;
double r109001 = sqrt(r109000);
double r109002 = 1.0;
double r109003 = r108990 * r108995;
double r109004 = 3.0;
double r109005 = pow(r108995, r109004);
double r109006 = cbrt(r109005);
double r109007 = r109003 * r109006;
double r109008 = r108987 + r109007;
double r109009 = r109002 - r109008;
double r109010 = sqrt(r109009);
double r109011 = atan2(r109001, r109010);
double r109012 = r108975 * r109011;
double r109013 = r108974 * r109012;
return r109013;
}



Bits error versus R



Bits error versus lambda1



Bits error versus lambda2



Bits error versus phi1



Bits error versus phi2
Results
Initial program 24.6
rmApplied div-sub24.6
Applied sin-diff24.0
rmApplied div-sub24.0
Applied sin-diff14.6
rmApplied log1p-expm1-u14.6
rmApplied add-cbrt-cube14.6
Simplified14.6
Final simplification14.6
herbie shell --seed 2019353 +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))))))))))