This is the mail archive of the
gsl-discuss@sourceware.org
mailing list for the GSL project.
haversine in specfunc
- From: Richard Mathar <mathar at strw dot leidenuniv dot nl>
- To: gsl-discuss at sourceware dot org
- Cc: mathar at mail dot strw dot leidenuniv dot nl
- Date: Fri, 24 Mar 2006 22:53:51 +0100
- Subject: haversine in specfunc
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));