
(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 10 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 (<= phi1 -0.00033) (* (hypot phi1 (* (- lambda1 lambda2) (cos (* phi1 0.5)))) R) (* R (hypot phi2 (* (- lambda1 lambda2) (cos (* 0.5 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 (phi1 <= -0.00033) {
tmp = hypot(phi1, ((lambda1 - lambda2) * cos((phi1 * 0.5)))) * R;
} else {
tmp = R * hypot(phi2, ((lambda1 - lambda2) * cos((0.5 * phi2))));
}
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 (phi1 <= -0.00033) {
tmp = Math.hypot(phi1, ((lambda1 - lambda2) * Math.cos((phi1 * 0.5)))) * R;
} else {
tmp = R * Math.hypot(phi2, ((lambda1 - lambda2) * Math.cos((0.5 * phi2))));
}
return tmp;
}
[R, lambda1, lambda2, phi1, phi2] = sort([R, lambda1, lambda2, phi1, phi2]) def code(R, lambda1, lambda2, phi1, phi2): tmp = 0 if phi1 <= -0.00033: tmp = math.hypot(phi1, ((lambda1 - lambda2) * math.cos((phi1 * 0.5)))) * R else: tmp = R * math.hypot(phi2, ((lambda1 - lambda2) * math.cos((0.5 * 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 (phi1 <= -0.00033) tmp = Float64(hypot(phi1, Float64(Float64(lambda1 - lambda2) * cos(Float64(phi1 * 0.5)))) * R); else tmp = Float64(R * hypot(phi2, Float64(Float64(lambda1 - lambda2) * cos(Float64(0.5 * phi2))))); 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 (phi1 <= -0.00033)
tmp = hypot(phi1, ((lambda1 - lambda2) * cos((phi1 * 0.5)))) * R;
else
tmp = R * hypot(phi2, ((lambda1 - lambda2) * cos((0.5 * phi2))));
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[phi1, -0.00033], N[(N[Sqrt[phi1 ^ 2 + N[(N[(lambda1 - lambda2), $MachinePrecision] * N[Cos[N[(phi1 * 0.5), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] ^ 2], $MachinePrecision] * R), $MachinePrecision], N[(R * N[Sqrt[phi2 ^ 2 + N[(N[(lambda1 - lambda2), $MachinePrecision] * N[Cos[N[(0.5 * phi2), $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_1 \leq -0.00033:\\
\;\;\;\;\mathsf{hypot}\left(\phi_1, \left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\phi_1 \cdot 0.5\right)\right) \cdot R\\
\mathbf{else}:\\
\;\;\;\;R \cdot \mathsf{hypot}\left(\phi_2, \left(\lambda_1 - \lambda_2\right) \cdot \cos \left(0.5 \cdot \phi_2\right)\right)\\
\end{array}
\end{array}
if phi1 < -3.3e-4Initial program 53.7%
Taylor expanded in phi2 around 0
*-commutativeN/A
lower-*.f64N/A
+-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-*.f6490.5
Simplified90.5%
if -3.3e-4 < phi1 Initial program 64.5%
Taylor expanded in phi1 around 0
*-commutativeN/A
lower-*.f64N/A
+-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--.f6496.9
Simplified96.9%
Final simplification94.0%
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 4.6e+84) (* (hypot phi1 (* (- lambda1 lambda2) (cos (* phi1 0.5)))) R) (* R (hypot phi2 (* lambda2 (cos (* 0.5 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 <= 4.6e+84) {
tmp = hypot(phi1, ((lambda1 - lambda2) * cos((phi1 * 0.5)))) * R;
} else {
tmp = R * hypot(phi2, (lambda2 * cos((0.5 * phi2))));
}
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 <= 4.6e+84) {
tmp = Math.hypot(phi1, ((lambda1 - lambda2) * Math.cos((phi1 * 0.5)))) * R;
} else {
tmp = R * Math.hypot(phi2, (lambda2 * Math.cos((0.5 * phi2))));
}
return tmp;
}
[R, lambda1, lambda2, phi1, phi2] = sort([R, lambda1, lambda2, phi1, phi2]) def code(R, lambda1, lambda2, phi1, phi2): tmp = 0 if phi2 <= 4.6e+84: tmp = math.hypot(phi1, ((lambda1 - lambda2) * math.cos((phi1 * 0.5)))) * R else: tmp = R * math.hypot(phi2, (lambda2 * math.cos((0.5 * 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 <= 4.6e+84) tmp = Float64(hypot(phi1, Float64(Float64(lambda1 - lambda2) * cos(Float64(phi1 * 0.5)))) * R); else tmp = Float64(R * hypot(phi2, Float64(lambda2 * cos(Float64(0.5 * phi2))))); 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 <= 4.6e+84)
tmp = hypot(phi1, ((lambda1 - lambda2) * cos((phi1 * 0.5)))) * R;
else
tmp = R * hypot(phi2, (lambda2 * cos((0.5 * phi2))));
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, 4.6e+84], N[(N[Sqrt[phi1 ^ 2 + N[(N[(lambda1 - lambda2), $MachinePrecision] * N[Cos[N[(phi1 * 0.5), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] ^ 2], $MachinePrecision] * R), $MachinePrecision], N[(R * N[Sqrt[phi2 ^ 2 + N[(lambda2 * N[Cos[N[(0.5 * phi2), $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 4.6 \cdot 10^{+84}:\\
\;\;\;\;\mathsf{hypot}\left(\phi_1, \left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\phi_1 \cdot 0.5\right)\right) \cdot R\\
\mathbf{else}:\\
\;\;\;\;R \cdot \mathsf{hypot}\left(\phi_2, \lambda_2 \cdot \cos \left(0.5 \cdot \phi_2\right)\right)\\
\end{array}
\end{array}
if phi2 < 4.5999999999999998e84Initial program 65.0%
Taylor expanded in phi2 around 0
*-commutativeN/A
lower-*.f64N/A
+-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-*.f6487.4
Simplified87.4%
if 4.5999999999999998e84 < phi2 Initial program 52.8%
Taylor expanded in phi1 around 0
*-commutativeN/A
lower-*.f64N/A
+-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.7
Simplified88.7%
Taylor expanded in lambda1 around 0
+-commutativeN/A
unpow2N/A
unpow2N/A
unpow2N/A
unswap-sqrN/A
lower-hypot.f64N/A
*-commutativeN/A
lower-*.f64N/A
lower-cos.f64N/A
lower-*.f6482.3
Simplified82.3%
Final simplification85.7%
herbie shell --seed 2024218
(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))))))