Equirectangular approximation to distance on a great circle

Percentage Accurate: 59.3% → 90.6%
Time: 8.9s
Alternatives: 15
Speedup: 8.5×

Specification

?
\[\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} \]
(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))));
}
module fmin_fmax_functions
    implicit none
    private
    public fmax
    public fmin

    interface fmax
        module procedure fmax88
        module procedure fmax44
        module procedure fmax84
        module procedure fmax48
    end interface
    interface fmin
        module procedure fmin88
        module procedure fmin44
        module procedure fmin84
        module procedure fmin48
    end interface
contains
    real(8) function fmax88(x, y) result (res)
        real(8), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(x, max(x, y), y /= y), x /= x)
    end function
    real(4) function fmax44(x, y) result (res)
        real(4), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(y, merge(x, max(x, y), y /= y), x /= x)
    end function
    real(8) function fmax84(x, y) result(res)
        real(8), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
    end function
    real(8) function fmax48(x, y) result(res)
        real(4), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
    end function
    real(8) function fmin88(x, y) result (res)
        real(8), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(x, min(x, y), y /= y), x /= x)
    end function
    real(4) function fmin44(x, y) result (res)
        real(4), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(y, merge(x, min(x, y), y /= y), x /= x)
    end function
    real(8) function fmin84(x, y) result(res)
        real(8), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
    end function
    real(8) function fmin48(x, y) result(res)
        real(4), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
    end function
end module

real(8) function code(r, lambda1, lambda2, phi1, phi2)
use fmin_fmax_functions
    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:

Local Percentage Accuracy vs ?

The average percentage accuracy by input value. Horizontal axis shows value of an input variable; the variable is choosen in the title. Vertical axis is accuracy; higher is better. Red represent the original program, while blue represents Herbie's suggestion. These can be toggled with buttons below the plot. The line is an average while dots represent individual samples.

Accuracy vs Speed?

Herbie found 15 alternatives:

AlternativeAccuracySpeedup
The accuracy (vertical axis) and speed (horizontal axis) of each alternatives. Up and to the right is better. The red square shows the initial program, and each blue circle shows an alternative.The line shows the best available speed-accuracy tradeoffs.

Initial Program: 59.3% accurate, 1.0× speedup?

\[\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} \]
(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))));
}
module fmin_fmax_functions
    implicit none
    private
    public fmax
    public fmin

    interface fmax
        module procedure fmax88
        module procedure fmax44
        module procedure fmax84
        module procedure fmax48
    end interface
    interface fmin
        module procedure fmin88
        module procedure fmin44
        module procedure fmin84
        module procedure fmin48
    end interface
contains
    real(8) function fmax88(x, y) result (res)
        real(8), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(x, max(x, y), y /= y), x /= x)
    end function
    real(4) function fmax44(x, y) result (res)
        real(4), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(y, merge(x, max(x, y), y /= y), x /= x)
    end function
    real(8) function fmax84(x, y) result(res)
        real(8), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
    end function
    real(8) function fmax48(x, y) result(res)
        real(4), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
    end function
    real(8) function fmin88(x, y) result (res)
        real(8), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(x, min(x, y), y /= y), x /= x)
    end function
    real(4) function fmin44(x, y) result (res)
        real(4), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(y, merge(x, min(x, y), y /= y), x /= x)
    end function
    real(8) function fmin84(x, y) result(res)
        real(8), intent (in) :: x
        real(4), intent (in) :: y
        res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
    end function
    real(8) function fmin48(x, y) result(res)
        real(4), intent (in) :: x
        real(8), intent (in) :: y
        res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
    end function
end module

real(8) function code(r, lambda1, lambda2, phi1, phi2)
use fmin_fmax_functions
    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}

Alternative 1: 90.6% accurate, 1.2× speedup?

\[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \begin{array}{l} \mathbf{if}\;\phi_1 \leq -1 \cdot 10^{-34}:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_1\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_1\right)\\ \mathbf{else}:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_2\right) \cdot \left(\lambda_1 - \lambda_2\right), \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 -1e-34)
   (* R (hypot (* (cos (* -0.5 phi1)) (- lambda1 lambda2)) phi1))
   (* R (hypot (* (cos (* -0.5 phi2)) (- lambda1 lambda2)) 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 <= -1e-34) {
		tmp = R * hypot((cos((-0.5 * phi1)) * (lambda1 - lambda2)), phi1);
	} else {
		tmp = R * hypot((cos((-0.5 * phi2)) * (lambda1 - lambda2)), 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 <= -1e-34) {
		tmp = R * Math.hypot((Math.cos((-0.5 * phi1)) * (lambda1 - lambda2)), phi1);
	} else {
		tmp = R * Math.hypot((Math.cos((-0.5 * phi2)) * (lambda1 - lambda2)), 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 <= -1e-34:
		tmp = R * math.hypot((math.cos((-0.5 * phi1)) * (lambda1 - lambda2)), phi1)
	else:
		tmp = R * math.hypot((math.cos((-0.5 * phi2)) * (lambda1 - lambda2)), 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 <= -1e-34)
		tmp = Float64(R * hypot(Float64(cos(Float64(-0.5 * phi1)) * Float64(lambda1 - lambda2)), phi1));
	else
		tmp = Float64(R * hypot(Float64(cos(Float64(-0.5 * phi2)) * Float64(lambda1 - lambda2)), 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 <= -1e-34)
		tmp = R * hypot((cos((-0.5 * phi1)) * (lambda1 - lambda2)), phi1);
	else
		tmp = R * hypot((cos((-0.5 * phi2)) * (lambda1 - lambda2)), 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, -1e-34], N[(R * N[Sqrt[N[(N[Cos[N[(-0.5 * phi1), $MachinePrecision]], $MachinePrecision] * N[(lambda1 - lambda2), $MachinePrecision]), $MachinePrecision] ^ 2 + phi1 ^ 2], $MachinePrecision]), $MachinePrecision], N[(R * N[Sqrt[N[(N[Cos[N[(-0.5 * phi2), $MachinePrecision]], $MachinePrecision] * N[(lambda1 - lambda2), $MachinePrecision]), $MachinePrecision] ^ 2 + phi2 ^ 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 -1 \cdot 10^{-34}:\\
\;\;\;\;R \cdot \mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_1\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_1\right)\\

\mathbf{else}:\\
\;\;\;\;R \cdot \mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_2\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_2\right)\\


\end{array}
\end{array}
Derivation
  1. Split input into 2 regimes
  2. if phi1 < -9.99999999999999928e-35

    1. Initial program 58.0%

      \[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)} \]
    2. Add Preprocessing
    3. Taylor expanded in phi2 around 0

      \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_1\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_1}^{2}}} \]
    4. Step-by-step derivation
      1. Applied rewrites83.7%

        \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_1\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_1\right)} \]

      if -9.99999999999999928e-35 < phi1

      1. Initial program 64.6%

        \[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)} \]
      2. Add Preprocessing
      3. Taylor expanded in phi1 around 0

        \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_2\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_2}^{2}}} \]
      4. Step-by-step derivation
        1. Applied rewrites82.6%

          \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_2\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_2\right)} \]
      5. Recombined 2 regimes into one program.
      6. Add Preprocessing

      Alternative 2: 81.0% accurate, 1.2× speedup?

      \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \begin{array}{l} t_0 := \cos \left(0.5 \cdot \phi_1\right)\\ \mathbf{if}\;\phi_1 \leq -1.18 \cdot 10^{+107}:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(t\_0 \cdot \lambda_1, \phi_1\right)\\ \mathbf{elif}\;\phi_1 \leq -215000:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(t\_0 \cdot \lambda_2, \phi_1\right)\\ \mathbf{else}:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \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
       (let* ((t_0 (cos (* 0.5 phi1))))
         (if (<= phi1 -1.18e+107)
           (* R (hypot (* t_0 lambda1) phi1))
           (if (<= phi1 -215000.0)
             (* R (hypot (* t_0 lambda2) phi1))
             (* R (hypot (- lambda1 lambda2) phi2))))))
      assert(R < lambda1 && lambda1 < lambda2 && lambda2 < phi1 && phi1 < phi2);
      double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
      	double t_0 = cos((0.5 * phi1));
      	double tmp;
      	if (phi1 <= -1.18e+107) {
      		tmp = R * hypot((t_0 * lambda1), phi1);
      	} else if (phi1 <= -215000.0) {
      		tmp = R * hypot((t_0 * lambda2), phi1);
      	} else {
      		tmp = R * hypot((lambda1 - lambda2), 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 t_0 = Math.cos((0.5 * phi1));
      	double tmp;
      	if (phi1 <= -1.18e+107) {
      		tmp = R * Math.hypot((t_0 * lambda1), phi1);
      	} else if (phi1 <= -215000.0) {
      		tmp = R * Math.hypot((t_0 * lambda2), phi1);
      	} else {
      		tmp = R * Math.hypot((lambda1 - lambda2), phi2);
      	}
      	return tmp;
      }
      
      [R, lambda1, lambda2, phi1, phi2] = sort([R, lambda1, lambda2, phi1, phi2])
      def code(R, lambda1, lambda2, phi1, phi2):
      	t_0 = math.cos((0.5 * phi1))
      	tmp = 0
      	if phi1 <= -1.18e+107:
      		tmp = R * math.hypot((t_0 * lambda1), phi1)
      	elif phi1 <= -215000.0:
      		tmp = R * math.hypot((t_0 * lambda2), phi1)
      	else:
      		tmp = R * math.hypot((lambda1 - lambda2), phi2)
      	return tmp
      
      R, lambda1, lambda2, phi1, phi2 = sort([R, lambda1, lambda2, phi1, phi2])
      function code(R, lambda1, lambda2, phi1, phi2)
      	t_0 = cos(Float64(0.5 * phi1))
      	tmp = 0.0
      	if (phi1 <= -1.18e+107)
      		tmp = Float64(R * hypot(Float64(t_0 * lambda1), phi1));
      	elseif (phi1 <= -215000.0)
      		tmp = Float64(R * hypot(Float64(t_0 * lambda2), phi1));
      	else
      		tmp = Float64(R * hypot(Float64(lambda1 - lambda2), 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)
      	t_0 = cos((0.5 * phi1));
      	tmp = 0.0;
      	if (phi1 <= -1.18e+107)
      		tmp = R * hypot((t_0 * lambda1), phi1);
      	elseif (phi1 <= -215000.0)
      		tmp = R * hypot((t_0 * lambda2), phi1);
      	else
      		tmp = R * hypot((lambda1 - lambda2), 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_] := Block[{t$95$0 = N[Cos[N[(0.5 * phi1), $MachinePrecision]], $MachinePrecision]}, If[LessEqual[phi1, -1.18e+107], N[(R * N[Sqrt[N[(t$95$0 * lambda1), $MachinePrecision] ^ 2 + phi1 ^ 2], $MachinePrecision]), $MachinePrecision], If[LessEqual[phi1, -215000.0], N[(R * N[Sqrt[N[(t$95$0 * lambda2), $MachinePrecision] ^ 2 + phi1 ^ 2], $MachinePrecision]), $MachinePrecision], N[(R * N[Sqrt[N[(lambda1 - lambda2), $MachinePrecision] ^ 2 + phi2 ^ 2], $MachinePrecision]), $MachinePrecision]]]]
      
      \begin{array}{l}
      [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\
      \\
      \begin{array}{l}
      t_0 := \cos \left(0.5 \cdot \phi_1\right)\\
      \mathbf{if}\;\phi_1 \leq -1.18 \cdot 10^{+107}:\\
      \;\;\;\;R \cdot \mathsf{hypot}\left(t\_0 \cdot \lambda_1, \phi_1\right)\\
      
      \mathbf{elif}\;\phi_1 \leq -215000:\\
      \;\;\;\;R \cdot \mathsf{hypot}\left(t\_0 \cdot \lambda_2, \phi_1\right)\\
      
      \mathbf{else}:\\
      \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_2\right)\\
      
      
      \end{array}
      \end{array}
      
      Derivation
      1. Split input into 3 regimes
      2. if phi1 < -1.18000000000000005e107

        1. Initial program 45.1%

          \[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)} \]
        2. Add Preprocessing
        3. Taylor expanded in phi2 around 0

          \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_1\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_1}^{2}}} \]
        4. Step-by-step derivation
          1. Applied rewrites87.6%

            \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_1\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_1\right)} \]
          2. Taylor expanded in lambda1 around inf

            \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 \cdot \cos \left(\frac{-1}{2} \cdot \phi_1\right), \phi_1\right) \]
          3. Step-by-step derivation
            1. Applied rewrites82.0%

              \[\leadsto R \cdot \mathsf{hypot}\left(\cos \left(0.5 \cdot \phi_1\right) \cdot \lambda_1, \phi_1\right) \]

            if -1.18000000000000005e107 < phi1 < -215000

            1. Initial program 72.2%

              \[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)} \]
            2. Add Preprocessing
            3. Taylor expanded in phi2 around 0

              \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_1\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_1}^{2}}} \]
            4. Step-by-step derivation
              1. Applied rewrites74.0%

                \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_1\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_1\right)} \]
              2. Taylor expanded in lambda1 around 0

                \[\leadsto R \cdot \sqrt{{\lambda_2}^{2} \cdot {\cos \left(\frac{-1}{2} \cdot \phi_1\right)}^{2} + {\phi_1}^{2}} \]
              3. Step-by-step derivation
                1. Applied rewrites68.5%

                  \[\leadsto R \cdot \mathsf{hypot}\left(\cos \left(0.5 \cdot \phi_1\right) \cdot \lambda_2, \color{blue}{\phi_1}\right) \]

                if -215000 < phi1

                1. Initial program 65.9%

                  \[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)} \]
                2. Add Preprocessing
                3. Taylor expanded in phi1 around 0

                  \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_2\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_2}^{2}}} \]
                4. Step-by-step derivation
                  1. Applied rewrites82.6%

                    \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_2\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_2\right)} \]
                  2. Taylor expanded in phi2 around 0

                    \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_2\right) \]
                  3. Step-by-step derivation
                    1. Applied rewrites77.1%

                      \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_2\right) \]
                  4. Recombined 3 regimes into one program.
                  5. Add Preprocessing

                  Alternative 3: 86.0% accurate, 1.2× speedup?

                  \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \begin{array}{l} \mathbf{if}\;\phi_2 \leq 5.8 \cdot 10^{+68}:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_1\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_1\right)\\ \mathbf{else}:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(\cos \left(0.5 \cdot \phi_2\right) \cdot \lambda_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 5.8e+68)
                     (* R (hypot (* (cos (* -0.5 phi1)) (- lambda1 lambda2)) phi1))
                     (* R (hypot (* (cos (* 0.5 phi2)) lambda1) 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 <= 5.8e+68) {
                  		tmp = R * hypot((cos((-0.5 * phi1)) * (lambda1 - lambda2)), phi1);
                  	} else {
                  		tmp = R * hypot((cos((0.5 * phi2)) * lambda1), 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 <= 5.8e+68) {
                  		tmp = R * Math.hypot((Math.cos((-0.5 * phi1)) * (lambda1 - lambda2)), phi1);
                  	} else {
                  		tmp = R * Math.hypot((Math.cos((0.5 * phi2)) * lambda1), 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 <= 5.8e+68:
                  		tmp = R * math.hypot((math.cos((-0.5 * phi1)) * (lambda1 - lambda2)), phi1)
                  	else:
                  		tmp = R * math.hypot((math.cos((0.5 * phi2)) * lambda1), 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 <= 5.8e+68)
                  		tmp = Float64(R * hypot(Float64(cos(Float64(-0.5 * phi1)) * Float64(lambda1 - lambda2)), phi1));
                  	else
                  		tmp = Float64(R * hypot(Float64(cos(Float64(0.5 * phi2)) * lambda1), 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 <= 5.8e+68)
                  		tmp = R * hypot((cos((-0.5 * phi1)) * (lambda1 - lambda2)), phi1);
                  	else
                  		tmp = R * hypot((cos((0.5 * phi2)) * lambda1), 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, 5.8e+68], N[(R * N[Sqrt[N[(N[Cos[N[(-0.5 * phi1), $MachinePrecision]], $MachinePrecision] * N[(lambda1 - lambda2), $MachinePrecision]), $MachinePrecision] ^ 2 + phi1 ^ 2], $MachinePrecision]), $MachinePrecision], N[(R * N[Sqrt[N[(N[Cos[N[(0.5 * phi2), $MachinePrecision]], $MachinePrecision] * lambda1), $MachinePrecision] ^ 2 + phi2 ^ 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 5.8 \cdot 10^{+68}:\\
                  \;\;\;\;R \cdot \mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_1\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_1\right)\\
                  
                  \mathbf{else}:\\
                  \;\;\;\;R \cdot \mathsf{hypot}\left(\cos \left(0.5 \cdot \phi_2\right) \cdot \lambda_1, \phi_2\right)\\
                  
                  
                  \end{array}
                  \end{array}
                  
                  Derivation
                  1. Split input into 2 regimes
                  2. if phi2 < 5.80000000000000023e68

                    1. Initial program 64.1%

                      \[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)} \]
                    2. Add Preprocessing
                    3. Taylor expanded in phi2 around 0

                      \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_1\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_1}^{2}}} \]
                    4. Step-by-step derivation
                      1. Applied rewrites84.2%

                        \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_1\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_1\right)} \]

                      if 5.80000000000000023e68 < phi2

                      1. Initial program 56.5%

                        \[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)} \]
                      2. Add Preprocessing
                      3. Taylor expanded in phi1 around 0

                        \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_2\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_2}^{2}}} \]
                      4. Step-by-step derivation
                        1. Applied rewrites82.2%

                          \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_2\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_2\right)} \]
                        2. Taylor expanded in lambda1 around inf

                          \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 \cdot \cos \left(\frac{-1}{2} \cdot \phi_2\right), \phi_2\right) \]
                        3. Step-by-step derivation
                          1. Applied rewrites79.1%

                            \[\leadsto R \cdot \mathsf{hypot}\left(\cos \left(0.5 \cdot \phi_2\right) \cdot \lambda_1, \phi_2\right) \]
                        4. Recombined 2 regimes into one program.
                        5. Add Preprocessing

                        Alternative 4: 80.6% accurate, 1.3× speedup?

                        \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \begin{array}{l} \mathbf{if}\;\phi_2 \leq 5.8 \cdot 10^{+68}:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_1\right)\\ \mathbf{else}:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(\cos \left(0.5 \cdot \phi_2\right) \cdot \lambda_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 5.8e+68)
                           (* R (hypot (- lambda1 lambda2) phi1))
                           (* R (hypot (* (cos (* 0.5 phi2)) lambda1) 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 <= 5.8e+68) {
                        		tmp = R * hypot((lambda1 - lambda2), phi1);
                        	} else {
                        		tmp = R * hypot((cos((0.5 * phi2)) * lambda1), 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 <= 5.8e+68) {
                        		tmp = R * Math.hypot((lambda1 - lambda2), phi1);
                        	} else {
                        		tmp = R * Math.hypot((Math.cos((0.5 * phi2)) * lambda1), 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 <= 5.8e+68:
                        		tmp = R * math.hypot((lambda1 - lambda2), phi1)
                        	else:
                        		tmp = R * math.hypot((math.cos((0.5 * phi2)) * lambda1), 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 <= 5.8e+68)
                        		tmp = Float64(R * hypot(Float64(lambda1 - lambda2), phi1));
                        	else
                        		tmp = Float64(R * hypot(Float64(cos(Float64(0.5 * phi2)) * lambda1), 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 <= 5.8e+68)
                        		tmp = R * hypot((lambda1 - lambda2), phi1);
                        	else
                        		tmp = R * hypot((cos((0.5 * phi2)) * lambda1), 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, 5.8e+68], N[(R * N[Sqrt[N[(lambda1 - lambda2), $MachinePrecision] ^ 2 + phi1 ^ 2], $MachinePrecision]), $MachinePrecision], N[(R * N[Sqrt[N[(N[Cos[N[(0.5 * phi2), $MachinePrecision]], $MachinePrecision] * lambda1), $MachinePrecision] ^ 2 + phi2 ^ 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 5.8 \cdot 10^{+68}:\\
                        \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_1\right)\\
                        
                        \mathbf{else}:\\
                        \;\;\;\;R \cdot \mathsf{hypot}\left(\cos \left(0.5 \cdot \phi_2\right) \cdot \lambda_1, \phi_2\right)\\
                        
                        
                        \end{array}
                        \end{array}
                        
                        Derivation
                        1. Split input into 2 regimes
                        2. if phi2 < 5.80000000000000023e68

                          1. Initial program 64.1%

                            \[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)} \]
                          2. Add Preprocessing
                          3. Taylor expanded in phi2 around 0

                            \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_1\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_1}^{2}}} \]
                          4. Step-by-step derivation
                            1. Applied rewrites84.2%

                              \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_1\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_1\right)} \]
                            2. Taylor expanded in phi1 around 0

                              \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_1\right) \]
                            3. Step-by-step derivation
                              1. Applied rewrites78.4%

                                \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_1\right) \]

                              if 5.80000000000000023e68 < phi2

                              1. Initial program 56.5%

                                \[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)} \]
                              2. Add Preprocessing
                              3. Taylor expanded in phi1 around 0

                                \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_2\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_2}^{2}}} \]
                              4. Step-by-step derivation
                                1. Applied rewrites82.2%

                                  \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_2\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_2\right)} \]
                                2. Taylor expanded in lambda1 around inf

                                  \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 \cdot \cos \left(\frac{-1}{2} \cdot \phi_2\right), \phi_2\right) \]
                                3. Step-by-step derivation
                                  1. Applied rewrites79.1%

                                    \[\leadsto R \cdot \mathsf{hypot}\left(\cos \left(0.5 \cdot \phi_2\right) \cdot \lambda_1, \phi_2\right) \]
                                4. Recombined 2 regimes into one program.
                                5. Add Preprocessing

                                Alternative 5: 81.1% accurate, 1.3× speedup?

                                \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \begin{array}{l} \mathbf{if}\;\phi_1 \leq -215000:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(\cos \left(0.5 \cdot \phi_1\right) \cdot \lambda_2, \phi_1\right)\\ \mathbf{else}:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \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 -215000.0)
                                   (* R (hypot (* (cos (* 0.5 phi1)) lambda2) phi1))
                                   (* R (hypot (- lambda1 lambda2) 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 <= -215000.0) {
                                		tmp = R * hypot((cos((0.5 * phi1)) * lambda2), phi1);
                                	} else {
                                		tmp = R * hypot((lambda1 - lambda2), 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 <= -215000.0) {
                                		tmp = R * Math.hypot((Math.cos((0.5 * phi1)) * lambda2), phi1);
                                	} else {
                                		tmp = R * Math.hypot((lambda1 - lambda2), 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 <= -215000.0:
                                		tmp = R * math.hypot((math.cos((0.5 * phi1)) * lambda2), phi1)
                                	else:
                                		tmp = R * math.hypot((lambda1 - lambda2), 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 <= -215000.0)
                                		tmp = Float64(R * hypot(Float64(cos(Float64(0.5 * phi1)) * lambda2), phi1));
                                	else
                                		tmp = Float64(R * hypot(Float64(lambda1 - lambda2), 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 <= -215000.0)
                                		tmp = R * hypot((cos((0.5 * phi1)) * lambda2), phi1);
                                	else
                                		tmp = R * hypot((lambda1 - lambda2), 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, -215000.0], N[(R * N[Sqrt[N[(N[Cos[N[(0.5 * phi1), $MachinePrecision]], $MachinePrecision] * lambda2), $MachinePrecision] ^ 2 + phi1 ^ 2], $MachinePrecision]), $MachinePrecision], N[(R * N[Sqrt[N[(lambda1 - lambda2), $MachinePrecision] ^ 2 + phi2 ^ 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 -215000:\\
                                \;\;\;\;R \cdot \mathsf{hypot}\left(\cos \left(0.5 \cdot \phi_1\right) \cdot \lambda_2, \phi_1\right)\\
                                
                                \mathbf{else}:\\
                                \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_2\right)\\
                                
                                
                                \end{array}
                                \end{array}
                                
                                Derivation
                                1. Split input into 2 regimes
                                2. if phi1 < -215000

                                  1. Initial program 53.5%

                                    \[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)} \]
                                  2. Add Preprocessing
                                  3. Taylor expanded in phi2 around 0

                                    \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_1\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_1}^{2}}} \]
                                  4. Step-by-step derivation
                                    1. Applied rewrites83.4%

                                      \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_1\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_1\right)} \]
                                    2. Taylor expanded in lambda1 around 0

                                      \[\leadsto R \cdot \sqrt{{\lambda_2}^{2} \cdot {\cos \left(\frac{-1}{2} \cdot \phi_1\right)}^{2} + {\phi_1}^{2}} \]
                                    3. Step-by-step derivation
                                      1. Applied rewrites76.2%

                                        \[\leadsto R \cdot \mathsf{hypot}\left(\cos \left(0.5 \cdot \phi_1\right) \cdot \lambda_2, \color{blue}{\phi_1}\right) \]

                                      if -215000 < phi1

                                      1. Initial program 65.9%

                                        \[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)} \]
                                      2. Add Preprocessing
                                      3. Taylor expanded in phi1 around 0

                                        \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_2\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_2}^{2}}} \]
                                      4. Step-by-step derivation
                                        1. Applied rewrites82.6%

                                          \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_2\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_2\right)} \]
                                        2. Taylor expanded in phi2 around 0

                                          \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_2\right) \]
                                        3. Step-by-step derivation
                                          1. Applied rewrites77.1%

                                            \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_2\right) \]
                                        4. Recombined 2 regimes into one program.
                                        5. Add Preprocessing

                                        Alternative 6: 80.2% accurate, 2.4× speedup?

                                        \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \begin{array}{l} \mathbf{if}\;\phi_2 \leq 5.2 \cdot 10^{+78}:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_1\right)\\ \mathbf{else}:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \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 5.2e+78)
                                           (* R (hypot (- lambda1 lambda2) phi1))
                                           (* R (hypot (- lambda1 lambda2) 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 <= 5.2e+78) {
                                        		tmp = R * hypot((lambda1 - lambda2), phi1);
                                        	} else {
                                        		tmp = R * hypot((lambda1 - lambda2), 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 <= 5.2e+78) {
                                        		tmp = R * Math.hypot((lambda1 - lambda2), phi1);
                                        	} else {
                                        		tmp = R * Math.hypot((lambda1 - lambda2), 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 <= 5.2e+78:
                                        		tmp = R * math.hypot((lambda1 - lambda2), phi1)
                                        	else:
                                        		tmp = R * math.hypot((lambda1 - lambda2), 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 <= 5.2e+78)
                                        		tmp = Float64(R * hypot(Float64(lambda1 - lambda2), phi1));
                                        	else
                                        		tmp = Float64(R * hypot(Float64(lambda1 - lambda2), 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 <= 5.2e+78)
                                        		tmp = R * hypot((lambda1 - lambda2), phi1);
                                        	else
                                        		tmp = R * hypot((lambda1 - lambda2), 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, 5.2e+78], N[(R * N[Sqrt[N[(lambda1 - lambda2), $MachinePrecision] ^ 2 + phi1 ^ 2], $MachinePrecision]), $MachinePrecision], N[(R * N[Sqrt[N[(lambda1 - lambda2), $MachinePrecision] ^ 2 + phi2 ^ 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 5.2 \cdot 10^{+78}:\\
                                        \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_1\right)\\
                                        
                                        \mathbf{else}:\\
                                        \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_2\right)\\
                                        
                                        
                                        \end{array}
                                        \end{array}
                                        
                                        Derivation
                                        1. Split input into 2 regimes
                                        2. if phi2 < 5.2e78

                                          1. Initial program 63.8%

                                            \[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)} \]
                                          2. Add Preprocessing
                                          3. Taylor expanded in phi2 around 0

                                            \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_1\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_1}^{2}}} \]
                                          4. Step-by-step derivation
                                            1. Applied rewrites83.9%

                                              \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_1\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_1\right)} \]
                                            2. Taylor expanded in phi1 around 0

                                              \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_1\right) \]
                                            3. Step-by-step derivation
                                              1. Applied rewrites78.1%

                                                \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_1\right) \]

                                              if 5.2e78 < phi2

                                              1. Initial program 57.7%

                                                \[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)} \]
                                              2. Add Preprocessing
                                              3. Taylor expanded in phi1 around 0

                                                \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_2\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_2}^{2}}} \]
                                              4. Step-by-step derivation
                                                1. Applied rewrites83.7%

                                                  \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_2\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_2\right)} \]
                                                2. Taylor expanded in phi2 around 0

                                                  \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_2\right) \]
                                                3. Step-by-step derivation
                                                  1. Applied rewrites82.0%

                                                    \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_2\right) \]
                                                4. Recombined 2 regimes into one program.
                                                5. Add Preprocessing

                                                Alternative 7: 79.8% accurate, 2.4× speedup?

                                                \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \begin{array}{l} \mathbf{if}\;\phi_2 \leq 2.2 \cdot 10^{+83}:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_1\right)\\ \mathbf{else}:\\ \;\;\;\;\left(\left(-\phi_1\right) + \phi_2\right) \cdot R\\ \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 2.2e+83)
                                                   (* R (hypot (- lambda1 lambda2) phi1))
                                                   (* (+ (- phi1) phi2) R)))
                                                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.2e+83) {
                                                		tmp = R * hypot((lambda1 - lambda2), phi1);
                                                	} else {
                                                		tmp = (-phi1 + phi2) * R;
                                                	}
                                                	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 <= 2.2e+83) {
                                                		tmp = R * Math.hypot((lambda1 - lambda2), phi1);
                                                	} else {
                                                		tmp = (-phi1 + phi2) * R;
                                                	}
                                                	return tmp;
                                                }
                                                
                                                [R, lambda1, lambda2, phi1, phi2] = sort([R, lambda1, lambda2, phi1, phi2])
                                                def code(R, lambda1, lambda2, phi1, phi2):
                                                	tmp = 0
                                                	if phi2 <= 2.2e+83:
                                                		tmp = R * math.hypot((lambda1 - lambda2), phi1)
                                                	else:
                                                		tmp = (-phi1 + phi2) * R
                                                	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.2e+83)
                                                		tmp = Float64(R * hypot(Float64(lambda1 - lambda2), phi1));
                                                	else
                                                		tmp = Float64(Float64(Float64(-phi1) + phi2) * R);
                                                	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 <= 2.2e+83)
                                                		tmp = R * hypot((lambda1 - lambda2), phi1);
                                                	else
                                                		tmp = (-phi1 + phi2) * R;
                                                	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, 2.2e+83], N[(R * N[Sqrt[N[(lambda1 - lambda2), $MachinePrecision] ^ 2 + phi1 ^ 2], $MachinePrecision]), $MachinePrecision], N[(N[((-phi1) + phi2), $MachinePrecision] * R), $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.2 \cdot 10^{+83}:\\
                                                \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_1\right)\\
                                                
                                                \mathbf{else}:\\
                                                \;\;\;\;\left(\left(-\phi_1\right) + \phi_2\right) \cdot R\\
                                                
                                                
                                                \end{array}
                                                \end{array}
                                                
                                                Derivation
                                                1. Split input into 2 regimes
                                                2. if phi2 < 2.19999999999999999e83

                                                  1. Initial program 63.6%

                                                    \[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)} \]
                                                  2. Add Preprocessing
                                                  3. Taylor expanded in phi2 around 0

                                                    \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_1\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_1}^{2}}} \]
                                                  4. Step-by-step derivation
                                                    1. Applied rewrites83.6%

                                                      \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_1\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_1\right)} \]
                                                    2. Taylor expanded in phi1 around 0

                                                      \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_1\right) \]
                                                    3. Step-by-step derivation
                                                      1. Applied rewrites77.8%

                                                        \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 - \lambda_2, \phi_1\right) \]

                                                      if 2.19999999999999999e83 < phi2

                                                      1. Initial program 59.0%

                                                        \[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)} \]
                                                      2. Add Preprocessing
                                                      3. Taylor expanded in phi1 around -inf

                                                        \[\leadsto \color{blue}{-1 \cdot \left(\phi_1 \cdot \left(R + -1 \cdot \frac{R \cdot \phi_2}{\phi_1}\right)\right)} \]
                                                      4. Step-by-step derivation
                                                        1. Applied rewrites67.1%

                                                          \[\leadsto \color{blue}{\left(-\phi_1\right) \cdot \mathsf{fma}\left(-R, \frac{\phi_2}{\phi_1}, R\right)} \]
                                                        2. Taylor expanded in phi1 around 0

                                                          \[\leadsto \left(-\phi_1\right) \cdot \frac{-1 \cdot \left(R \cdot \phi_2\right) + R \cdot \phi_1}{\color{blue}{\phi_1}} \]
                                                        3. Step-by-step derivation
                                                          1. Applied rewrites80.4%

                                                            \[\leadsto \left(-\phi_1\right) \cdot \frac{\mathsf{fma}\left(-R, \phi_2, \phi_1 \cdot R\right)}{\color{blue}{\phi_1}} \]
                                                          2. Taylor expanded in phi1 around 0

                                                            \[\leadsto -1 \cdot \left(R \cdot \phi_1\right) + \color{blue}{R \cdot \phi_2} \]
                                                          3. Step-by-step derivation
                                                            1. Applied rewrites82.7%

                                                              \[\leadsto \mathsf{fma}\left(-1, \phi_1, \phi_2\right) \cdot \color{blue}{R} \]
                                                            2. Step-by-step derivation
                                                              1. Applied rewrites82.7%

                                                                \[\leadsto \left(\left(-\phi_1\right) + \phi_2\right) \cdot R \]
                                                            3. Recombined 2 regimes into one program.
                                                            4. Add Preprocessing

                                                            Alternative 8: 68.6% accurate, 2.5× speedup?

                                                            \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \begin{array}{l} \mathbf{if}\;\phi_1 \leq -2.25 \cdot 10^{+51}:\\ \;\;\;\;\left(\left(-\phi_1\right) + \phi_2\right) \cdot R\\ \mathbf{else}:\\ \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_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 -2.25e+51) (* (+ (- phi1) phi2) R) (* R (hypot lambda1 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 <= -2.25e+51) {
                                                            		tmp = (-phi1 + phi2) * R;
                                                            	} else {
                                                            		tmp = R * hypot(lambda1, 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 <= -2.25e+51) {
                                                            		tmp = (-phi1 + phi2) * R;
                                                            	} else {
                                                            		tmp = R * Math.hypot(lambda1, 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 <= -2.25e+51:
                                                            		tmp = (-phi1 + phi2) * R
                                                            	else:
                                                            		tmp = R * math.hypot(lambda1, 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 <= -2.25e+51)
                                                            		tmp = Float64(Float64(Float64(-phi1) + phi2) * R);
                                                            	else
                                                            		tmp = Float64(R * hypot(lambda1, 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 <= -2.25e+51)
                                                            		tmp = (-phi1 + phi2) * R;
                                                            	else
                                                            		tmp = R * hypot(lambda1, 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, -2.25e+51], N[(N[((-phi1) + phi2), $MachinePrecision] * R), $MachinePrecision], N[(R * N[Sqrt[lambda1 ^ 2 + phi2 ^ 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 -2.25 \cdot 10^{+51}:\\
                                                            \;\;\;\;\left(\left(-\phi_1\right) + \phi_2\right) \cdot R\\
                                                            
                                                            \mathbf{else}:\\
                                                            \;\;\;\;R \cdot \mathsf{hypot}\left(\lambda_1, \phi_2\right)\\
                                                            
                                                            
                                                            \end{array}
                                                            \end{array}
                                                            
                                                            Derivation
                                                            1. Split input into 2 regimes
                                                            2. if phi1 < -2.25e51

                                                              1. Initial program 54.9%

                                                                \[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)} \]
                                                              2. Add Preprocessing
                                                              3. Taylor expanded in phi1 around -inf

                                                                \[\leadsto \color{blue}{-1 \cdot \left(\phi_1 \cdot \left(R + -1 \cdot \frac{R \cdot \phi_2}{\phi_1}\right)\right)} \]
                                                              4. Step-by-step derivation
                                                                1. Applied rewrites72.8%

                                                                  \[\leadsto \color{blue}{\left(-\phi_1\right) \cdot \mathsf{fma}\left(-R, \frac{\phi_2}{\phi_1}, R\right)} \]
                                                                2. Taylor expanded in phi1 around 0

                                                                  \[\leadsto \left(-\phi_1\right) \cdot \frac{-1 \cdot \left(R \cdot \phi_2\right) + R \cdot \phi_1}{\color{blue}{\phi_1}} \]
                                                                3. Step-by-step derivation
                                                                  1. Applied rewrites76.1%

                                                                    \[\leadsto \left(-\phi_1\right) \cdot \frac{\mathsf{fma}\left(-R, \phi_2, \phi_1 \cdot R\right)}{\color{blue}{\phi_1}} \]
                                                                  2. Taylor expanded in phi1 around 0

                                                                    \[\leadsto -1 \cdot \left(R \cdot \phi_1\right) + \color{blue}{R \cdot \phi_2} \]
                                                                  3. Step-by-step derivation
                                                                    1. Applied rewrites72.8%

                                                                      \[\leadsto \mathsf{fma}\left(-1, \phi_1, \phi_2\right) \cdot \color{blue}{R} \]
                                                                    2. Step-by-step derivation
                                                                      1. Applied rewrites72.8%

                                                                        \[\leadsto \left(\left(-\phi_1\right) + \phi_2\right) \cdot R \]

                                                                      if -2.25e51 < phi1

                                                                      1. Initial program 65.1%

                                                                        \[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)} \]
                                                                      2. Add Preprocessing
                                                                      3. Taylor expanded in phi1 around 0

                                                                        \[\leadsto R \cdot \color{blue}{\sqrt{{\cos \left(\frac{1}{2} \cdot \phi_2\right)}^{2} \cdot {\left(\lambda_1 - \lambda_2\right)}^{2} + {\phi_2}^{2}}} \]
                                                                      4. Step-by-step derivation
                                                                        1. Applied rewrites81.9%

                                                                          \[\leadsto R \cdot \color{blue}{\mathsf{hypot}\left(\cos \left(-0.5 \cdot \phi_2\right) \cdot \left(\lambda_1 - \lambda_2\right), \phi_2\right)} \]
                                                                        2. Taylor expanded in lambda1 around inf

                                                                          \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1 \cdot \cos \left(\frac{-1}{2} \cdot \phi_2\right), \phi_2\right) \]
                                                                        3. Step-by-step derivation
                                                                          1. Applied rewrites57.8%

                                                                            \[\leadsto R \cdot \mathsf{hypot}\left(\cos \left(0.5 \cdot \phi_2\right) \cdot \lambda_1, \phi_2\right) \]
                                                                          2. Taylor expanded in phi2 around 0

                                                                            \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1, \phi_2\right) \]
                                                                          3. Step-by-step derivation
                                                                            1. Applied rewrites54.4%

                                                                              \[\leadsto R \cdot \mathsf{hypot}\left(\lambda_1, \phi_2\right) \]
                                                                          4. Recombined 2 regimes into one program.
                                                                          5. Add Preprocessing

                                                                          Alternative 9: 58.9% accurate, 5.9× speedup?

                                                                          \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \begin{array}{l} \mathbf{if}\;R \leq 4.25 \cdot 10^{+80}:\\ \;\;\;\;\left(\left(-\phi_1\right) + \phi_2\right) \cdot R\\ \mathbf{else}:\\ \;\;\;\;\left(-\phi_1\right) \cdot \left(\mathsf{fma}\left(\frac{R}{\phi_1}, -1, \frac{R}{\phi_2}\right) \cdot \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 (<= R 4.25e+80)
                                                                             (* (+ (- phi1) phi2) R)
                                                                             (* (- phi1) (* (fma (/ R phi1) -1.0 (/ R 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 (R <= 4.25e+80) {
                                                                          		tmp = (-phi1 + phi2) * R;
                                                                          	} else {
                                                                          		tmp = -phi1 * (fma((R / phi1), -1.0, (R / 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 (R <= 4.25e+80)
                                                                          		tmp = Float64(Float64(Float64(-phi1) + phi2) * R);
                                                                          	else
                                                                          		tmp = Float64(Float64(-phi1) * Float64(fma(Float64(R / phi1), -1.0, Float64(R / 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[R, 4.25e+80], N[(N[((-phi1) + phi2), $MachinePrecision] * R), $MachinePrecision], N[((-phi1) * N[(N[(N[(R / phi1), $MachinePrecision] * -1.0 + N[(R / phi2), $MachinePrecision]), $MachinePrecision] * phi2), $MachinePrecision]), $MachinePrecision]]
                                                                          
                                                                          \begin{array}{l}
                                                                          [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\
                                                                          \\
                                                                          \begin{array}{l}
                                                                          \mathbf{if}\;R \leq 4.25 \cdot 10^{+80}:\\
                                                                          \;\;\;\;\left(\left(-\phi_1\right) + \phi_2\right) \cdot R\\
                                                                          
                                                                          \mathbf{else}:\\
                                                                          \;\;\;\;\left(-\phi_1\right) \cdot \left(\mathsf{fma}\left(\frac{R}{\phi_1}, -1, \frac{R}{\phi_2}\right) \cdot \phi_2\right)\\
                                                                          
                                                                          
                                                                          \end{array}
                                                                          \end{array}
                                                                          
                                                                          Derivation
                                                                          1. Split input into 2 regimes
                                                                          2. if R < 4.25000000000000003e80

                                                                            1. Initial program 52.6%

                                                                              \[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)} \]
                                                                            2. Add Preprocessing
                                                                            3. Taylor expanded in phi1 around -inf

                                                                              \[\leadsto \color{blue}{-1 \cdot \left(\phi_1 \cdot \left(R + -1 \cdot \frac{R \cdot \phi_2}{\phi_1}\right)\right)} \]
                                                                            4. Step-by-step derivation
                                                                              1. Applied rewrites26.0%

                                                                                \[\leadsto \color{blue}{\left(-\phi_1\right) \cdot \mathsf{fma}\left(-R, \frac{\phi_2}{\phi_1}, R\right)} \]
                                                                              2. Taylor expanded in phi1 around 0

                                                                                \[\leadsto \left(-\phi_1\right) \cdot \frac{-1 \cdot \left(R \cdot \phi_2\right) + R \cdot \phi_1}{\color{blue}{\phi_1}} \]
                                                                              3. Step-by-step derivation
                                                                                1. Applied rewrites29.5%

                                                                                  \[\leadsto \left(-\phi_1\right) \cdot \frac{\mathsf{fma}\left(-R, \phi_2, \phi_1 \cdot R\right)}{\color{blue}{\phi_1}} \]
                                                                                2. Taylor expanded in phi1 around 0

                                                                                  \[\leadsto -1 \cdot \left(R \cdot \phi_1\right) + \color{blue}{R \cdot \phi_2} \]
                                                                                3. Step-by-step derivation
                                                                                  1. Applied rewrites28.1%

                                                                                    \[\leadsto \mathsf{fma}\left(-1, \phi_1, \phi_2\right) \cdot \color{blue}{R} \]
                                                                                  2. Step-by-step derivation
                                                                                    1. Applied rewrites28.1%

                                                                                      \[\leadsto \left(\left(-\phi_1\right) + \phi_2\right) \cdot R \]

                                                                                    if 4.25000000000000003e80 < R

                                                                                    1. Initial program 98.3%

                                                                                      \[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)} \]
                                                                                    2. Add Preprocessing
                                                                                    3. Taylor expanded in phi1 around -inf

                                                                                      \[\leadsto \color{blue}{-1 \cdot \left(\phi_1 \cdot \left(R + -1 \cdot \frac{R \cdot \phi_2}{\phi_1}\right)\right)} \]
                                                                                    4. Step-by-step derivation
                                                                                      1. Applied rewrites40.1%

                                                                                        \[\leadsto \color{blue}{\left(-\phi_1\right) \cdot \mathsf{fma}\left(-R, \frac{\phi_2}{\phi_1}, R\right)} \]
                                                                                      2. Taylor expanded in phi2 around inf

                                                                                        \[\leadsto \left(-\phi_1\right) \cdot \left(\phi_2 \cdot \color{blue}{\left(-1 \cdot \frac{R}{\phi_1} + \frac{R}{\phi_2}\right)}\right) \]
                                                                                      3. Step-by-step derivation
                                                                                        1. Applied rewrites46.5%

                                                                                          \[\leadsto \left(-\phi_1\right) \cdot \left(\mathsf{fma}\left(\frac{R}{\phi_1}, -1, \frac{R}{\phi_2}\right) \cdot \color{blue}{\phi_2}\right) \]
                                                                                      4. Recombined 2 regimes into one program.
                                                                                      5. Add Preprocessing

                                                                                      Alternative 10: 59.9% accurate, 8.5× speedup?

                                                                                      \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \begin{array}{l} \mathbf{if}\;\phi_2 \leq 2 \cdot 10^{+83}:\\ \;\;\;\;\left(-\phi_1\right) \cdot \mathsf{fma}\left(\phi_2, \frac{-R}{\phi_1}, R\right)\\ \mathbf{else}:\\ \;\;\;\;\left(\left(-\phi_1\right) + \phi_2\right) \cdot R\\ \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 2e+83)
                                                                                         (* (- phi1) (fma phi2 (/ (- R) phi1) R))
                                                                                         (* (+ (- phi1) phi2) R)))
                                                                                      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 <= 2e+83) {
                                                                                      		tmp = -phi1 * fma(phi2, (-R / phi1), R);
                                                                                      	} else {
                                                                                      		tmp = (-phi1 + phi2) * R;
                                                                                      	}
                                                                                      	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 <= 2e+83)
                                                                                      		tmp = Float64(Float64(-phi1) * fma(phi2, Float64(Float64(-R) / phi1), R));
                                                                                      	else
                                                                                      		tmp = Float64(Float64(Float64(-phi1) + phi2) * R);
                                                                                      	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, 2e+83], N[((-phi1) * N[(phi2 * N[((-R) / phi1), $MachinePrecision] + R), $MachinePrecision]), $MachinePrecision], N[(N[((-phi1) + phi2), $MachinePrecision] * R), $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 \cdot 10^{+83}:\\
                                                                                      \;\;\;\;\left(-\phi_1\right) \cdot \mathsf{fma}\left(\phi_2, \frac{-R}{\phi_1}, R\right)\\
                                                                                      
                                                                                      \mathbf{else}:\\
                                                                                      \;\;\;\;\left(\left(-\phi_1\right) + \phi_2\right) \cdot R\\
                                                                                      
                                                                                      
                                                                                      \end{array}
                                                                                      \end{array}
                                                                                      
                                                                                      Derivation
                                                                                      1. Split input into 2 regimes
                                                                                      2. if phi2 < 2.00000000000000006e83

                                                                                        1. Initial program 63.6%

                                                                                          \[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)} \]
                                                                                        2. Add Preprocessing
                                                                                        3. Taylor expanded in phi1 around -inf

                                                                                          \[\leadsto \color{blue}{-1 \cdot \left(\phi_1 \cdot \left(R + -1 \cdot \frac{R \cdot \phi_2}{\phi_1}\right)\right)} \]
                                                                                        4. Step-by-step derivation
                                                                                          1. Applied rewrites21.3%

                                                                                            \[\leadsto \color{blue}{\left(-\phi_1\right) \cdot \mathsf{fma}\left(-R, \frac{\phi_2}{\phi_1}, R\right)} \]
                                                                                          2. Taylor expanded in phi1 around 0

                                                                                            \[\leadsto \left(-\phi_1\right) \cdot \frac{-1 \cdot \left(R \cdot \phi_2\right) + R \cdot \phi_1}{\color{blue}{\phi_1}} \]
                                                                                          3. Step-by-step derivation
                                                                                            1. Applied rewrites22.3%

                                                                                              \[\leadsto \left(-\phi_1\right) \cdot \frac{\mathsf{fma}\left(-R, \phi_2, \phi_1 \cdot R\right)}{\color{blue}{\phi_1}} \]
                                                                                            2. Step-by-step derivation
                                                                                              1. Applied rewrites26.3%

                                                                                                \[\leadsto \left(-\phi_1\right) \cdot \mathsf{fma}\left(\phi_2, \frac{-R}{\color{blue}{\phi_1}}, \frac{\phi_1 \cdot R}{\phi_1}\right) \]
                                                                                              2. Taylor expanded in R around 0

                                                                                                \[\leadsto \left(-\phi_1\right) \cdot \mathsf{fma}\left(\phi_2, \frac{-R}{\phi_1}, R\right) \]
                                                                                              3. Step-by-step derivation
                                                                                                1. Applied rewrites25.4%

                                                                                                  \[\leadsto \left(-\phi_1\right) \cdot \mathsf{fma}\left(\phi_2, \frac{-R}{\phi_1}, R\right) \]

                                                                                                if 2.00000000000000006e83 < phi2

                                                                                                1. Initial program 59.0%

                                                                                                  \[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)} \]
                                                                                                2. Add Preprocessing
                                                                                                3. Taylor expanded in phi1 around -inf

                                                                                                  \[\leadsto \color{blue}{-1 \cdot \left(\phi_1 \cdot \left(R + -1 \cdot \frac{R \cdot \phi_2}{\phi_1}\right)\right)} \]
                                                                                                4. Step-by-step derivation
                                                                                                  1. Applied rewrites67.1%

                                                                                                    \[\leadsto \color{blue}{\left(-\phi_1\right) \cdot \mathsf{fma}\left(-R, \frac{\phi_2}{\phi_1}, R\right)} \]
                                                                                                  2. Taylor expanded in phi1 around 0

                                                                                                    \[\leadsto \left(-\phi_1\right) \cdot \frac{-1 \cdot \left(R \cdot \phi_2\right) + R \cdot \phi_1}{\color{blue}{\phi_1}} \]
                                                                                                  3. Step-by-step derivation
                                                                                                    1. Applied rewrites80.4%

                                                                                                      \[\leadsto \left(-\phi_1\right) \cdot \frac{\mathsf{fma}\left(-R, \phi_2, \phi_1 \cdot R\right)}{\color{blue}{\phi_1}} \]
                                                                                                    2. Taylor expanded in phi1 around 0

                                                                                                      \[\leadsto -1 \cdot \left(R \cdot \phi_1\right) + \color{blue}{R \cdot \phi_2} \]
                                                                                                    3. Step-by-step derivation
                                                                                                      1. Applied rewrites82.7%

                                                                                                        \[\leadsto \mathsf{fma}\left(-1, \phi_1, \phi_2\right) \cdot \color{blue}{R} \]
                                                                                                      2. Step-by-step derivation
                                                                                                        1. Applied rewrites82.7%

                                                                                                          \[\leadsto \left(\left(-\phi_1\right) + \phi_2\right) \cdot R \]
                                                                                                      3. Recombined 2 regimes into one program.
                                                                                                      4. Add Preprocessing

                                                                                                      Alternative 11: 58.4% accurate, 9.0× speedup?

                                                                                                      \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \begin{array}{l} \mathbf{if}\;R \leq 2.4 \cdot 10^{+104}:\\ \;\;\;\;\left(\left(-\phi_1\right) + \phi_2\right) \cdot R\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(-R, \frac{\phi_1}{\phi_2}, R\right) \cdot \phi_2\\ \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 (<= R 2.4e+104)
                                                                                                         (* (+ (- phi1) phi2) R)
                                                                                                         (* (fma (- R) (/ phi1 phi2) R) 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 (R <= 2.4e+104) {
                                                                                                      		tmp = (-phi1 + phi2) * R;
                                                                                                      	} else {
                                                                                                      		tmp = fma(-R, (phi1 / phi2), R) * 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 (R <= 2.4e+104)
                                                                                                      		tmp = Float64(Float64(Float64(-phi1) + phi2) * R);
                                                                                                      	else
                                                                                                      		tmp = Float64(fma(Float64(-R), Float64(phi1 / phi2), R) * 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[R, 2.4e+104], N[(N[((-phi1) + phi2), $MachinePrecision] * R), $MachinePrecision], N[(N[((-R) * N[(phi1 / phi2), $MachinePrecision] + R), $MachinePrecision] * phi2), $MachinePrecision]]
                                                                                                      
                                                                                                      \begin{array}{l}
                                                                                                      [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\
                                                                                                      \\
                                                                                                      \begin{array}{l}
                                                                                                      \mathbf{if}\;R \leq 2.4 \cdot 10^{+104}:\\
                                                                                                      \;\;\;\;\left(\left(-\phi_1\right) + \phi_2\right) \cdot R\\
                                                                                                      
                                                                                                      \mathbf{else}:\\
                                                                                                      \;\;\;\;\mathsf{fma}\left(-R, \frac{\phi_1}{\phi_2}, R\right) \cdot \phi_2\\
                                                                                                      
                                                                                                      
                                                                                                      \end{array}
                                                                                                      \end{array}
                                                                                                      
                                                                                                      Derivation
                                                                                                      1. Split input into 2 regimes
                                                                                                      2. if R < 2.4e104

                                                                                                        1. Initial program 53.5%

                                                                                                          \[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)} \]
                                                                                                        2. Add Preprocessing
                                                                                                        3. Taylor expanded in phi1 around -inf

                                                                                                          \[\leadsto \color{blue}{-1 \cdot \left(\phi_1 \cdot \left(R + -1 \cdot \frac{R \cdot \phi_2}{\phi_1}\right)\right)} \]
                                                                                                        4. Step-by-step derivation
                                                                                                          1. Applied rewrites26.0%

                                                                                                            \[\leadsto \color{blue}{\left(-\phi_1\right) \cdot \mathsf{fma}\left(-R, \frac{\phi_2}{\phi_1}, R\right)} \]
                                                                                                          2. Taylor expanded in phi1 around 0

                                                                                                            \[\leadsto \left(-\phi_1\right) \cdot \frac{-1 \cdot \left(R \cdot \phi_2\right) + R \cdot \phi_1}{\color{blue}{\phi_1}} \]
                                                                                                          3. Step-by-step derivation
                                                                                                            1. Applied rewrites29.4%

                                                                                                              \[\leadsto \left(-\phi_1\right) \cdot \frac{\mathsf{fma}\left(-R, \phi_2, \phi_1 \cdot R\right)}{\color{blue}{\phi_1}} \]
                                                                                                            2. Taylor expanded in phi1 around 0

                                                                                                              \[\leadsto -1 \cdot \left(R \cdot \phi_1\right) + \color{blue}{R \cdot \phi_2} \]
                                                                                                            3. Step-by-step derivation
                                                                                                              1. Applied rewrites28.1%

                                                                                                                \[\leadsto \mathsf{fma}\left(-1, \phi_1, \phi_2\right) \cdot \color{blue}{R} \]
                                                                                                              2. Step-by-step derivation
                                                                                                                1. Applied rewrites28.1%

                                                                                                                  \[\leadsto \left(\left(-\phi_1\right) + \phi_2\right) \cdot R \]

                                                                                                                if 2.4e104 < R

                                                                                                                1. Initial program 98.2%

                                                                                                                  \[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)} \]
                                                                                                                2. Add Preprocessing
                                                                                                                3. Taylor expanded in phi2 around inf

                                                                                                                  \[\leadsto \color{blue}{\phi_2 \cdot \left(R + -1 \cdot \frac{R \cdot \phi_1}{\phi_2}\right)} \]
                                                                                                                4. Step-by-step derivation
                                                                                                                  1. Applied rewrites39.3%

                                                                                                                    \[\leadsto \color{blue}{\mathsf{fma}\left(-R, \frac{\phi_1}{\phi_2}, R\right) \cdot \phi_2} \]
                                                                                                                5. Recombined 2 regimes into one program.
                                                                                                                6. Add Preprocessing

                                                                                                                Alternative 12: 58.2% accurate, 9.0× speedup?

                                                                                                                \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \begin{array}{l} \mathbf{if}\;R \leq 2 \cdot 10^{+98}:\\ \;\;\;\;\left(\left(-\phi_1\right) + \phi_2\right) \cdot R\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(R, \frac{\phi_2}{\phi_1}, -R\right) \cdot \phi_1\\ \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 (<= R 2e+98)
                                                                                                                   (* (+ (- phi1) phi2) R)
                                                                                                                   (* (fma R (/ phi2 phi1) (- R)) phi1)))
                                                                                                                assert(R < lambda1 && lambda1 < lambda2 && lambda2 < phi1 && phi1 < phi2);
                                                                                                                double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
                                                                                                                	double tmp;
                                                                                                                	if (R <= 2e+98) {
                                                                                                                		tmp = (-phi1 + phi2) * R;
                                                                                                                	} else {
                                                                                                                		tmp = fma(R, (phi2 / phi1), -R) * phi1;
                                                                                                                	}
                                                                                                                	return tmp;
                                                                                                                }
                                                                                                                
                                                                                                                R, lambda1, lambda2, phi1, phi2 = sort([R, lambda1, lambda2, phi1, phi2])
                                                                                                                function code(R, lambda1, lambda2, phi1, phi2)
                                                                                                                	tmp = 0.0
                                                                                                                	if (R <= 2e+98)
                                                                                                                		tmp = Float64(Float64(Float64(-phi1) + phi2) * R);
                                                                                                                	else
                                                                                                                		tmp = Float64(fma(R, Float64(phi2 / phi1), Float64(-R)) * phi1);
                                                                                                                	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[R, 2e+98], N[(N[((-phi1) + phi2), $MachinePrecision] * R), $MachinePrecision], N[(N[(R * N[(phi2 / phi1), $MachinePrecision] + (-R)), $MachinePrecision] * phi1), $MachinePrecision]]
                                                                                                                
                                                                                                                \begin{array}{l}
                                                                                                                [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\
                                                                                                                \\
                                                                                                                \begin{array}{l}
                                                                                                                \mathbf{if}\;R \leq 2 \cdot 10^{+98}:\\
                                                                                                                \;\;\;\;\left(\left(-\phi_1\right) + \phi_2\right) \cdot R\\
                                                                                                                
                                                                                                                \mathbf{else}:\\
                                                                                                                \;\;\;\;\mathsf{fma}\left(R, \frac{\phi_2}{\phi_1}, -R\right) \cdot \phi_1\\
                                                                                                                
                                                                                                                
                                                                                                                \end{array}
                                                                                                                \end{array}
                                                                                                                
                                                                                                                Derivation
                                                                                                                1. Split input into 2 regimes
                                                                                                                2. if R < 2e98

                                                                                                                  1. Initial program 53.5%

                                                                                                                    \[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)} \]
                                                                                                                  2. Add Preprocessing
                                                                                                                  3. Taylor expanded in phi1 around -inf

                                                                                                                    \[\leadsto \color{blue}{-1 \cdot \left(\phi_1 \cdot \left(R + -1 \cdot \frac{R \cdot \phi_2}{\phi_1}\right)\right)} \]
                                                                                                                  4. Step-by-step derivation
                                                                                                                    1. Applied rewrites26.0%

                                                                                                                      \[\leadsto \color{blue}{\left(-\phi_1\right) \cdot \mathsf{fma}\left(-R, \frac{\phi_2}{\phi_1}, R\right)} \]
                                                                                                                    2. Taylor expanded in phi1 around 0

                                                                                                                      \[\leadsto \left(-\phi_1\right) \cdot \frac{-1 \cdot \left(R \cdot \phi_2\right) + R \cdot \phi_1}{\color{blue}{\phi_1}} \]
                                                                                                                    3. Step-by-step derivation
                                                                                                                      1. Applied rewrites29.4%

                                                                                                                        \[\leadsto \left(-\phi_1\right) \cdot \frac{\mathsf{fma}\left(-R, \phi_2, \phi_1 \cdot R\right)}{\color{blue}{\phi_1}} \]
                                                                                                                      2. Taylor expanded in phi1 around 0

                                                                                                                        \[\leadsto -1 \cdot \left(R \cdot \phi_1\right) + \color{blue}{R \cdot \phi_2} \]
                                                                                                                      3. Step-by-step derivation
                                                                                                                        1. Applied rewrites28.1%

                                                                                                                          \[\leadsto \mathsf{fma}\left(-1, \phi_1, \phi_2\right) \cdot \color{blue}{R} \]
                                                                                                                        2. Step-by-step derivation
                                                                                                                          1. Applied rewrites28.1%

                                                                                                                            \[\leadsto \left(\left(-\phi_1\right) + \phi_2\right) \cdot R \]

                                                                                                                          if 2e98 < R

                                                                                                                          1. Initial program 98.2%

                                                                                                                            \[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)} \]
                                                                                                                          2. Add Preprocessing
                                                                                                                          3. Taylor expanded in phi1 around -inf

                                                                                                                            \[\leadsto \color{blue}{-1 \cdot \left(\phi_1 \cdot \left(R + -1 \cdot \frac{R \cdot \phi_2}{\phi_1}\right)\right)} \]
                                                                                                                          4. Step-by-step derivation
                                                                                                                            1. Applied rewrites41.1%

                                                                                                                              \[\leadsto \color{blue}{\left(-\phi_1\right) \cdot \mathsf{fma}\left(-R, \frac{\phi_2}{\phi_1}, R\right)} \]
                                                                                                                            2. Taylor expanded in phi1 around inf

                                                                                                                              \[\leadsto \phi_1 \cdot \color{blue}{\left(-1 \cdot R + \frac{R \cdot \phi_2}{\phi_1}\right)} \]
                                                                                                                            3. Step-by-step derivation
                                                                                                                              1. Applied rewrites41.1%

                                                                                                                                \[\leadsto \mathsf{fma}\left(R, \frac{\phi_2}{\phi_1}, -R\right) \cdot \color{blue}{\phi_1} \]
                                                                                                                            4. Recombined 2 regimes into one program.
                                                                                                                            5. Add Preprocessing

                                                                                                                            Alternative 13: 52.3% accurate, 19.9× speedup?

                                                                                                                            \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \begin{array}{l} \mathbf{if}\;\phi_2 \leq 8.8 \cdot 10^{+63}:\\ \;\;\;\;R \cdot \left(-\phi_1\right)\\ \mathbf{else}:\\ \;\;\;\;R \cdot \phi_2\\ \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 8.8e+63) (* R (- phi1)) (* R 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 <= 8.8e+63) {
                                                                                                                            		tmp = R * -phi1;
                                                                                                                            	} else {
                                                                                                                            		tmp = R * phi2;
                                                                                                                            	}
                                                                                                                            	return tmp;
                                                                                                                            }
                                                                                                                            
                                                                                                                            NOTE: R, lambda1, lambda2, phi1, and phi2 should be sorted in increasing order before calling this function.
                                                                                                                            module fmin_fmax_functions
                                                                                                                                implicit none
                                                                                                                                private
                                                                                                                                public fmax
                                                                                                                                public fmin
                                                                                                                            
                                                                                                                                interface fmax
                                                                                                                                    module procedure fmax88
                                                                                                                                    module procedure fmax44
                                                                                                                                    module procedure fmax84
                                                                                                                                    module procedure fmax48
                                                                                                                                end interface
                                                                                                                                interface fmin
                                                                                                                                    module procedure fmin88
                                                                                                                                    module procedure fmin44
                                                                                                                                    module procedure fmin84
                                                                                                                                    module procedure fmin48
                                                                                                                                end interface
                                                                                                                            contains
                                                                                                                                real(8) function fmax88(x, y) result (res)
                                                                                                                                    real(8), intent (in) :: x
                                                                                                                                    real(8), intent (in) :: y
                                                                                                                                    res = merge(y, merge(x, max(x, y), y /= y), x /= x)
                                                                                                                                end function
                                                                                                                                real(4) function fmax44(x, y) result (res)
                                                                                                                                    real(4), intent (in) :: x
                                                                                                                                    real(4), intent (in) :: y
                                                                                                                                    res = merge(y, merge(x, max(x, y), y /= y), x /= x)
                                                                                                                                end function
                                                                                                                                real(8) function fmax84(x, y) result(res)
                                                                                                                                    real(8), intent (in) :: x
                                                                                                                                    real(4), intent (in) :: y
                                                                                                                                    res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
                                                                                                                                end function
                                                                                                                                real(8) function fmax48(x, y) result(res)
                                                                                                                                    real(4), intent (in) :: x
                                                                                                                                    real(8), intent (in) :: y
                                                                                                                                    res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
                                                                                                                                end function
                                                                                                                                real(8) function fmin88(x, y) result (res)
                                                                                                                                    real(8), intent (in) :: x
                                                                                                                                    real(8), intent (in) :: y
                                                                                                                                    res = merge(y, merge(x, min(x, y), y /= y), x /= x)
                                                                                                                                end function
                                                                                                                                real(4) function fmin44(x, y) result (res)
                                                                                                                                    real(4), intent (in) :: x
                                                                                                                                    real(4), intent (in) :: y
                                                                                                                                    res = merge(y, merge(x, min(x, y), y /= y), x /= x)
                                                                                                                                end function
                                                                                                                                real(8) function fmin84(x, y) result(res)
                                                                                                                                    real(8), intent (in) :: x
                                                                                                                                    real(4), intent (in) :: y
                                                                                                                                    res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
                                                                                                                                end function
                                                                                                                                real(8) function fmin48(x, y) result(res)
                                                                                                                                    real(4), intent (in) :: x
                                                                                                                                    real(8), intent (in) :: y
                                                                                                                                    res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
                                                                                                                                end function
                                                                                                                            end module
                                                                                                                            
                                                                                                                            real(8) function code(r, lambda1, lambda2, phi1, phi2)
                                                                                                                            use fmin_fmax_functions
                                                                                                                                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) :: tmp
                                                                                                                                if (phi2 <= 8.8d+63) then
                                                                                                                                    tmp = r * -phi1
                                                                                                                                else
                                                                                                                                    tmp = r * phi2
                                                                                                                                end if
                                                                                                                                code = tmp
                                                                                                                            end function
                                                                                                                            
                                                                                                                            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 <= 8.8e+63) {
                                                                                                                            		tmp = R * -phi1;
                                                                                                                            	} else {
                                                                                                                            		tmp = R * 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 <= 8.8e+63:
                                                                                                                            		tmp = R * -phi1
                                                                                                                            	else:
                                                                                                                            		tmp = R * 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 <= 8.8e+63)
                                                                                                                            		tmp = Float64(R * Float64(-phi1));
                                                                                                                            	else
                                                                                                                            		tmp = Float64(R * 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 <= 8.8e+63)
                                                                                                                            		tmp = R * -phi1;
                                                                                                                            	else
                                                                                                                            		tmp = R * 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, 8.8e+63], N[(R * (-phi1)), $MachinePrecision], N[(R * phi2), $MachinePrecision]]
                                                                                                                            
                                                                                                                            \begin{array}{l}
                                                                                                                            [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\
                                                                                                                            \\
                                                                                                                            \begin{array}{l}
                                                                                                                            \mathbf{if}\;\phi_2 \leq 8.8 \cdot 10^{+63}:\\
                                                                                                                            \;\;\;\;R \cdot \left(-\phi_1\right)\\
                                                                                                                            
                                                                                                                            \mathbf{else}:\\
                                                                                                                            \;\;\;\;R \cdot \phi_2\\
                                                                                                                            
                                                                                                                            
                                                                                                                            \end{array}
                                                                                                                            \end{array}
                                                                                                                            
                                                                                                                            Derivation
                                                                                                                            1. Split input into 2 regimes
                                                                                                                            2. if phi2 < 8.7999999999999995e63

                                                                                                                              1. Initial program 64.1%

                                                                                                                                \[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)} \]
                                                                                                                              2. Add Preprocessing
                                                                                                                              3. Taylor expanded in phi1 around -inf

                                                                                                                                \[\leadsto R \cdot \color{blue}{\left(-1 \cdot \phi_1\right)} \]
                                                                                                                              4. Step-by-step derivation
                                                                                                                                1. Applied rewrites19.9%

                                                                                                                                  \[\leadsto R \cdot \color{blue}{\left(-\phi_1\right)} \]

                                                                                                                                if 8.7999999999999995e63 < phi2

                                                                                                                                1. Initial program 56.5%

                                                                                                                                  \[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)} \]
                                                                                                                                2. Add Preprocessing
                                                                                                                                3. Taylor expanded in phi2 around inf

                                                                                                                                  \[\leadsto R \cdot \color{blue}{\phi_2} \]
                                                                                                                                4. Step-by-step derivation
                                                                                                                                  1. Applied rewrites74.6%

                                                                                                                                    \[\leadsto R \cdot \color{blue}{\phi_2} \]
                                                                                                                                5. Recombined 2 regimes into one program.
                                                                                                                                6. Add Preprocessing

                                                                                                                                Alternative 14: 57.9% accurate, 25.4× speedup?

                                                                                                                                \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ \left(\left(-\phi_1\right) + \phi_2\right) \cdot R \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
                                                                                                                                 (* (+ (- phi1) phi2) R))
                                                                                                                                assert(R < lambda1 && lambda1 < lambda2 && lambda2 < phi1 && phi1 < phi2);
                                                                                                                                double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
                                                                                                                                	return (-phi1 + phi2) * R;
                                                                                                                                }
                                                                                                                                
                                                                                                                                NOTE: R, lambda1, lambda2, phi1, and phi2 should be sorted in increasing order before calling this function.
                                                                                                                                module fmin_fmax_functions
                                                                                                                                    implicit none
                                                                                                                                    private
                                                                                                                                    public fmax
                                                                                                                                    public fmin
                                                                                                                                
                                                                                                                                    interface fmax
                                                                                                                                        module procedure fmax88
                                                                                                                                        module procedure fmax44
                                                                                                                                        module procedure fmax84
                                                                                                                                        module procedure fmax48
                                                                                                                                    end interface
                                                                                                                                    interface fmin
                                                                                                                                        module procedure fmin88
                                                                                                                                        module procedure fmin44
                                                                                                                                        module procedure fmin84
                                                                                                                                        module procedure fmin48
                                                                                                                                    end interface
                                                                                                                                contains
                                                                                                                                    real(8) function fmax88(x, y) result (res)
                                                                                                                                        real(8), intent (in) :: x
                                                                                                                                        real(8), intent (in) :: y
                                                                                                                                        res = merge(y, merge(x, max(x, y), y /= y), x /= x)
                                                                                                                                    end function
                                                                                                                                    real(4) function fmax44(x, y) result (res)
                                                                                                                                        real(4), intent (in) :: x
                                                                                                                                        real(4), intent (in) :: y
                                                                                                                                        res = merge(y, merge(x, max(x, y), y /= y), x /= x)
                                                                                                                                    end function
                                                                                                                                    real(8) function fmax84(x, y) result(res)
                                                                                                                                        real(8), intent (in) :: x
                                                                                                                                        real(4), intent (in) :: y
                                                                                                                                        res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
                                                                                                                                    end function
                                                                                                                                    real(8) function fmax48(x, y) result(res)
                                                                                                                                        real(4), intent (in) :: x
                                                                                                                                        real(8), intent (in) :: y
                                                                                                                                        res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
                                                                                                                                    end function
                                                                                                                                    real(8) function fmin88(x, y) result (res)
                                                                                                                                        real(8), intent (in) :: x
                                                                                                                                        real(8), intent (in) :: y
                                                                                                                                        res = merge(y, merge(x, min(x, y), y /= y), x /= x)
                                                                                                                                    end function
                                                                                                                                    real(4) function fmin44(x, y) result (res)
                                                                                                                                        real(4), intent (in) :: x
                                                                                                                                        real(4), intent (in) :: y
                                                                                                                                        res = merge(y, merge(x, min(x, y), y /= y), x /= x)
                                                                                                                                    end function
                                                                                                                                    real(8) function fmin84(x, y) result(res)
                                                                                                                                        real(8), intent (in) :: x
                                                                                                                                        real(4), intent (in) :: y
                                                                                                                                        res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
                                                                                                                                    end function
                                                                                                                                    real(8) function fmin48(x, y) result(res)
                                                                                                                                        real(4), intent (in) :: x
                                                                                                                                        real(8), intent (in) :: y
                                                                                                                                        res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
                                                                                                                                    end function
                                                                                                                                end module
                                                                                                                                
                                                                                                                                real(8) function code(r, lambda1, lambda2, phi1, phi2)
                                                                                                                                use fmin_fmax_functions
                                                                                                                                    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
                                                                                                                                    code = (-phi1 + phi2) * r
                                                                                                                                end function
                                                                                                                                
                                                                                                                                assert R < lambda1 && lambda1 < lambda2 && lambda2 < phi1 && phi1 < phi2;
                                                                                                                                public static double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
                                                                                                                                	return (-phi1 + phi2) * R;
                                                                                                                                }
                                                                                                                                
                                                                                                                                [R, lambda1, lambda2, phi1, phi2] = sort([R, lambda1, lambda2, phi1, phi2])
                                                                                                                                def code(R, lambda1, lambda2, phi1, phi2):
                                                                                                                                	return (-phi1 + phi2) * R
                                                                                                                                
                                                                                                                                R, lambda1, lambda2, phi1, phi2 = sort([R, lambda1, lambda2, phi1, phi2])
                                                                                                                                function code(R, lambda1, lambda2, phi1, phi2)
                                                                                                                                	return Float64(Float64(Float64(-phi1) + phi2) * R)
                                                                                                                                end
                                                                                                                                
                                                                                                                                R, lambda1, lambda2, phi1, phi2 = num2cell(sort([R, lambda1, lambda2, phi1, phi2])){:}
                                                                                                                                function tmp = code(R, lambda1, lambda2, phi1, phi2)
                                                                                                                                	tmp = (-phi1 + phi2) * R;
                                                                                                                                end
                                                                                                                                
                                                                                                                                NOTE: R, lambda1, lambda2, phi1, and phi2 should be sorted in increasing order before calling this function.
                                                                                                                                code[R_, lambda1_, lambda2_, phi1_, phi2_] := N[(N[((-phi1) + phi2), $MachinePrecision] * R), $MachinePrecision]
                                                                                                                                
                                                                                                                                \begin{array}{l}
                                                                                                                                [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\
                                                                                                                                \\
                                                                                                                                \left(\left(-\phi_1\right) + \phi_2\right) \cdot R
                                                                                                                                \end{array}
                                                                                                                                
                                                                                                                                Derivation
                                                                                                                                1. Initial program 62.8%

                                                                                                                                  \[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)} \]
                                                                                                                                2. Add Preprocessing
                                                                                                                                3. Taylor expanded in phi1 around -inf

                                                                                                                                  \[\leadsto \color{blue}{-1 \cdot \left(\phi_1 \cdot \left(R + -1 \cdot \frac{R \cdot \phi_2}{\phi_1}\right)\right)} \]
                                                                                                                                4. Step-by-step derivation
                                                                                                                                  1. Applied rewrites29.2%

                                                                                                                                    \[\leadsto \color{blue}{\left(-\phi_1\right) \cdot \mathsf{fma}\left(-R, \frac{\phi_2}{\phi_1}, R\right)} \]
                                                                                                                                  2. Taylor expanded in phi1 around 0

                                                                                                                                    \[\leadsto \left(-\phi_1\right) \cdot \frac{-1 \cdot \left(R \cdot \phi_2\right) + R \cdot \phi_1}{\color{blue}{\phi_1}} \]
                                                                                                                                  3. Step-by-step derivation
                                                                                                                                    1. Applied rewrites32.3%

                                                                                                                                      \[\leadsto \left(-\phi_1\right) \cdot \frac{\mathsf{fma}\left(-R, \phi_2, \phi_1 \cdot R\right)}{\color{blue}{\phi_1}} \]
                                                                                                                                    2. Taylor expanded in phi1 around 0

                                                                                                                                      \[\leadsto -1 \cdot \left(R \cdot \phi_1\right) + \color{blue}{R \cdot \phi_2} \]
                                                                                                                                    3. Step-by-step derivation
                                                                                                                                      1. Applied rewrites30.4%

                                                                                                                                        \[\leadsto \mathsf{fma}\left(-1, \phi_1, \phi_2\right) \cdot \color{blue}{R} \]
                                                                                                                                      2. Step-by-step derivation
                                                                                                                                        1. Applied rewrites30.4%

                                                                                                                                          \[\leadsto \left(\left(-\phi_1\right) + \phi_2\right) \cdot R \]
                                                                                                                                        2. Add Preprocessing

                                                                                                                                        Alternative 15: 32.0% accurate, 46.5× speedup?

                                                                                                                                        \[\begin{array}{l} [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\ \\ R \cdot \phi_2 \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 (* R phi2))
                                                                                                                                        assert(R < lambda1 && lambda1 < lambda2 && lambda2 < phi1 && phi1 < phi2);
                                                                                                                                        double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
                                                                                                                                        	return R * phi2;
                                                                                                                                        }
                                                                                                                                        
                                                                                                                                        NOTE: R, lambda1, lambda2, phi1, and phi2 should be sorted in increasing order before calling this function.
                                                                                                                                        module fmin_fmax_functions
                                                                                                                                            implicit none
                                                                                                                                            private
                                                                                                                                            public fmax
                                                                                                                                            public fmin
                                                                                                                                        
                                                                                                                                            interface fmax
                                                                                                                                                module procedure fmax88
                                                                                                                                                module procedure fmax44
                                                                                                                                                module procedure fmax84
                                                                                                                                                module procedure fmax48
                                                                                                                                            end interface
                                                                                                                                            interface fmin
                                                                                                                                                module procedure fmin88
                                                                                                                                                module procedure fmin44
                                                                                                                                                module procedure fmin84
                                                                                                                                                module procedure fmin48
                                                                                                                                            end interface
                                                                                                                                        contains
                                                                                                                                            real(8) function fmax88(x, y) result (res)
                                                                                                                                                real(8), intent (in) :: x
                                                                                                                                                real(8), intent (in) :: y
                                                                                                                                                res = merge(y, merge(x, max(x, y), y /= y), x /= x)
                                                                                                                                            end function
                                                                                                                                            real(4) function fmax44(x, y) result (res)
                                                                                                                                                real(4), intent (in) :: x
                                                                                                                                                real(4), intent (in) :: y
                                                                                                                                                res = merge(y, merge(x, max(x, y), y /= y), x /= x)
                                                                                                                                            end function
                                                                                                                                            real(8) function fmax84(x, y) result(res)
                                                                                                                                                real(8), intent (in) :: x
                                                                                                                                                real(4), intent (in) :: y
                                                                                                                                                res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
                                                                                                                                            end function
                                                                                                                                            real(8) function fmax48(x, y) result(res)
                                                                                                                                                real(4), intent (in) :: x
                                                                                                                                                real(8), intent (in) :: y
                                                                                                                                                res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
                                                                                                                                            end function
                                                                                                                                            real(8) function fmin88(x, y) result (res)
                                                                                                                                                real(8), intent (in) :: x
                                                                                                                                                real(8), intent (in) :: y
                                                                                                                                                res = merge(y, merge(x, min(x, y), y /= y), x /= x)
                                                                                                                                            end function
                                                                                                                                            real(4) function fmin44(x, y) result (res)
                                                                                                                                                real(4), intent (in) :: x
                                                                                                                                                real(4), intent (in) :: y
                                                                                                                                                res = merge(y, merge(x, min(x, y), y /= y), x /= x)
                                                                                                                                            end function
                                                                                                                                            real(8) function fmin84(x, y) result(res)
                                                                                                                                                real(8), intent (in) :: x
                                                                                                                                                real(4), intent (in) :: y
                                                                                                                                                res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
                                                                                                                                            end function
                                                                                                                                            real(8) function fmin48(x, y) result(res)
                                                                                                                                                real(4), intent (in) :: x
                                                                                                                                                real(8), intent (in) :: y
                                                                                                                                                res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
                                                                                                                                            end function
                                                                                                                                        end module
                                                                                                                                        
                                                                                                                                        real(8) function code(r, lambda1, lambda2, phi1, phi2)
                                                                                                                                        use fmin_fmax_functions
                                                                                                                                            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
                                                                                                                                            code = r * phi2
                                                                                                                                        end function
                                                                                                                                        
                                                                                                                                        assert R < lambda1 && lambda1 < lambda2 && lambda2 < phi1 && phi1 < phi2;
                                                                                                                                        public static double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
                                                                                                                                        	return R * phi2;
                                                                                                                                        }
                                                                                                                                        
                                                                                                                                        [R, lambda1, lambda2, phi1, phi2] = sort([R, lambda1, lambda2, phi1, phi2])
                                                                                                                                        def code(R, lambda1, lambda2, phi1, phi2):
                                                                                                                                        	return R * phi2
                                                                                                                                        
                                                                                                                                        R, lambda1, lambda2, phi1, phi2 = sort([R, lambda1, lambda2, phi1, phi2])
                                                                                                                                        function code(R, lambda1, lambda2, phi1, phi2)
                                                                                                                                        	return Float64(R * phi2)
                                                                                                                                        end
                                                                                                                                        
                                                                                                                                        R, lambda1, lambda2, phi1, phi2 = num2cell(sort([R, lambda1, lambda2, phi1, phi2])){:}
                                                                                                                                        function tmp = code(R, lambda1, lambda2, phi1, phi2)
                                                                                                                                        	tmp = R * phi2;
                                                                                                                                        end
                                                                                                                                        
                                                                                                                                        NOTE: R, lambda1, lambda2, phi1, and phi2 should be sorted in increasing order before calling this function.
                                                                                                                                        code[R_, lambda1_, lambda2_, phi1_, phi2_] := N[(R * phi2), $MachinePrecision]
                                                                                                                                        
                                                                                                                                        \begin{array}{l}
                                                                                                                                        [R, lambda1, lambda2, phi1, phi2] = \mathsf{sort}([R, lambda1, lambda2, phi1, phi2])\\
                                                                                                                                        \\
                                                                                                                                        R \cdot \phi_2
                                                                                                                                        \end{array}
                                                                                                                                        
                                                                                                                                        Derivation
                                                                                                                                        1. Initial program 62.8%

                                                                                                                                          \[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)} \]
                                                                                                                                        2. Add Preprocessing
                                                                                                                                        3. Taylor expanded in phi2 around inf

                                                                                                                                          \[\leadsto R \cdot \color{blue}{\phi_2} \]
                                                                                                                                        4. Step-by-step derivation
                                                                                                                                          1. Applied rewrites17.3%

                                                                                                                                            \[\leadsto R \cdot \color{blue}{\phi_2} \]
                                                                                                                                          2. Add Preprocessing

                                                                                                                                          Reproduce

                                                                                                                                          ?
                                                                                                                                          herbie shell --seed 2025019 
                                                                                                                                          (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))))))