
(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 15 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 9.6e-30) (* 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 <= 9.6e-30) {
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 <= 9.6e-30) {
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 <= 9.6e-30: 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 <= 9.6e-30) 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 <= 9.6e-30)
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, 9.6e-30], 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 9.6 \cdot 10^{-30}:\\
\;\;\;\;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 < 9.5999999999999994e-30Initial program 69.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 9.5999999999999994e-30 < phi2 Initial program 60.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--.f6488.2
Applied rewrites88.2%
Final simplification92.9%
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 9.6e-30) (* R (hypot phi1 (* (- lambda1 lambda2) (cos (* phi1 0.5))))) (* R (hypot phi2 (- lambda1 lambda2)))))
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 <= 9.6e-30) {
tmp = R * hypot(phi1, ((lambda1 - lambda2) * cos((phi1 * 0.5))));
} else {
tmp = R * hypot(phi2, (lambda1 - lambda2));
}
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 <= 9.6e-30) {
tmp = R * Math.hypot(phi1, ((lambda1 - lambda2) * Math.cos((phi1 * 0.5))));
} else {
tmp = R * Math.hypot(phi2, (lambda1 - lambda2));
}
return tmp;
}
[R, lambda1, lambda2, phi1, phi2] = sort([R, lambda1, lambda2, phi1, phi2]) def code(R, lambda1, lambda2, phi1, phi2): tmp = 0 if phi2 <= 9.6e-30: tmp = R * math.hypot(phi1, ((lambda1 - lambda2) * math.cos((phi1 * 0.5)))) else: tmp = R * math.hypot(phi2, (lambda1 - lambda2)) 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 <= 9.6e-30) tmp = Float64(R * hypot(phi1, Float64(Float64(lambda1 - lambda2) * cos(Float64(phi1 * 0.5))))); else tmp = Float64(R * hypot(phi2, Float64(lambda1 - lambda2))); 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 <= 9.6e-30)
tmp = R * hypot(phi1, ((lambda1 - lambda2) * cos((phi1 * 0.5))));
else
tmp = R * hypot(phi2, (lambda1 - lambda2));
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, 9.6e-30], 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[(lambda1 - lambda2), $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 9.6 \cdot 10^{-30}:\\
\;\;\;\;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, \lambda_1 - \lambda_2\right)\\
\end{array}
\end{array}
if phi2 < 9.5999999999999994e-30Initial program 63.2%
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-*.f6496.5
Applied rewrites96.5%
if 9.5999999999999994e-30 < phi2 Initial program 54.9%
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.6
Applied rewrites84.6%
Taylor expanded in phi2 around 0
Applied rewrites73.3%
Final simplification85.5%
herbie shell --seed 2024230
(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))))))