Equirectangular approximation to distance on a great circle

Average Error: 38.8 → 28.6

Time: 12.6s

Precision: binary64

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}\;\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) \leq 3.0573610316182304 \cdot 10^{+306}:\\ \;\;\;\;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)}\\ \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) (cos (/ (+ phi1 phi2) 2.0)))
        (* (- lambda1 lambda2) (cos (/ (+ phi1 phi2) 2.0))))
       (* (- phi1 phi2) (- phi1 phi2)))
      3.0573610316182304e+306)
   (*
    R
    (sqrt
     (+
      (*
       (* (- lambda1 lambda2) (cos (/ (+ phi1 phi2) 2.0)))
       (* (- lambda1 lambda2) (cos (/ (+ phi1 phi2) 2.0))))
      (* (- phi1 phi2) (- phi1 phi2)))))
   (* R (- phi2 phi1))))

C

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

↓

double code(double R, double lambda1, double lambda2, double phi1, double phi2) {
	double tmp;
	if ((((double) (((double) (((double) (((double) (lambda1 - lambda2)) * ((double) cos((((double) (phi1 + phi2)) / 2.0))))) * ((double) (((double) (lambda1 - lambda2)) * ((double) cos((((double) (phi1 + phi2)) / 2.0))))))) + ((double) (((double) (phi1 - phi2)) * ((double) (phi1 - phi2)))))) <= 3.0573610316182304e+306)) {
		tmp = ((double) (R * ((double) sqrt(((double) (((double) (((double) (((double) (lambda1 - lambda2)) * ((double) cos((((double) (phi1 + phi2)) / 2.0))))) * ((double) (((double) (lambda1 - lambda2)) * ((double) cos((((double) (phi1 + phi2)) / 2.0))))))) + ((double) (((double) (phi1 - phi2)) * ((double) (phi1 - phi2))))))))));
	} else {
		tmp = ((double) (R * ((double) (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 2 regimes

2. if (+.f64 (*.f64 (*.f64 (-.f64 lambda1 lambda2) (cos.f64 (/.f64 (+.f64 phi1 phi2) 2))) (*.f64 (-.f64 lambda1 lambda2) (cos.f64 (/.f64 (+.f64 phi1 phi2) 2)))) (*.f64 (-.f64 phi1 phi2) (-.f64 phi1 phi2))) < 3.0573610316182304e306

   1. Initial program 1.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)}\]

   if 3.0573610316182304e306 < (+.f64 (*.f64 (*.f64 (-.f64 lambda1 lambda2) (cos.f64 (/.f64 (+.f64 phi1 phi2) 2))) (*.f64 (-.f64 lambda1 lambda2) (cos.f64 (/.f64 (+.f64 phi1 phi2) 2)))) (*.f64 (-.f64 phi1 phi2) (-.f64 phi1 phi2)))

   1. Initial program 63.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. Taylor expanded around 0 46.6

      \[\leadsto R \cdot \color{blue}{\left(\phi_2 - \phi_1\right)}\]
3. Recombined 2 regimes into one program.

4. Final simplification 28.6

   \[\leadsto \begin{array}{l} \mathbf{if}\;\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) \leq 3.0573610316182304 \cdot 10^{+306}:\\ \;\;\;\;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)}\\ \mathbf{else}:\\ \;\;\;\;R \cdot \left(\phi_2 - \phi_1\right)\\ \end{array}\]

Reproduce

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