
(FPCore (R lambda1 lambda2 phi1 phi2) :precision binary64 (let* ((t_0 (* (- lambda1 lambda2) (cos (/ (+ phi1 phi2) 2.0))))) (* R (sqrt (+ (* t_0 t_0) (* (- phi1 phi2) (- phi1 phi2)))))))
double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
double t_0 = (lambda1 - lambda2) * cos(((phi1 + phi2) / 2.0));
return R * sqrt(((t_0 * t_0) + ((phi1 - phi2) * (phi1 - phi2))));
}
real(8) function code(r, lambda1, lambda2, phi1, phi2)
real(8), intent (in) :: r
real(8), intent (in) :: lambda1
real(8), intent (in) :: lambda2
real(8), intent (in) :: phi1
real(8), intent (in) :: phi2
real(8) :: t_0
t_0 = (lambda1 - lambda2) * cos(((phi1 + phi2) / 2.0d0))
code = r * sqrt(((t_0 * t_0) + ((phi1 - phi2) * (phi1 - phi2))))
end function
public static double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
double t_0 = (lambda1 - lambda2) * Math.cos(((phi1 + phi2) / 2.0));
return R * Math.sqrt(((t_0 * t_0) + ((phi1 - phi2) * (phi1 - phi2))));
}
def code(R, lambda1, lambda2, phi1, phi2): t_0 = (lambda1 - lambda2) * math.cos(((phi1 + phi2) / 2.0)) return R * math.sqrt(((t_0 * t_0) + ((phi1 - phi2) * (phi1 - phi2))))
function code(R, lambda1, lambda2, phi1, phi2) t_0 = Float64(Float64(lambda1 - lambda2) * cos(Float64(Float64(phi1 + phi2) / 2.0))) return Float64(R * sqrt(Float64(Float64(t_0 * t_0) + Float64(Float64(phi1 - phi2) * Float64(phi1 - phi2))))) end
function tmp = code(R, lambda1, lambda2, phi1, phi2) t_0 = (lambda1 - lambda2) * cos(((phi1 + phi2) / 2.0)); tmp = R * sqrt(((t_0 * t_0) + ((phi1 - phi2) * (phi1 - phi2)))); end
code[R_, lambda1_, lambda2_, phi1_, phi2_] := Block[{t$95$0 = N[(N[(lambda1 - lambda2), $MachinePrecision] * N[Cos[N[(N[(phi1 + phi2), $MachinePrecision] / 2.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]}, N[(R * N[Sqrt[N[(N[(t$95$0 * t$95$0), $MachinePrecision] + N[(N[(phi1 - phi2), $MachinePrecision] * N[(phi1 - phi2), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\frac{\phi_1 + \phi_2}{2}\right)\\
R \cdot \sqrt{t\_0 \cdot t\_0 + \left(\phi_1 - \phi_2\right) \cdot \left(\phi_1 - \phi_2\right)}
\end{array}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 12 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (R lambda1 lambda2 phi1 phi2) :precision binary64 (let* ((t_0 (* (- lambda1 lambda2) (cos (/ (+ phi1 phi2) 2.0))))) (* R (sqrt (+ (* t_0 t_0) (* (- phi1 phi2) (- phi1 phi2)))))))
double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
double t_0 = (lambda1 - lambda2) * cos(((phi1 + phi2) / 2.0));
return R * sqrt(((t_0 * t_0) + ((phi1 - phi2) * (phi1 - phi2))));
}
real(8) function code(r, lambda1, lambda2, phi1, phi2)
real(8), intent (in) :: r
real(8), intent (in) :: lambda1
real(8), intent (in) :: lambda2
real(8), intent (in) :: phi1
real(8), intent (in) :: phi2
real(8) :: t_0
t_0 = (lambda1 - lambda2) * cos(((phi1 + phi2) / 2.0d0))
code = r * sqrt(((t_0 * t_0) + ((phi1 - phi2) * (phi1 - phi2))))
end function
public static double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
double t_0 = (lambda1 - lambda2) * Math.cos(((phi1 + phi2) / 2.0));
return R * Math.sqrt(((t_0 * t_0) + ((phi1 - phi2) * (phi1 - phi2))));
}
def code(R, lambda1, lambda2, phi1, phi2): t_0 = (lambda1 - lambda2) * math.cos(((phi1 + phi2) / 2.0)) return R * math.sqrt(((t_0 * t_0) + ((phi1 - phi2) * (phi1 - phi2))))
function code(R, lambda1, lambda2, phi1, phi2) t_0 = Float64(Float64(lambda1 - lambda2) * cos(Float64(Float64(phi1 + phi2) / 2.0))) return Float64(R * sqrt(Float64(Float64(t_0 * t_0) + Float64(Float64(phi1 - phi2) * Float64(phi1 - phi2))))) end
function tmp = code(R, lambda1, lambda2, phi1, phi2) t_0 = (lambda1 - lambda2) * cos(((phi1 + phi2) / 2.0)); tmp = R * sqrt(((t_0 * t_0) + ((phi1 - phi2) * (phi1 - phi2)))); end
code[R_, lambda1_, lambda2_, phi1_, phi2_] := Block[{t$95$0 = N[(N[(lambda1 - lambda2), $MachinePrecision] * N[Cos[N[(N[(phi1 + phi2), $MachinePrecision] / 2.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]}, N[(R * N[Sqrt[N[(N[(t$95$0 * t$95$0), $MachinePrecision] + N[(N[(phi1 - phi2), $MachinePrecision] * N[(phi1 - phi2), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\frac{\phi_1 + \phi_2}{2}\right)\\
R \cdot \sqrt{t\_0 \cdot t\_0 + \left(\phi_1 - \phi_2\right) \cdot \left(\phi_1 - \phi_2\right)}
\end{array}
\end{array}
NOTE: R, lambda1, lambda2, phi1, and phi2 should be sorted in increasing order before calling this function. (FPCore (R lambda1 lambda2 phi1 phi2) :precision binary64 (if (<= phi2 0.029) (* R (hypot phi1 (* (- lambda1 lambda2) (cos (* phi1 0.5))))) (* R (hypot phi2 (* (- lambda1 lambda2) (cos (* phi2 0.5)))))))
assert(R < lambda1 && lambda1 < lambda2 && lambda2 < phi1 && phi1 < phi2);
double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
double tmp;
if (phi2 <= 0.029) {
tmp = R * hypot(phi1, ((lambda1 - lambda2) * cos((phi1 * 0.5))));
} else {
tmp = R * hypot(phi2, ((lambda1 - lambda2) * cos((phi2 * 0.5))));
}
return tmp;
}
assert R < lambda1 && lambda1 < lambda2 && lambda2 < phi1 && phi1 < phi2;
public static double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
double tmp;
if (phi2 <= 0.029) {
tmp = R * Math.hypot(phi1, ((lambda1 - lambda2) * Math.cos((phi1 * 0.5))));
} else {
tmp = R * Math.hypot(phi2, ((lambda1 - lambda2) * Math.cos((phi2 * 0.5))));
}
return tmp;
}
[R, lambda1, lambda2, phi1, phi2] = sort([R, lambda1, lambda2, phi1, phi2]) def code(R, lambda1, lambda2, phi1, phi2): tmp = 0 if phi2 <= 0.029: tmp = R * math.hypot(phi1, ((lambda1 - lambda2) * math.cos((phi1 * 0.5)))) else: tmp = R * math.hypot(phi2, ((lambda1 - lambda2) * math.cos((phi2 * 0.5)))) return tmp
R, lambda1, lambda2, phi1, phi2 = sort([R, lambda1, lambda2, phi1, phi2]) function code(R, lambda1, lambda2, phi1, phi2) tmp = 0.0 if (phi2 <= 0.029) tmp = Float64(R * hypot(phi1, Float64(Float64(lambda1 - lambda2) * cos(Float64(phi1 * 0.5))))); else tmp = Float64(R * hypot(phi2, Float64(Float64(lambda1 - lambda2) * cos(Float64(phi2 * 0.5))))); end return tmp end
R, lambda1, lambda2, phi1, phi2 = num2cell(sort([R, lambda1, lambda2, phi1, phi2])){:}
function tmp_2 = code(R, lambda1, lambda2, phi1, phi2)
tmp = 0.0;
if (phi2 <= 0.029)
tmp = R * hypot(phi1, ((lambda1 - lambda2) * cos((phi1 * 0.5))));
else
tmp = R * hypot(phi2, ((lambda1 - lambda2) * cos((phi2 * 0.5))));
end
tmp_2 = tmp;
end
NOTE: R, lambda1, lambda2, phi1, and phi2 should be sorted in increasing order before calling this function. code[R_, lambda1_, lambda2_, phi1_, phi2_] := If[LessEqual[phi2, 0.029], N[(R * N[Sqrt[phi1 ^ 2 + N[(N[(lambda1 - lambda2), $MachinePrecision] * N[Cos[N[(phi1 * 0.5), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] ^ 2], $MachinePrecision]), $MachinePrecision], N[(R * N[Sqrt[phi2 ^ 2 + N[(N[(lambda1 - lambda2), $MachinePrecision] * N[Cos[N[(phi2 * 0.5), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] ^ 2], $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
[R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\
\\
\begin{array}{l}
\mathbf{if}\;\phi_2 \leq 0.029:\\
\;\;\;\;R \cdot \mathsf{hypot}\left(\phi_1, \left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\phi_1 \cdot 0.5\right)\right)\\
\mathbf{else}:\\
\;\;\;\;R \cdot \mathsf{hypot}\left(\phi_2, \left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\phi_2 \cdot 0.5\right)\right)\\
\end{array}
\end{array}
if phi2 < 0.0290000000000000015Initial program 72.4%
Taylor expanded in phi2 around 0
+-commutativeN/A
unpow2N/A
unpow2N/A
unpow2N/A
unswap-sqrN/A
lower-hypot.f64N/A
*-commutativeN/A
lower-*.f64N/A
lower--.f64N/A
lower-cos.f64N/A
lower-*.f6497.1
Applied rewrites97.1%
if 0.0290000000000000015 < phi2 Initial program 53.3%
Taylor expanded in phi1 around 0
+-commutativeN/A
unpow2N/A
unpow2N/A
unpow2N/A
unswap-sqrN/A
lower-hypot.f64N/A
lower-*.f64N/A
lower-cos.f64N/A
lower-*.f64N/A
lower--.f6484.5
Applied rewrites84.5%
Final simplification91.6%
NOTE: R, lambda1, lambda2, phi1, and phi2 should be sorted in increasing order before calling this function. (FPCore (R lambda1 lambda2 phi1 phi2) :precision binary64 (if (<= phi2 2.25e+38) (* R (hypot phi1 (* (- lambda1 lambda2) (cos (* phi1 0.5))))) (* R (fma phi2 (/ phi1 (- phi2)) phi2))))
assert(R < lambda1 && lambda1 < lambda2 && lambda2 < phi1 && phi1 < phi2);
double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
double tmp;
if (phi2 <= 2.25e+38) {
tmp = R * hypot(phi1, ((lambda1 - lambda2) * cos((phi1 * 0.5))));
} else {
tmp = R * fma(phi2, (phi1 / -phi2), phi2);
}
return tmp;
}
R, lambda1, lambda2, phi1, phi2 = sort([R, lambda1, lambda2, phi1, phi2]) function code(R, lambda1, lambda2, phi1, phi2) tmp = 0.0 if (phi2 <= 2.25e+38) tmp = Float64(R * hypot(phi1, Float64(Float64(lambda1 - lambda2) * cos(Float64(phi1 * 0.5))))); else tmp = Float64(R * fma(phi2, Float64(phi1 / Float64(-phi2)), phi2)); end return tmp end
NOTE: R, lambda1, lambda2, phi1, and phi2 should be sorted in increasing order before calling this function. code[R_, lambda1_, lambda2_, phi1_, phi2_] := If[LessEqual[phi2, 2.25e+38], N[(R * N[Sqrt[phi1 ^ 2 + N[(N[(lambda1 - lambda2), $MachinePrecision] * N[Cos[N[(phi1 * 0.5), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] ^ 2], $MachinePrecision]), $MachinePrecision], N[(R * N[(phi2 * N[(phi1 / (-phi2)), $MachinePrecision] + phi2), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
[R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\
\\
\begin{array}{l}
\mathbf{if}\;\phi_2 \leq 2.25 \cdot 10^{+38}:\\
\;\;\;\;R \cdot \mathsf{hypot}\left(\phi_1, \left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\phi_1 \cdot 0.5\right)\right)\\
\mathbf{else}:\\
\;\;\;\;R \cdot \mathsf{fma}\left(\phi_2, \frac{\phi_1}{-\phi_2}, \phi_2\right)\\
\end{array}
\end{array}
if phi2 < 2.2499999999999999e38Initial program 63.3%
Taylor expanded in phi2 around 0
+-commutativeN/A
unpow2N/A
unpow2N/A
unpow2N/A
unswap-sqrN/A
lower-hypot.f64N/A
*-commutativeN/A
lower-*.f64N/A
lower--.f64N/A
lower-cos.f64N/A
lower-*.f6491.4
Applied rewrites91.4%
if 2.2499999999999999e38 < phi2 Initial program 51.8%
Taylor expanded in phi2 around inf
+-commutativeN/A
distribute-lft-inN/A
*-rgt-identityN/A
lower-fma.f64N/A
mul-1-negN/A
distribute-neg-frac2N/A
mul-1-negN/A
lower-/.f64N/A
mul-1-negN/A
lower-neg.f6475.4
Applied rewrites75.4%
Final simplification85.1%
herbie shell --seed 2024228
(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.0))) (* (- lambda1 lambda2) (cos (/ (+ phi1 phi2) 2.0)))) (* (- phi1 phi2) (- phi1 phi2))))))