haversine in specfunc

haversine in specfunc

[Richard Mathar]

Since the evaluation of 1-cos(x) close to x=0, 2*pi etc
suffers from numerical instability in the standard implementations,
it would be convenient to have a haversine implementation,
where the haversine is defined according to section 4.3.147 of
the Abramowitz/Stegun book. This could be generated from
the existing Chebyshev implementation of the cosine as proposed
in the patches of gsl_sf_trig.h and trig.c below, which implement
hav(x) for real-valued x.
[I also think that the error estimate for cos(x) for x around zero
is a factor of 2 too large, by comparison with the standard
Taylor series that is used. In addition, 0.25*M_PI would generally
be replaced by M_PI_4 .]

Richard J. Mathar, http://www.strw.leidenuniv.nl/~mathar
 
--- gsl_sf_trig.h.org	2006-03-24 21:45:20.000000000 +0100
+++ gsl_sf_trig.h	2006-03-24 21:46:04.000000000 +0100
@@ -51,6 +51,11 @@
 int gsl_sf_cos_e(double x, gsl_sf_result * result);
 double gsl_sf_cos(const double x);
 
+/* Haversine(x) with GSL semantics.
+ */
+int gsl_sf_hav_e(double x, gsl_sf_result * result);
+double gsl_sf_hav(const double x);
+
 
 /* Hypot(x,y) with GSL semantics.
  */
--- trig.c.org	2006-03-24 20:51:56.000000000 +0100
+++ trig.c	2006-03-24 22:29:07.000000000 +0100
@@ -180,7 +180,7 @@
     }
     else {
       double sgn_result = sgn_x;
-      double y = floor(abs_x/(0.25*M_PI));
+      double y = floor(abs_x/M_PI_4);
       int octant = y - ldexp(floor(ldexp(y,-3)),3);
       int stat_cs;
       double z;
@@ -247,12 +247,12 @@
     if(abs_x < GSL_ROOT4_DBL_EPSILON) {
       const double x2 = x*x;
       result->val = 1.0 - 0.5*x2;
-      result->err = fabs(x2*x2/12.0);
+      result->err = fabs(x2*x2/24.0);
       return GSL_SUCCESS;
     }
     else {
       double sgn_result = 1.0;
-      double y = floor(abs_x/(0.25*M_PI));
+      double y = floor(abs_x/M_PI_4);
       int octant = y - ldexp(floor(ldexp(y,-3)),3);
       int stat_cs;
       double z;
@@ -307,6 +307,83 @@
   }
 }
 
+/* Haversine added by Richard J. Mathar, 2006-03-24 */
+int
+gsl_sf_hav_e(double x, gsl_sf_result * result)
+{
+  /* CHECK_POINTER(result) */
+
+  {
+    const double P1 = 7.85398125648498535156e-1;
+    const double P2 = 3.77489470793079817668e-8;
+    const double P3 = 2.69515142907905952645e-15;
+
+    const double abs_x = fabs(x);
+
+    if(abs_x < GSL_ROOT4_DBL_EPSILON) {
+      const double x2 = x*x;
+      result->val = 0.25*x2*(1.0 - x2/12.0) ;
+      result->err = fabs(x2*x2*x2/1440.0);
+      return GSL_SUCCESS;
+    }
+    else {
+      double sgn_result = 1.0;
+      double y = floor(abs_x/M_PI_4);
+      int octant = y - ldexp(floor(ldexp(y,-3)),3);
+      int stat_cs;
+      double z;
+
+      if(GSL_IS_ODD(octant)) {
+        octant += 1;
+        octant &= 07;
+        y += 1.0;
+      }
+
+      if(octant > 3) {
+        sgn_result = -sgn_result;
+        octant -= 4;
+      }
+
+      if(octant > 1) {
+        sgn_result = -sgn_result;
+      }
+
+      z = ((abs_x - y * P1) - y * P2) - y * P3;
+
+      if(octant == 0) {
+        gsl_sf_result cos_cs_result;
+        const double t = 8.0*fabs(z)/M_PI - 1.0;
+        stat_cs = cheb_eval_e(&cos_cs, t, &cos_cs_result);
+        if ( sgn_result > 0. )
+            result->val = 0.25*z*z * (1.0 - z*z * cos_cs_result.val);
+        else
+            result->val = 1.0 - 0.25*z*z * (1.0 - z*z * cos_cs_result.val);
+      }
+      else { /* octant == 2 */
+        gsl_sf_result sin_cs_result;
+        const double t = 8.0*fabs(z)/M_PI - 1.0;
+        stat_cs = cheb_eval_e(&sin_cs, t, &sin_cs_result);
+        result->val = 0.5 - 0.5 * sgn_result * z *(1.0 + z*z * sin_cs_result.val) ;
+      }
+
+      if(abs_x > 1.0/GSL_DBL_EPSILON) {
+        result->err = result->val ;
+      }
+      else if(abs_x > 100.0/GSL_SQRT_DBL_EPSILON) {
+        result->err = 2.0 * abs_x * GSL_DBL_EPSILON * result->val ;
+      }
+      else if(abs_x > 0.1/GSL_SQRT_DBL_EPSILON) {
+        result->err = 2.0 * GSL_SQRT_DBL_EPSILON * result->val ;
+      }
+      else {
+        result->err = 2.0 * GSL_DBL_EPSILON * result->val ;
+      }
+
+      return stat_cs;
+    }
+  }
+}
+
 
 int
 gsl_sf_hypot_e(const double x, const double y, gsl_sf_result * result)
@@ -413,13 +490,13 @@
 
   if(zi > 60.0) {
     lszr->val = -M_LN2 + zi;
-    lszi->val =  0.5*M_PI - zr;
+    lszi->val =  M_PI_2 - zr;
     lszr->err = 2.0 * GSL_DBL_EPSILON * fabs(lszr->val);
     lszi->err = 2.0 * GSL_DBL_EPSILON * fabs(lszi->val);
   }
   else if(zi < -60.0) {
     lszr->val = -M_LN2 - zi;
-    lszi->val = -0.5*M_PI + zr;
+    lszi->val = -M_PI_2 + zr;
     lszr->err = 2.0 * GSL_DBL_EPSILON * fabs(lszr->val);
     lszi->err = 2.0 * GSL_DBL_EPSILON * fabs(lszi->val);
   }
@@ -718,6 +795,12 @@
   EVAL_RESULT(gsl_sf_cos_e(x, &result));
 }
 
+/* haversine added by Richard J. Mathar, 2006-03-24 */
+double gsl_sf_hav(const double x)
+{
+  EVAL_RESULT(gsl_sf_hav_e(x, &result));
+}
+
 double gsl_sf_hypot(const double x, const double y)
 {
   EVAL_RESULT(gsl_sf_hypot_e(x, y, &result));

Re: haversine in specfunc

[Brian Gough]

Richard Mathar writes:
 > Since the evaluation of 1-cos(x) close to x=0, 2*pi etc
 > suffers from numerical instability in the standard implementations,
 > it would be convenient to have a haversine implementation,
 > where the haversine is defined according to section 4.3.147 of
 > the Abramowitz/Stegun book.

Hello,

Thanks for the suggestions.  For this case I think a double angle
formula will do the job: 1-cos(x) = 2 sin^2(x/2)

 > I also think that the error estimate for cos(x) for x around zero
 > is a factor of 2 too large, by comparison with the standard
 > Taylor series that is used.

Generally there is an extra factor of 2 for safety in the error terms.
I'll add a note about that to the design document.

-- 
Brian Gough