R \cdot \sqrt{\left(\left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\frac{\phi_1 + \phi_2}{2}\right)\right) \cdot \left(\left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\frac{\phi_1 + \phi_2}{2}\right)\right) + \left(\phi_1 - \phi_2\right) \cdot \left(\phi_1 - \phi_2\right)}\begin{array}{l}
\mathbf{if}\;\phi_1 \le -7.2710650707068061 \cdot 10^{71}:\\
\;\;\;\;R \cdot \left(\phi_2 - \phi_1\right)\\
\mathbf{else}:\\
\;\;\;\;R \cdot \sqrt{\left(\left(\lambda_1 - \lambda_2\right) \cdot \sqrt[3]{{\left(\cos \left(\frac{\phi_1 + \phi_2}{2}\right)\right)}^{3}}\right) \cdot \left(\left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\frac{\phi_1 + \phi_2}{2}\right)\right) + \left(\phi_1 - \phi_2\right) \cdot \left(\phi_1 - \phi_2\right)}\\
\end{array}double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r70928 = R;
double r70929 = lambda1;
double r70930 = lambda2;
double r70931 = r70929 - r70930;
double r70932 = phi1;
double r70933 = phi2;
double r70934 = r70932 + r70933;
double r70935 = 2.0;
double r70936 = r70934 / r70935;
double r70937 = cos(r70936);
double r70938 = r70931 * r70937;
double r70939 = r70938 * r70938;
double r70940 = r70932 - r70933;
double r70941 = r70940 * r70940;
double r70942 = r70939 + r70941;
double r70943 = sqrt(r70942);
double r70944 = r70928 * r70943;
return r70944;
}
double f(double R, double lambda1, double lambda2, double phi1, double phi2) {
double r70945 = phi1;
double r70946 = -7.271065070706806e+71;
bool r70947 = r70945 <= r70946;
double r70948 = R;
double r70949 = phi2;
double r70950 = r70949 - r70945;
double r70951 = r70948 * r70950;
double r70952 = lambda1;
double r70953 = lambda2;
double r70954 = r70952 - r70953;
double r70955 = r70945 + r70949;
double r70956 = 2.0;
double r70957 = r70955 / r70956;
double r70958 = cos(r70957);
double r70959 = 3.0;
double r70960 = pow(r70958, r70959);
double r70961 = cbrt(r70960);
double r70962 = r70954 * r70961;
double r70963 = r70954 * r70958;
double r70964 = r70962 * r70963;
double r70965 = r70945 - r70949;
double r70966 = r70965 * r70965;
double r70967 = r70964 + r70966;
double r70968 = sqrt(r70967);
double r70969 = r70948 * r70968;
double r70970 = r70947 ? r70951 : r70969;
return r70970;
}



Bits error versus R



Bits error versus lambda1



Bits error versus lambda2



Bits error versus phi1



Bits error versus phi2
Results
if phi1 < -7.271065070706806e+71Initial program 52.6
Taylor expanded around 0 21.0
if -7.271065070706806e+71 < phi1 Initial program 36.4
rmApplied add-cbrt-cube36.4
Simplified36.4
Final simplification33.5
herbie shell --seed 2020065
(FPCore (R lambda1 lambda2 phi1 phi2)
:name "Equirectangular approximation to distance on a great circle"
:precision binary64
(* R (sqrt (+ (* (* (- lambda1 lambda2) (cos (/ (+ phi1 phi2) 2))) (* (- lambda1 lambda2) (cos (/ (+ phi1 phi2) 2)))) (* (- phi1 phi2) (- phi1 phi2))))))