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 \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 r106244 = R;
double r106245 = 2.0;
double r106246 = phi1;
double r106247 = phi2;
double r106248 = r106246 - r106247;
double r106249 = r106248 / r106245;
double r106250 = sin(r106249);
double r106251 = pow(r106250, r106245);
double r106252 = cos(r106246);
double r106253 = cos(r106247);
double r106254 = r106252 * r106253;
double r106255 = lambda1;
double r106256 = lambda2;
double r106257 = r106255 - r106256;
double r106258 = r106257 / r106245;
double r106259 = sin(r106258);
double r106260 = r106254 * r106259;
double r106261 = r106260 * r106259;
double r106262 = r106251 + r106261;
double r106263 = sqrt(r106262);
double r106264 = 1.0;
double r106265 = r106264 - r106262;
double r106266 = sqrt(r106265);
double r106267 = atan2(r106263, r106266);
double r106268 = r106245 * r106267;
double r106269 = r106244 * r106268;
return r106269;
}
double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r106270 = R;
double r106271 = 2.0;
double r106272 = phi1;
double r106273 = r106272 / r106271;
double r106274 = sin(r106273);
double r106275 = phi2;
double r106276 = r106275 / r106271;
double r106277 = cos(r106276);
double r106278 = r106274 * r106277;
double r106279 = cos(r106273);
double r106280 = sin(r106276);
double r106281 = r106279 * r106280;
double r106282 = r106278 - r106281;
double r106283 = pow(r106282, r106271);
double r106284 = cos(r106272);
double r106285 = cos(r106275);
double r106286 = r106284 * r106285;
double r106287 = lambda1;
double r106288 = lambda2;
double r106289 = r106287 - r106288;
double r106290 = r106289 / r106271;
double r106291 = sin(r106290);
double r106292 = r106286 * r106291;
double r106293 = r106292 * r106291;
double r106294 = r106283 + r106293;
double r106295 = sqrt(r106294);
double r106296 = 1.0;
double r106297 = expm1(r106291);
double r106298 = log1p(r106297);
double r106299 = r106286 * r106298;
double r106300 = 3.0;
double r106301 = pow(r106291, r106300);
double r106302 = cbrt(r106301);
double r106303 = r106299 * r106302;
double r106304 = r106283 + r106303;
double r106305 = r106296 - r106304;
double r106306 = sqrt(r106305);
double r106307 = atan2(r106295, r106306);
double r106308 = r106271 * r106307;
double r106309 = r106270 * r106308;
return r106309;
}



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.3
rmApplied log1p-expm1-u14.3
rmApplied add-cbrt-cube14.3
Simplified14.3
Final simplification14.3
herbie shell --seed 2020056 +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))))))))))