htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/trig.c

   1 /* specfunc/trig.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman */
  21
  22 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_log.h>
  26 #include <gsl/gsl_sf_trig.h>
  27
  28 #include "error.h"
  29
  30 #include "chebyshev.h"
  31 #include "cheb_eval.c"
  32
  33 /* sinh(x) series
  34  * double-precision for |x| < 1.0
  35  */
  36 inline
  37 static
  38 int
  39 sinh_series(const double x, double * result)
  40 {
  41   const double y = x*x;
  42   const double c0 = 1.0/6.0;
  43   const double c1 = 1.0/120.0;
  44   const double c2 = 1.0/5040.0;
  45   const double c3 = 1.0/362880.0;
  46   const double c4 = 1.0/39916800.0;
  47   const double c5 = 1.0/6227020800.0;
  48   const double c6 = 1.0/1307674368000.0;
  49   const double c7 = 1.0/355687428096000.0;
  50   *result = x*(1.0 + y*(c0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*(c6+y*c7))))))));
  51   return GSL_SUCCESS;
  52 }
  53
  54
  55 /* cosh(x)-1 series
  56  * double-precision for |x| < 1.0
  57  */
  58 inline
  59 static
  60 int
  61 cosh_m1_series(const double x, double * result)
  62 {
  63   const double y = x*x;
  64   const double c0 = 0.5;
  65   const double c1 = 1.0/24.0;
  66   const double c2 = 1.0/720.0;
  67   const double c3 = 1.0/40320.0;
  68   const double c4 = 1.0/3628800.0;
  69   const double c5 = 1.0/479001600.0;
  70   const double c6 = 1.0/87178291200.0;
  71   const double c7 = 1.0/20922789888000.0;
  72   const double c8 = 1.0/6402373705728000.0;
  73   *result = y*(c0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*(c6+y*(c7+y*c8))))))));
  74   return GSL_SUCCESS;
  75 }
  76
  77
  78 /* Chebyshev expansion for f(t) = sinc((t+1)/2), -1 < t < 1
  79  */
  80 static double sinc_data[17] = {
  81   1.133648177811747875422,
  82  -0.532677564732557348781,
  83  -0.068293048346633177859,
  84   0.033403684226353715020,
  85   0.001485679893925747818,
  86  -0.000734421305768455295,
  87  -0.000016837282388837229,
  88   0.000008359950146618018,
  89   0.000000117382095601192,
  90  -0.000000058413665922724,
  91  -0.000000000554763755743,
  92   0.000000000276434190426,
  93   0.000000000001895374892,
  94  -0.000000000000945237101,
  95  -0.000000000000004900690,
  96   0.000000000000002445383,
  97   0.000000000000000009925
  98 };
  99 static cheb_series sinc_cs = {
 100   sinc_data,
 101   16,
 102   -1, 1,
 103   10
 104 };
 105
 106
 107 /* Chebyshev expansion for f(t) = g((t+1)Pi/8), -1<t<1
 108  * g(x) = (sin(x)/x - 1)/(x*x)
 109  */
 110 static double sin_data[12] = {
 111   -0.3295190160663511504173,
 112    0.0025374284671667991990,
 113    0.0006261928782647355874,
 114   -4.6495547521854042157541e-06,
 115   -5.6917531549379706526677e-07,
 116    3.7283335140973803627866e-09,
 117    3.0267376484747473727186e-10,
 118   -1.7400875016436622322022e-12,
 119   -1.0554678305790849834462e-13,
 120    5.3701981409132410797062e-16,
 121    2.5984137983099020336115e-17,
 122   -1.1821555255364833468288e-19
 123 };
 124 static cheb_series sin_cs = {
 125   sin_data,
 126   11,
 127   -1, 1,
 128   11
 129 };
 130
 131 /* Chebyshev expansion for f(t) = g((t+1)Pi/8), -1<t<1
 132  * g(x) = (2(cos(x) - 1)/(x^2) + 1) / x^2
 133  */
 134 static double cos_data[11] = {
 135   0.165391825637921473505668118136,
 136  -0.00084852883845000173671196530195,
 137  -0.000210086507222940730213625768083,
 138   1.16582269619760204299639757584e-6,
 139   1.43319375856259870334412701165e-7,
 140  -7.4770883429007141617951330184e-10,
 141  -6.0969994944584252706997438007e-11,
 142   2.90748249201909353949854872638e-13,
 143   1.77126739876261435667156490461e-14,
 144  -7.6896421502815579078577263149e-17,
 145  -3.7363121133079412079201377318e-18
 146 };
 147 static cheb_series cos_cs = {
 148   cos_data,
 149   10,
 150   -1, 1,
 151   10
 152 };
 153
 154
 155 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 156
 157 /* I would have prefered just using the library sin() function.
 158  * But after some experimentation I decided that there was
 159  * no good way to understand the error; library sin() is just a black box.
 160  * So we have to roll our own.
 161  */
 162 int
 163 gsl_sf_sin_e(double x, gsl_sf_result * result)
 164 {
 165   /* CHECK_POINTER(result) */
 166
 167   {
 168     const double P1 = 7.85398125648498535156e-1;
 169     const double P2 = 3.77489470793079817668e-8;
 170     const double P3 = 2.69515142907905952645e-15;
 171
 172     const double sgn_x = GSL_SIGN(x);
 173     const double abs_x = fabs(x);
 174
 175     if(abs_x < GSL_ROOT4_DBL_EPSILON) {
 176       const double x2 = x*x;
 177       result->val = x * (1.0 - x2/6.0);
 178       result->err = fabs(x*x2*x2 / 100.0);
 179       return GSL_SUCCESS;
 180     }
 181     else {
 182       double sgn_result = sgn_x;
 183       double y = floor(abs_x/(0.25*M_PI));
 184       int octant = y - ldexp(floor(ldexp(y,-3)),3);
 185       int stat_cs;
 186       double z;
 187
 188       if(GSL_IS_ODD(octant)) {
 189         octant += 1;
 190         octant &= 07;
 191         y += 1.0;
 192       }
 193
 194       if(octant > 3) {
 195         octant -= 4;
 196         sgn_result = -sgn_result;
 197       }
 198
 199       z = ((abs_x - y * P1) - y * P2) - y * P3;
 200
 201       if(octant == 0) {
 202         gsl_sf_result sin_cs_result;
 203         const double t = 8.0*fabs(z)/M_PI - 1.0;
 204         stat_cs = cheb_eval_e(&sin_cs, t, &sin_cs_result);
 205         result->val = z * (1.0 + z*z * sin_cs_result.val);
 206       }
 207       else { /* octant == 2 */
 208         gsl_sf_result cos_cs_result;
 209         const double t = 8.0*fabs(z)/M_PI - 1.0;
 210         stat_cs = cheb_eval_e(&cos_cs, t, &cos_cs_result);
 211         result->val = 1.0 - 0.5*z*z * (1.0 - z*z * cos_cs_result.val);
 212       }
 213
 214       result->val *= sgn_result;
 215
 216       if(abs_x > 1.0/GSL_DBL_EPSILON) {
 217         result->err = fabs(result->val);
 218       }
 219       else if(abs_x > 100.0/GSL_SQRT_DBL_EPSILON) {
 220         result->err = 2.0 * abs_x * GSL_DBL_EPSILON * fabs(result->val);
 221       }
 222       else if(abs_x > 0.1/GSL_SQRT_DBL_EPSILON) {
 223         result->err = 2.0 * GSL_SQRT_DBL_EPSILON * fabs(result->val);
 224       }
 225       else {
 226         result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 227       }
 228
 229       return stat_cs;
 230     }
 231   }
 232 }
 233
 234
 235 int
 236 gsl_sf_cos_e(double x, gsl_sf_result * result)
 237 {
 238   /* CHECK_POINTER(result) */
 239
 240   {
 241     const double P1 = 7.85398125648498535156e-1;
 242     const double P2 = 3.77489470793079817668e-8;
 243     const double P3 = 2.69515142907905952645e-15;
 244
 245     const double abs_x = fabs(x);
 246
 247     if(abs_x < GSL_ROOT4_DBL_EPSILON) {
 248       const double x2 = x*x;
 249       result->val = 1.0 - 0.5*x2;
 250       result->err = fabs(x2*x2/12.0);
 251       return GSL_SUCCESS;
 252     }
 253     else {
 254       double sgn_result = 1.0;
 255       double y = floor(abs_x/(0.25*M_PI));
 256       int octant = y - ldexp(floor(ldexp(y,-3)),3);
 257       int stat_cs;
 258       double z;
 259
 260       if(GSL_IS_ODD(octant)) {
 261         octant += 1;
 262         octant &= 07;
 263         y += 1.0;
 264       }
 265
 266       if(octant > 3) {
 267         octant -= 4;
 268         sgn_result = -sgn_result;
 269       }
 270
 271       if(octant > 1) {
 272         sgn_result = -sgn_result;
 273       }
 274
 275       z = ((abs_x - y * P1) - y * P2) - y * P3;
 276
 277       if(octant == 0) {
 278         gsl_sf_result cos_cs_result;
 279         const double t = 8.0*fabs(z)/M_PI - 1.0;
 280         stat_cs = cheb_eval_e(&cos_cs, t, &cos_cs_result);
 281         result->val = 1.0 - 0.5*z*z * (1.0 - z*z * cos_cs_result.val);
 282       }
 283       else { /* octant == 2 */
 284         gsl_sf_result sin_cs_result;
 285         const double t = 8.0*fabs(z)/M_PI - 1.0;
 286         stat_cs = cheb_eval_e(&sin_cs, t, &sin_cs_result);
 287         result->val = z * (1.0 + z*z * sin_cs_result.val);
 288       }
 289
 290       result->val *= sgn_result;
 291
 292       if(abs_x > 1.0/GSL_DBL_EPSILON) {
 293         result->err = fabs(result->val);
 294       }
 295       else if(abs_x > 100.0/GSL_SQRT_DBL_EPSILON) {
 296         result->err = 2.0 * abs_x * GSL_DBL_EPSILON * fabs(result->val);
 297       }
 298       else if(abs_x > 0.1/GSL_SQRT_DBL_EPSILON) {
 299         result->err = 2.0 * GSL_SQRT_DBL_EPSILON * fabs(result->val);
 300       }
 301       else {
 302         result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 303       }
 304
 305       return stat_cs;
 306     }
 307   }
 308 }
 309
 310
 311 int
 312 gsl_sf_hypot_e(const double x, const double y, gsl_sf_result * result)
 313 {
 314   /* CHECK_POINTER(result) */
 315
 316   if(x == 0.0 && y == 0.0) {
 317     result->val = 0.0;
 318     result->err = 0.0;
 319     return GSL_SUCCESS;
 320   }
 321   else {
 322     const double a = fabs(x);
 323     const double b = fabs(y);
 324     const double min = GSL_MIN_DBL(a,b);
 325     const double max = GSL_MAX_DBL(a,b);
 326     const double rat = min/max;
 327     const double root_term = sqrt(1.0 + rat*rat);
 328
 329     if(max < GSL_DBL_MAX/root_term) {
 330       result->val = max * root_term;
 331       result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 332       return GSL_SUCCESS;
 333     }
 334     else {
 335       OVERFLOW_ERROR(result);
 336     }
 337   }
 338 }
 339
 340
 341 int
 342 gsl_sf_complex_sin_e(const double zr, const double zi,
 343                         gsl_sf_result * szr, gsl_sf_result * szi)
 344 {
 345   /* CHECK_POINTER(szr) */
 346   /* CHECK_POINTER(szi) */
 347
 348   if(fabs(zi) < 1.0) {
 349     double ch_m1, sh;
 350     sinh_series(zi, &sh);
 351     cosh_m1_series(zi, &ch_m1);
 352     szr->val = sin(zr)*(ch_m1 + 1.0);
 353     szi->val = cos(zr)*sh;
 354     szr->err = 2.0 * GSL_DBL_EPSILON * fabs(szr->val);
 355     szi->err = 2.0 * GSL_DBL_EPSILON * fabs(szi->val);
 356     return GSL_SUCCESS;
 357   }
 358   else if(fabs(zi) < GSL_LOG_DBL_MAX) {
 359     double ex = exp(zi);
 360     double ch = 0.5*(ex+1.0/ex);
 361     double sh = 0.5*(ex-1.0/ex);
 362     szr->val = sin(zr)*ch;
 363     szi->val = cos(zr)*sh;
 364     szr->err = 2.0 * GSL_DBL_EPSILON * fabs(szr->val);
 365     szi->err = 2.0 * GSL_DBL_EPSILON * fabs(szi->val);
 366     return GSL_SUCCESS;
 367   }
 368   else {
 369     OVERFLOW_ERROR_2(szr, szi);
 370   }
 371 }
 372
 373
 374 int
 375 gsl_sf_complex_cos_e(const double zr, const double zi,
 376                         gsl_sf_result * czr, gsl_sf_result * czi)
 377 {
 378   /* CHECK_POINTER(czr) */
 379   /* CHECK_POINTER(czi) */
 380
 381   if(fabs(zi) < 1.0) {
 382     double ch_m1, sh;
 383     sinh_series(zi, &sh);
 384     cosh_m1_series(zi, &ch_m1);
 385     czr->val =  cos(zr)*(ch_m1 + 1.0);
 386     czi->val = -sin(zr)*sh;
 387     czr->err = 2.0 * GSL_DBL_EPSILON * fabs(czr->val);
 388     czi->err = 2.0 * GSL_DBL_EPSILON * fabs(czi->val);
 389     return GSL_SUCCESS;
 390   }
 391   else if(fabs(zi) < GSL_LOG_DBL_MAX) {
 392     double ex = exp(zi);
 393     double ch = 0.5*(ex+1.0/ex);
 394     double sh = 0.5*(ex-1.0/ex);
 395     czr->val =  cos(zr)*ch;
 396     czi->val = -sin(zr)*sh;
 397     czr->err = 2.0 * GSL_DBL_EPSILON * fabs(czr->val);
 398     czi->err = 2.0 * GSL_DBL_EPSILON * fabs(czi->val);
 399     return GSL_SUCCESS;
 400   }
 401   else {
 402     OVERFLOW_ERROR_2(czr,czi);
 403   }
 404 }
 405
 406
 407 int
 408 gsl_sf_complex_logsin_e(const double zr, const double zi,
 409                            gsl_sf_result * lszr, gsl_sf_result * lszi)
 410 {
 411   /* CHECK_POINTER(lszr) */
 412   /* CHECK_POINTER(lszi) */
 413
 414   if(zi > 60.0) {
 415     lszr->val = -M_LN2 + zi;
 416     lszi->val =  0.5*M_PI - zr;
 417     lszr->err = 2.0 * GSL_DBL_EPSILON * fabs(lszr->val);
 418     lszi->err = 2.0 * GSL_DBL_EPSILON * fabs(lszi->val);
 419   }
 420   else if(zi < -60.0) {
 421     lszr->val = -M_LN2 - zi;
 422     lszi->val = -0.5*M_PI + zr;
 423     lszr->err = 2.0 * GSL_DBL_EPSILON * fabs(lszr->val);
 424     lszi->err = 2.0 * GSL_DBL_EPSILON * fabs(lszi->val);
 425   }
 426   else {
 427     gsl_sf_result sin_r, sin_i;
 428     int status;
 429     gsl_sf_complex_sin_e(zr, zi, &sin_r, &sin_i); /* ok by construction */
 430     status = gsl_sf_complex_log_e(sin_r.val, sin_i.val, lszr, lszi);
 431     if(status == GSL_EDOM) {
 432       DOMAIN_ERROR_2(lszr, lszi);
 433     }
 434   }
 435   return gsl_sf_angle_restrict_symm_e(&(lszi->val));
 436 }
 437
 438
 439 int
 440 gsl_sf_lnsinh_e(const double x, gsl_sf_result * result)
 441 {
 442   /* CHECK_POINTER(result) */
 443
 444   if(x <= 0.0) {
 445     DOMAIN_ERROR(result);
 446   }
 447   else if(fabs(x) < 1.0) {
 448     double eps;
 449     sinh_series(x, &eps);
 450     result->val = log(eps);
 451     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 452     return GSL_SUCCESS;
 453   }
 454   else if(x < -0.5*GSL_LOG_DBL_EPSILON) {
 455     result->val = x + log(0.5*(1.0 - exp(-2.0*x)));
 456     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 457     return GSL_SUCCESS;
 458   }
 459   else {
 460     result->val = -M_LN2 + x;
 461     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 462     return GSL_SUCCESS;
 463   }
 464 }
 465
 466
 467 int gsl_sf_lncosh_e(const double x, gsl_sf_result * result)
 468 {
 469   /* CHECK_POINTER(result) */
 470
 471   if(fabs(x) < 1.0) {
 472     double eps;
 473     cosh_m1_series(x, &eps);
 474     return gsl_sf_log_1plusx_e(eps, result);
 475   }
 476   else if(x < -0.5*GSL_LOG_DBL_EPSILON) {
 477     result->val = x + log(0.5*(1.0 + exp(-2.0*x)));
 478     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 479     return GSL_SUCCESS;
 480   }
 481   else {
 482     result->val = -M_LN2 + x;
 483     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 484     return GSL_SUCCESS;
 485   }
 486 }
 487
 488
 489 /*
 490 inline int gsl_sf_sincos_e(const double theta, double * s, double * c)
 491 {
 492   double tan_half = tan(0.5 * theta);
 493   double den = 1. + tan_half*tan_half;
 494   double cos_theta = (1.0 - tan_half*tan_half) / den;
 495   double sin_theta = 2.0 * tan_half / den;
 496 }
 497 */
 498
 499 int
 500 gsl_sf_polar_to_rect(const double r, const double theta,
 501                           gsl_sf_result * x, gsl_sf_result * y)
 502 {
 503   double t   = theta;
 504   int status = gsl_sf_angle_restrict_symm_e(&t);
 505   double c = cos(t);
 506   double s = sin(t);
 507   x->val = r * cos(t);
 508   y->val = r * sin(t);
 509   x->err  = r * fabs(s * GSL_DBL_EPSILON * t);
 510   x->err += 2.0 * GSL_DBL_EPSILON * fabs(x->val);
 511   y->err  = r * fabs(c * GSL_DBL_EPSILON * t);
 512   y->err += 2.0 * GSL_DBL_EPSILON * fabs(y->val);
 513   return status;
 514 }
 515
 516
 517 int
 518 gsl_sf_rect_to_polar(const double x, const double y,
 519                           gsl_sf_result * r, gsl_sf_result * theta)
 520 {
 521   int stat_h = gsl_sf_hypot_e(x, y, r);
 522   if(r->val > 0.0) {
 523     theta->val = atan2(y, x);
 524     theta->err = 2.0 * GSL_DBL_EPSILON * fabs(theta->val);
 525     return stat_h;
 526   }
 527   else {
 528     DOMAIN_ERROR(theta);
 529   }
 530 }
 531
 532
 533 int gsl_sf_angle_restrict_symm_err_e(const double theta, gsl_sf_result * result)
 534 {
 535   /* synthetic extended precision constants */
 536   const double P1 = 4 * 7.8539812564849853515625e-01;
 537   const double P2 = 4 * 3.7748947079307981766760e-08;
 538   const double P3 = 4 * 2.6951514290790594840552e-15;
 539   const double TwoPi = 2*(P1 + P2 + P3);
 540
 541   const double y = GSL_SIGN(theta) * 2 * floor(fabs(theta)/TwoPi);
 542   double r = ((theta - y*P1) - y*P2) - y*P3;
 543
 544   if(r >  M_PI) { r = (((r-2*P1)-2*P2)-2*P3); }  /* r-TwoPi */
 545   else if (r < -M_PI) r = (((r+2*P1)+2*P2)+2*P3); /* r+TwoPi */
 546
 547   result->val = r;
 548
 549   if(fabs(theta) > 0.0625/GSL_DBL_EPSILON) {
 550     result->val = GSL_NAN;
 551     result->err = GSL_NAN;
 552     GSL_ERROR ("error", GSL_ELOSS);
 553   }
 554   else if(fabs(theta) > 0.0625/GSL_SQRT_DBL_EPSILON) {
 555     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val - theta);
 556     return GSL_SUCCESS;
 557   }
 558   else {
 559     double delta = fabs(result->val - theta);
 560     result->err = 2.0 * GSL_DBL_EPSILON * ((delta < M_PI) ? delta : M_PI);
 561     return GSL_SUCCESS;
 562   }
 563 }
 564
 565
 566 int gsl_sf_angle_restrict_pos_err_e(const double theta, gsl_sf_result * result)
 567 {
 568   /* synthetic extended precision constants */
 569   const double P1 = 4 * 7.85398125648498535156e-01;
 570   const double P2 = 4 * 3.77489470793079817668e-08;
 571   const double P3 = 4 * 2.69515142907905952645e-15;
 572   const double TwoPi = 2*(P1 + P2 + P3);
 573
 574   const double y = 2*floor(theta/TwoPi);
 575
 576   double r = ((theta - y*P1) - y*P2) - y*P3;
 577
 578   if(r > TwoPi) {r = (((r-2*P1)-2*P2)-2*P3); }  /* r-TwoPi */
 579   else if (r < 0) { /* may happen due to FP rounding */
 580     r = (((r+2*P1)+2*P2)+2*P3); /* r+TwoPi */
 581   }
 582
 583   result->val = r;
 584
 585   if(fabs(theta) > 0.0625/GSL_DBL_EPSILON) {
 586     result->val = GSL_NAN;
 587     result->err = fabs(result->val);
 588     GSL_ERROR ("error", GSL_ELOSS);
 589   }
 590   else if(fabs(theta) > 0.0625/GSL_SQRT_DBL_EPSILON) {
 591     result->err = GSL_DBL_EPSILON * fabs(result->val - theta);
 592     return GSL_SUCCESS;
 593   }
 594   else {
 595     double delta = fabs(result->val - theta);
 596     result->err = 2.0 * GSL_DBL_EPSILON * ((delta < M_PI) ? delta : M_PI);
 597     return GSL_SUCCESS;
 598   }
 599 }
 600
 601
 602 int gsl_sf_angle_restrict_symm_e(double * theta)
 603 {
 604   gsl_sf_result r;
 605   int stat = gsl_sf_angle_restrict_symm_err_e(*theta, &r);
 606   *theta = r.val;
 607   return stat;
 608 }
 609
 610
 611 int gsl_sf_angle_restrict_pos_e(double * theta)
 612 {
 613   gsl_sf_result r;
 614   int stat = gsl_sf_angle_restrict_pos_err_e(*theta, &r);
 615   *theta = r.val;
 616   return stat;
 617 }
 618
 619
 620 int gsl_sf_sin_err_e(const double x, const double dx, gsl_sf_result * result)
 621 {
 622   int stat_s = gsl_sf_sin_e(x, result);
 623   result->err += fabs(cos(x) * dx);
 624   result->err += GSL_DBL_EPSILON * fabs(result->val);
 625   return stat_s;
 626 }
 627
 628
 629 int gsl_sf_cos_err_e(const double x, const double dx, gsl_sf_result * result)
 630 {
 631   int stat_c = gsl_sf_cos_e(x, result);
 632   result->err += fabs(sin(x) * dx);
 633   result->err += GSL_DBL_EPSILON * fabs(result->val);
 634   return stat_c;
 635 }
 636
 637
 638 #if 0
 639 int
 640 gsl_sf_sin_pi_x_e(const double x, gsl_sf_result * result)
 641 {
 642   /* CHECK_POINTER(result) */
 643
 644   if(-100.0 < x && x < 100.0) {
 645     result->val = sin(M_PI * x) / (M_PI * x);
 646     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 647     return GSL_SUCCESS;
 648   }
 649   else {
 650     const double N = floor(x + 0.5);
 651     const double f = x - N;
 652
 653     if(N < INT_MAX && N > INT_MIN) {
 654       /* Make it an integer if we can. Saves another
 655        * call to floor().
 656        */
 657       const int intN    = (int)N;
 658       const double sign = ( GSL_IS_ODD(intN) ? -1.0 : 1.0 );
 659       result->val = sign * sin(M_PI * f);
 660       result->err = GSL_DBL_EPSILON * fabs(result->val);
 661     }
 662     else if(N > 2.0/GSL_DBL_EPSILON || N < -2.0/GSL_DBL_EPSILON) {
 663       /* All integer-valued floating point numbers
 664        * bigger than 2/eps=2^53 are actually even.
 665        */
 666       result->val = 0.0;
 667       result->err = 0.0;
 668     }
 669     else {
 670       const double resN = N - 2.0*floor(0.5*N); /* 0 for even N, 1 for odd N */
 671       const double sign = ( fabs(resN) > 0.5 ? -1.0 : 1.0 );
 672       result->val = sign * sin(M_PI*f);
 673       result->err = GSL_DBL_EPSILON * fabs(result->val);
 674     }
 675
 676     return GSL_SUCCESS;
 677   }
 678 }
 679 #endif
 680
 681
 682 int gsl_sf_sinc_e(double x, gsl_sf_result * result)
 683 {
 684   /* CHECK_POINTER(result) */
 685
 686   {
 687     const double ax = fabs(x);
 688
 689     if(ax < 0.8) {
 690       /* Do not go to the limit of the fit since
 691        * there is a zero there and the Chebyshev
 692        * accuracy will go to zero.
 693        */
 694       return cheb_eval_e(&sinc_cs, 2.0*ax-1.0, result);
 695     }
 696     else if(ax < 100.0) {
 697       /* Small arguments are no problem.
 698        * We trust the library sin() to
 699        * roughly machine precision.
 700        */
 701       result->val = sin(M_PI * ax)/(M_PI * ax);
 702       result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 703       return GSL_SUCCESS;
 704     }
 705     else {
 706       /* Large arguments must be handled separately.
 707        */
 708       const double r = M_PI*ax;
 709       gsl_sf_result s;
 710       int stat_s = gsl_sf_sin_e(r, &s);
 711       result->val = s.val/r;
 712       result->err = s.err/r + 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 713       return stat_s;
 714     }
 715   }
 716 }
 717
 718
 719
 720 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 721
 722 #include "eval.h"
 723
 724 double gsl_sf_sin(const double x)
 725 {
 726   EVAL_RESULT(gsl_sf_sin_e(x, &result));
 727 }
 728
 729 double gsl_sf_cos(const double x)
 730 {
 731   EVAL_RESULT(gsl_sf_cos_e(x, &result));
 732 }
 733
 734 double gsl_sf_hypot(const double x, const double y)
 735 {
 736   EVAL_RESULT(gsl_sf_hypot_e(x, y, &result));
 737 }
 738
 739 double gsl_sf_lnsinh(const double x)
 740 {
 741   EVAL_RESULT(gsl_sf_lnsinh_e(x, &result));
 742 }
 743
 744 double gsl_sf_lncosh(const double x)
 745 {
 746   EVAL_RESULT(gsl_sf_lncosh_e(x, &result));
 747 }
 748
 749 double gsl_sf_angle_restrict_symm(const double theta)
 750 {
 751   double result = theta;
 752   EVAL_DOUBLE(gsl_sf_angle_restrict_symm_e(&result));
 753 }
 754
 755 double gsl_sf_angle_restrict_pos(const double theta)
 756 {
 757   double result = theta;
 758   EVAL_DOUBLE(gsl_sf_angle_restrict_pos_e(&result));
 759 }
 760
 761 #if 0
 762 double gsl_sf_sin_pi_x(const double x)
 763 {
 764   EVAL_RESULT(gsl_sf_sin_pi_x_e(x, &result));
 765 }
 766 #endif
 767
 768 double gsl_sf_sinc(const double x)
 769 {
 770   EVAL_RESULT(gsl_sf_sinc_e(x, &result));
 771 }