Equirectangular approximation to distance on a great circle

Average Error: 39.4 → 25.9

Time: 27.2s

Precision: binary64

[lambda1 lambda2]: =sort([lambda1 lambda2])

[phi1 phi2]: =sort([phi1 phi2])

Math

\[R \cdot \sqrt{\left(\left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\frac{\phi_1 + \phi_2}{2}\right)\right) \cdot \left(\left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\frac{\phi_1 + \phi_2}{2}\right)\right) + \left(\phi_1 - \phi_2\right) \cdot \left(\phi_1 - \phi_2\right)}\]

↓

\[\begin{array}{l} \mathbf{if}\;\lambda_1 - \lambda_2 \leq -1.2292846382901678 \cdot 10^{+289}:\\ \;\;\;\;R \cdot \left(\cos \left(0.5 \cdot \left(\phi_2 + \phi_1\right)\right) \cdot \left(\lambda_2 - \lambda_1\right)\right)\\ \mathbf{elif}\;\lambda_1 - \lambda_2 \leq -4.536558980401533 \cdot 10^{+249}:\\ \;\;\;\;R \cdot \left(\left(\lambda_1 - \lambda_2\right) \cdot \cos \left(0.5 \cdot \left(\phi_2 + \phi_1\right)\right)\right)\\ \mathbf{elif}\;\lambda_1 - \lambda_2 \leq -7.322681926544795 \cdot 10^{+217}:\\ \;\;\;\;R \cdot \left(\cos \left(0.5 \cdot \left(\phi_2 + \phi_1\right)\right) \cdot \left(\lambda_2 - \lambda_1\right)\right)\\ \mathbf{elif}\;\lambda_1 - \lambda_2 \leq -2.505868495132052 \cdot 10^{+207}:\\ \;\;\;\;R \cdot \left(\phi_2 - \phi_1\right)\\ \mathbf{elif}\;\lambda_1 - \lambda_2 \leq -4.328619352458468 \cdot 10^{+156}:\\ \;\;\;\;R \cdot \left(\cos \left(0.5 \cdot \left(\phi_2 + \phi_1\right)\right) \cdot \left(\lambda_2 - \lambda_1\right)\right)\\ \mathbf{elif}\;\lambda_1 - \lambda_2 \leq -7.325191687074718 \cdot 10^{+35}:\\ \;\;\;\;R \cdot \sqrt{\left(\left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\frac{\phi_2 + \phi_1}{2}\right)\right) \cdot \left(\left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\frac{\phi_2 + \phi_1}{2}\right)\right) + \left(\phi_1 - \phi_2\right) \cdot \left(\phi_1 - \phi_2\right)}\\ \mathbf{else}:\\ \;\;\;\;R \cdot \left(\phi_2 - \phi_1\right)\\ \end{array}\]

FPCore

(FPCore (R lambda1 lambda2 phi1 phi2)
 :precision binary64
 (*
  R
  (sqrt
   (+
    (*
     (* (- lambda1 lambda2) (cos (/ (+ phi1 phi2) 2.0)))
     (* (- lambda1 lambda2) (cos (/ (+ phi1 phi2) 2.0))))
    (* (- phi1 phi2) (- phi1 phi2))))))

↓

(FPCore (R lambda1 lambda2 phi1 phi2)
 :precision binary64
 (if (<= (- lambda1 lambda2) -1.2292846382901678e+289)
   (* R (* (cos (* 0.5 (+ phi2 phi1))) (- lambda2 lambda1)))
   (if (<= (- lambda1 lambda2) -4.536558980401533e+249)
     (* R (* (- lambda1 lambda2) (cos (* 0.5 (+ phi2 phi1)))))
     (if (<= (- lambda1 lambda2) -7.322681926544795e+217)
       (* R (* (cos (* 0.5 (+ phi2 phi1))) (- lambda2 lambda1)))
       (if (<= (- lambda1 lambda2) -2.505868495132052e+207)
         (* R (- phi2 phi1))
         (if (<= (- lambda1 lambda2) -4.328619352458468e+156)
           (* R (* (cos (* 0.5 (+ phi2 phi1))) (- lambda2 lambda1)))
           (if (<= (- lambda1 lambda2) -7.325191687074718e+35)
             (*
              R
              (sqrt
               (+
                (*
                 (* (- lambda1 lambda2) (cos (/ (+ phi2 phi1) 2.0)))
                 (* (- lambda1 lambda2) (cos (/ (+ phi2 phi1) 2.0))))
                (* (- phi1 phi2) (- phi1 phi2)))))
             (* R (- phi2 phi1)))))))))

C

double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
	return R * sqrt((((lambda1 - lambda2) * cos((phi1 + phi2) / 2.0)) * ((lambda1 - lambda2) * cos((phi1 + phi2) / 2.0))) + ((phi1 - phi2) * (phi1 - phi2)));
}

↓

double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
	double tmp;
	if ((lambda1 - lambda2) <= -1.2292846382901678e+289) {
		tmp = R * (cos(0.5 * (phi2 + phi1)) * (lambda2 - lambda1));
	} else if ((lambda1 - lambda2) <= -4.536558980401533e+249) {
		tmp = R * ((lambda1 - lambda2) * cos(0.5 * (phi2 + phi1)));
	} else if ((lambda1 - lambda2) <= -7.322681926544795e+217) {
		tmp = R * (cos(0.5 * (phi2 + phi1)) * (lambda2 - lambda1));
	} else if ((lambda1 - lambda2) <= -2.505868495132052e+207) {
		tmp = R * (phi2 - phi1);
	} else if ((lambda1 - lambda2) <= -4.328619352458468e+156) {
		tmp = R * (cos(0.5 * (phi2 + phi1)) * (lambda2 - lambda1));
	} else if ((lambda1 - lambda2) <= -7.325191687074718e+35) {
		tmp = R * sqrt((((lambda1 - lambda2) * cos((phi2 + phi1) / 2.0)) * ((lambda1 - lambda2) * cos((phi2 + phi1) / 2.0))) + ((phi1 - phi2) * (phi1 - phi2)));
	} else {
		tmp = R * (phi2 - phi1);
	}
	return tmp;
}

Error

Bits error versus R

Bits error versus lambda1

Bits error versus lambda2

Bits error versus phi1

Bits error versus phi2

Derivation

1. Split input into 4 regimes

2. if (-.f64 lambda1 lambda2) < -1.2292846382901678e289 or -4.5365589804015331e249 < (-.f64 lambda1 lambda2) < -7.32268192654479537e217 or -2.50586849513205195e207 < (-.f64 lambda1 lambda2) < -4.3286193524584683e156

   1. Initial program 64.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. Taylor expanded around -inf 34.3

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

   3. Simplified 34.3

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

   if -1.2292846382901678e289 < (-.f64 lambda1 lambda2) < -4.5365589804015331e249

   1. Initial program 64.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. Taylor expanded around inf 47.1

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

   3. Simplified 47.1

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

   if -7.32268192654479537e217 < (-.f64 lambda1 lambda2) < -2.50586849513205195e207 or -7.32519168707471779e35 < (-.f64 lambda1 lambda2)

   1. Initial program 25.4

      \[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. Taylor expanded around -inf 13.8

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

   if -4.3286193524584683e156 < (-.f64 lambda1 lambda2) < -7.32519168707471779e35

   1. Initial program 26.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)}\]
3. Recombined 4 regimes into one program.

4. Final simplification 25.9

   \[\leadsto \begin{array}{l} \mathbf{if}\;\lambda_1 - \lambda_2 \leq -1.2292846382901678 \cdot 10^{+289}:\\ \;\;\;\;R \cdot \left(\cos \left(0.5 \cdot \left(\phi_2 + \phi_1\right)\right) \cdot \left(\lambda_2 - \lambda_1\right)\right)\\ \mathbf{elif}\;\lambda_1 - \lambda_2 \leq -4.536558980401533 \cdot 10^{+249}:\\ \;\;\;\;R \cdot \left(\left(\lambda_1 - \lambda_2\right) \cdot \cos \left(0.5 \cdot \left(\phi_2 + \phi_1\right)\right)\right)\\ \mathbf{elif}\;\lambda_1 - \lambda_2 \leq -7.322681926544795 \cdot 10^{+217}:\\ \;\;\;\;R \cdot \left(\cos \left(0.5 \cdot \left(\phi_2 + \phi_1\right)\right) \cdot \left(\lambda_2 - \lambda_1\right)\right)\\ \mathbf{elif}\;\lambda_1 - \lambda_2 \leq -2.505868495132052 \cdot 10^{+207}:\\ \;\;\;\;R \cdot \left(\phi_2 - \phi_1\right)\\ \mathbf{elif}\;\lambda_1 - \lambda_2 \leq -4.328619352458468 \cdot 10^{+156}:\\ \;\;\;\;R \cdot \left(\cos \left(0.5 \cdot \left(\phi_2 + \phi_1\right)\right) \cdot \left(\lambda_2 - \lambda_1\right)\right)\\ \mathbf{elif}\;\lambda_1 - \lambda_2 \leq -7.325191687074718 \cdot 10^{+35}:\\ \;\;\;\;R \cdot \sqrt{\left(\left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\frac{\phi_2 + \phi_1}{2}\right)\right) \cdot \left(\left(\lambda_1 - \lambda_2\right) \cdot \cos \left(\frac{\phi_2 + \phi_1}{2}\right)\right) + \left(\phi_1 - \phi_2\right) \cdot \left(\phi_1 - \phi_2\right)}\\ \mathbf{else}:\\ \;\;\;\;R \cdot \left(\phi_2 - \phi_1\right)\\ \end{array}\]

Reproduce

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