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{log1p}\left(\mathsf{expm1}\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 \sqrt[3]{{\left(\sin \left(\frac{\lambda_1 - \lambda_2}{2}\right)\right)}^{3}}\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 r75267 = R;
double r75268 = 2.0;
double r75269 = phi1;
double r75270 = phi2;
double r75271 = r75269 - r75270;
double r75272 = r75271 / r75268;
double r75273 = sin(r75272);
double r75274 = pow(r75273, r75268);
double r75275 = cos(r75269);
double r75276 = cos(r75270);
double r75277 = r75275 * r75276;
double r75278 = lambda1;
double r75279 = lambda2;
double r75280 = r75278 - r75279;
double r75281 = r75280 / r75268;
double r75282 = sin(r75281);
double r75283 = r75277 * r75282;
double r75284 = r75283 * r75282;
double r75285 = r75274 + r75284;
double r75286 = sqrt(r75285);
double r75287 = 1.0;
double r75288 = r75287 - r75285;
double r75289 = sqrt(r75288);
double r75290 = atan2(r75286, r75289);
double r75291 = r75268 * r75290;
double r75292 = r75267 * r75291;
return r75292;
}
double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r75293 = R;
double r75294 = 2.0;
double r75295 = phi1;
double r75296 = r75295 / r75294;
double r75297 = sin(r75296);
double r75298 = phi2;
double r75299 = r75298 / r75294;
double r75300 = cos(r75299);
double r75301 = r75297 * r75300;
double r75302 = cos(r75296);
double r75303 = sin(r75299);
double r75304 = r75302 * r75303;
double r75305 = r75301 - r75304;
double r75306 = pow(r75305, r75294);
double r75307 = cos(r75295);
double r75308 = cos(r75298);
double r75309 = r75307 * r75308;
double r75310 = lambda1;
double r75311 = lambda2;
double r75312 = r75310 - r75311;
double r75313 = r75312 / r75294;
double r75314 = sin(r75313);
double r75315 = r75309 * r75314;
double r75316 = expm1(r75314);
double r75317 = log1p(r75316);
double r75318 = r75315 * r75317;
double r75319 = r75306 + r75318;
double r75320 = sqrt(r75319);
double r75321 = 1.0;
double r75322 = 3.0;
double r75323 = pow(r75314, r75322);
double r75324 = cbrt(r75323);
double r75325 = r75309 * r75324;
double r75326 = r75325 * r75314;
double r75327 = r75306 + r75326;
double r75328 = r75321 - r75327;
double r75329 = sqrt(r75328);
double r75330 = atan2(r75320, r75329);
double r75331 = r75294 * r75330;
double r75332 = r75293 * r75331;
return r75332;
}



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 add-cbrt-cube14.6
Simplified14.6
rmApplied log1p-expm1-u14.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))))))))))