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 \log \left(e^{\sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)}\right)\right)}}\right)double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r113221 = R;
double r113222 = 2.0;
double r113223 = phi1;
double r113224 = phi2;
double r113225 = r113223 - r113224;
double r113226 = r113225 / r113222;
double r113227 = sin(r113226);
double r113228 = pow(r113227, r113222);
double r113229 = cos(r113223);
double r113230 = cos(r113224);
double r113231 = r113229 * r113230;
double r113232 = lambda1;
double r113233 = lambda2;
double r113234 = r113232 - r113233;
double r113235 = r113234 / r113222;
double r113236 = sin(r113235);
double r113237 = r113231 * r113236;
double r113238 = r113237 * r113236;
double r113239 = r113228 + r113238;
double r113240 = sqrt(r113239);
double r113241 = 1.0;
double r113242 = r113241 - r113239;
double r113243 = sqrt(r113242);
double r113244 = atan2(r113240, r113243);
double r113245 = r113222 * r113244;
double r113246 = r113221 * r113245;
return r113246;
}
double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r113247 = R;
double r113248 = 2.0;
double r113249 = phi1;
double r113250 = r113249 / r113248;
double r113251 = sin(r113250);
double r113252 = phi2;
double r113253 = r113252 / r113248;
double r113254 = cos(r113253);
double r113255 = r113251 * r113254;
double r113256 = cos(r113250);
double r113257 = sin(r113253);
double r113258 = r113256 * r113257;
double r113259 = r113255 - r113258;
double r113260 = pow(r113259, r113248);
double r113261 = cos(r113249);
double r113262 = cos(r113252);
double r113263 = r113261 * r113262;
double r113264 = lambda1;
double r113265 = lambda2;
double r113266 = r113264 - r113265;
double r113267 = r113266 / r113248;
double r113268 = sin(r113267);
double r113269 = expm1(r113268);
double r113270 = log1p(r113269);
double r113271 = r113263 * r113270;
double r113272 = r113271 * r113268;
double r113273 = r113260 + r113272;
double r113274 = sqrt(r113273);
double r113275 = 1.0;
double r113276 = r113263 * r113268;
double r113277 = exp(r113268);
double r113278 = log(r113277);
double r113279 = r113276 * r113278;
double r113280 = r113260 + r113279;
double r113281 = r113275 - r113280;
double r113282 = sqrt(r113281);
double r113283 = atan2(r113274, r113282);
double r113284 = r113248 * r113283;
double r113285 = r113247 * r113284;
return r113285;
}



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-diff14.0
rmApplied add-log-exp14.0
rmApplied log1p-expm1-u14.0
Final simplification14.0
herbie shell --seed 2020033 +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))))))))))