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 \mathsf{log1p}\left(\mathsf{expm1}\left(\sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\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 r94347 = R;
double r94348 = 2.0;
double r94349 = phi1;
double r94350 = phi2;
double r94351 = r94349 - r94350;
double r94352 = r94351 / r94348;
double r94353 = sin(r94352);
double r94354 = pow(r94353, r94348);
double r94355 = cos(r94349);
double r94356 = cos(r94350);
double r94357 = r94355 * r94356;
double r94358 = lambda1;
double r94359 = lambda2;
double r94360 = r94358 - r94359;
double r94361 = r94360 / r94348;
double r94362 = sin(r94361);
double r94363 = r94357 * r94362;
double r94364 = r94363 * r94362;
double r94365 = r94354 + r94364;
double r94366 = sqrt(r94365);
double r94367 = 1.0;
double r94368 = r94367 - r94365;
double r94369 = sqrt(r94368);
double r94370 = atan2(r94366, r94369);
double r94371 = r94348 * r94370;
double r94372 = r94347 * r94371;
return r94372;
}
double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r94373 = R;
double r94374 = 2.0;
double r94375 = phi1;
double r94376 = r94375 / r94374;
double r94377 = sin(r94376);
double r94378 = phi2;
double r94379 = r94378 / r94374;
double r94380 = cos(r94379);
double r94381 = r94377 * r94380;
double r94382 = cos(r94376);
double r94383 = sin(r94379);
double r94384 = r94382 * r94383;
double r94385 = r94381 - r94384;
double r94386 = pow(r94385, r94374);
double r94387 = cos(r94375);
double r94388 = cos(r94378);
double r94389 = r94387 * r94388;
double r94390 = lambda1;
double r94391 = lambda2;
double r94392 = r94390 - r94391;
double r94393 = r94392 / r94374;
double r94394 = sin(r94393);
double r94395 = r94389 * r94394;
double r94396 = r94395 * r94394;
double r94397 = r94386 + r94396;
double r94398 = sqrt(r94397);
double r94399 = 1.0;
double r94400 = expm1(r94394);
double r94401 = log1p(r94400);
double r94402 = r94389 * r94401;
double r94403 = r94402 * r94401;
double r94404 = r94386 + r94403;
double r94405 = r94399 - r94404;
double r94406 = sqrt(r94405);
double r94407 = atan2(r94398, r94406);
double r94408 = r94374 * r94407;
double r94409 = r94373 * r94408;
return r94409;
}



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 log1p-expm1-u14.0
rmApplied log1p-expm1-u14.0
Final simplification14.0
herbie shell --seed 2020042 +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))))))))))