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

   1 /* specfunc/erfc.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003 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:  J. Theiler (modifications by G. Jungman) */
  21
  22 /*
  23  * See Hart et al, Computer Approximations, John Wiley and Sons, New York (1968)
  24  * (This applies only to the erfc8 stuff, which is the part
  25  *  of the original code that survives. I have replaced much of
  26  *  the other stuff with Chebyshev fits. These are simpler and
  27  *  more precise than the original approximations. [GJ])
  28  */
  29 #include <config.h>
  30 #include <gsl/gsl_math.h>
  31 #include <gsl/gsl_errno.h>
  32 #include <gsl/gsl_sf_exp.h>
  33 #include <gsl/gsl_sf_erf.h>
  34
  35 #include "check.h"
  36
  37 #include "chebyshev.h"
  38 #include "cheb_eval.c"
  39
  40 #define LogRootPi_  0.57236494292470008706
  41
  42
  43 static double erfc8_sum(double x)
  44 {
  45   /* estimates erfc(x) valid for 8 < x < 100 */
  46   /* This is based on index 5725 in Hart et al */
  47
  48   static double P[] = {
  49       2.97886562639399288862,
  50       7.409740605964741794425,
  51       6.1602098531096305440906,
  52       5.019049726784267463450058,
  53       1.275366644729965952479585264,
  54       0.5641895835477550741253201704
  55   };
  56   static double Q[] = {
  57       3.3690752069827527677,
  58       9.608965327192787870698,
  59       17.08144074746600431571095,
  60       12.0489519278551290360340491,
  61       9.396034016235054150430579648,
  62       2.260528520767326969591866945,
  63       1.0
  64   };
  65   double num=0.0, den=0.0;
  66   int i;
  67
  68   num = P[5];
  69   for (i=4; i>=0; --i) {
  70       num = x*num + P[i];
  71   }
  72   den = Q[6];
  73   for (i=5; i>=0; --i) {
  74       den = x*den + Q[i];
  75   }
  76
  77   return num/den;
  78 }
  79
  80 inline
  81 static double erfc8(double x)
  82 {
  83   double e;
  84   e = erfc8_sum(x);
  85   e *= exp(-x*x);
  86   return e;
  87 }
  88
  89 inline
  90 static double log_erfc8(double x)
  91 {
  92   double e;
  93   e = erfc8_sum(x);
  94   e = log(e) - x*x;
  95   return e;
  96 }
  97
  98 #if 0
  99 /* Abramowitz+Stegun, 7.2.14 */
 100 static double erfcasympsum(double x)
 101 {
 102   int i;
 103   double e = 1.;
 104   double coef = 1.;
 105   for (i=1; i<5; ++i) {
 106     /* coef *= -(2*i-1)/(2*x*x); ??? [GJ] */
 107     coef *= -(2*i+1)/(i*(4*x*x*x*x));
 108     e += coef;
 109     /*
 110     if (fabs(coef) < 1.0e-15) break;
 111     if (fabs(coef) > 1.0e10) break;
 112
 113     [GJ]: These tests are not useful. This function is only
 114     used below. Took them out; they gum up the pipeline.
 115     */
 116   }
 117   return e;
 118 }
 119 #endif /* 0 */
 120
 121
 122 /* Abramowitz+Stegun, 7.1.5 */
 123 static int erfseries(double x, gsl_sf_result * result)
 124 {
 125   double coef = x;
 126   double e    = coef;
 127   double del;
 128   int k;
 129   for (k=1; k<30; ++k) {
 130     coef *= -x*x/k;
 131     del   = coef/(2.0*k+1.0);
 132     e += del;
 133   }
 134   result->val = 2.0 / M_SQRTPI * e;
 135   result->err = 2.0 / M_SQRTPI * (fabs(del) + GSL_DBL_EPSILON);
 136   return GSL_SUCCESS;
 137 }
 138
 139
 140 /* Chebyshev fit for erfc((t+1)/2), -1 < t < 1
 141  */
 142 static double erfc_xlt1_data[20] = {
 143   1.06073416421769980345174155056,
 144  -0.42582445804381043569204735291,
 145   0.04955262679620434040357683080,
 146   0.00449293488768382749558001242,
 147  -0.00129194104658496953494224761,
 148  -0.00001836389292149396270416979,
 149   0.00002211114704099526291538556,
 150  -5.23337485234257134673693179020e-7,
 151  -2.78184788833537885382530989578e-7,
 152   1.41158092748813114560316684249e-8,
 153   2.72571296330561699984539141865e-9,
 154  -2.06343904872070629406401492476e-10,
 155  -2.14273991996785367924201401812e-11,
 156   2.22990255539358204580285098119e-12,
 157   1.36250074650698280575807934155e-13,
 158  -1.95144010922293091898995913038e-14,
 159  -6.85627169231704599442806370690e-16,
 160   1.44506492869699938239521607493e-16,
 161   2.45935306460536488037576200030e-18,
 162  -9.29599561220523396007359328540e-19
 163 };
 164 static cheb_series erfc_xlt1_cs = {
 165   erfc_xlt1_data,
 166   19,
 167   -1, 1,
 168   12
 169 };
 170
 171 /* Chebyshev fit for erfc(x) exp(x^2), 1 < x < 5, x = 2t + 3, -1 < t < 1
 172  */
 173 static double erfc_x15_data[25] = {
 174   0.44045832024338111077637466616,
 175  -0.143958836762168335790826895326,
 176   0.044786499817939267247056666937,
 177  -0.013343124200271211203618353102,
 178   0.003824682739750469767692372556,
 179  -0.001058699227195126547306482530,
 180   0.000283859419210073742736310108,
 181  -0.000073906170662206760483959432,
 182   0.000018725312521489179015872934,
 183  -4.62530981164919445131297264430e-6,
 184   1.11558657244432857487884006422e-6,
 185  -2.63098662650834130067808832725e-7,
 186   6.07462122724551777372119408710e-8,
 187  -1.37460865539865444777251011793e-8,
 188   3.05157051905475145520096717210e-9,
 189  -6.65174789720310713757307724790e-10,
 190   1.42483346273207784489792999706e-10,
 191  -3.00141127395323902092018744545e-11,
 192   6.22171792645348091472914001250e-12,
 193  -1.26994639225668496876152836555e-12,
 194   2.55385883033257575402681845385e-13,
 195  -5.06258237507038698392265499770e-14,
 196   9.89705409478327321641264227110e-15,
 197  -1.90685978789192181051961024995e-15,
 198   3.50826648032737849245113757340e-16
 199 };
 200 static cheb_series erfc_x15_cs = {
 201   erfc_x15_data,
 202   24,
 203   -1, 1,
 204   16
 205 };
 206
 207 /* Chebyshev fit for erfc(x) x exp(x^2), 5 < x < 10, x = (5t + 15)/2, -1 < t < 1
 208  */
 209 static double erfc_x510_data[20] = {
 210   1.11684990123545698684297865808,
 211   0.003736240359381998520654927536,
 212  -0.000916623948045470238763619870,
 213   0.000199094325044940833965078819,
 214  -0.000040276384918650072591781859,
 215   7.76515264697061049477127605790e-6,
 216  -1.44464794206689070402099225301e-6,
 217   2.61311930343463958393485241947e-7,
 218  -4.61833026634844152345304095560e-8,
 219   8.00253111512943601598732144340e-9,
 220  -1.36291114862793031395712122089e-9,
 221   2.28570483090160869607683087722e-10,
 222  -3.78022521563251805044056974560e-11,
 223   6.17253683874528285729910462130e-12,
 224  -9.96019290955316888445830597430e-13,
 225   1.58953143706980770269506726000e-13,
 226  -2.51045971047162509999527428316e-14,
 227   3.92607828989125810013581287560e-15,
 228  -6.07970619384160374392535453420e-16,
 229   9.12600607264794717315507477670e-17
 230 };
 231 static cheb_series erfc_x510_cs = {
 232   erfc_x510_data,
 233   19,
 234   -1, 1,
 235   12
 236 };
 237
 238 #if 0
 239 inline
 240 static double
 241 erfc_asymptotic(double x)
 242 {
 243   return exp(-x*x)/x * erfcasympsum(x) / M_SQRTPI;
 244 }
 245 inline
 246 static double
 247 log_erfc_asymptotic(double x)
 248 {
 249   return log(erfcasympsum(x)/x) - x*x - LogRootPi_;
 250 }
 251 #endif /* 0 */
 252
 253
 254 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 255
 256 int gsl_sf_erfc_e(double x, gsl_sf_result * result)
 257 {
 258   const double ax = fabs(x);
 259   double e_val, e_err;
 260
 261   /* CHECK_POINTER(result) */
 262
 263   if(ax <= 1.0) {
 264     double t = 2.0*ax - 1.0;
 265     gsl_sf_result c;
 266     cheb_eval_e(&erfc_xlt1_cs, t, &c);
 267     e_val = c.val;
 268     e_err = c.err;
 269   }
 270   else if(ax <= 5.0) {
 271     double ex2 = exp(-x*x);
 272     double t = 0.5*(ax-3.0);
 273     gsl_sf_result c;
 274     cheb_eval_e(&erfc_x15_cs, t, &c);
 275     e_val = ex2 * c.val;
 276     e_err = ex2 * (c.err + 2.0*fabs(x)*GSL_DBL_EPSILON);
 277   }
 278   else if(ax < 10.0) {
 279     double exterm = exp(-x*x) / ax;
 280     double t = (2.0*ax - 15.0)/5.0;
 281     gsl_sf_result c;
 282     cheb_eval_e(&erfc_x510_cs, t, &c);
 283     e_val = exterm * c.val;
 284     e_err = exterm * (c.err + 2.0*fabs(x)*GSL_DBL_EPSILON + GSL_DBL_EPSILON);
 285   }
 286   else {
 287     e_val = erfc8(ax);
 288     e_err = (x*x + 1.0) * GSL_DBL_EPSILON * fabs(e_val);
 289   }
 290
 291   if(x < 0.0) {
 292     result->val  = 2.0 - e_val;
 293     result->err  = e_err;
 294     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 295   }
 296   else {
 297     result->val  = e_val;
 298     result->err  = e_err;
 299     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 300   }
 301
 302   return GSL_SUCCESS;
 303 }
 304
 305
 306 int gsl_sf_log_erfc_e(double x, gsl_sf_result * result)
 307 {
 308   /* CHECK_POINTER(result) */
 309
 310   if(x*x < 10.0*GSL_ROOT6_DBL_EPSILON) {
 311     const double y = x / M_SQRTPI;
 312     /* series for -1/2 Log[Erfc[Sqrt[Pi] y]] */
 313     const double c3 = (4.0 - M_PI)/3.0;
 314     const double c4 = 2.0*(1.0 - M_PI/3.0);
 315     const double c5 = -0.001829764677455021;  /* (96.0 - 40.0*M_PI + 3.0*M_PI*M_PI)/30.0  */
 316     const double c6 =  0.02629651521057465;   /* 2.0*(120.0 - 60.0*M_PI + 7.0*M_PI*M_PI)/45.0 */
 317     const double c7 = -0.01621575378835404;
 318     const double c8 =  0.00125993961762116;
 319     const double c9 =  0.00556964649138;
 320     const double c10 = -0.0045563339802;
 321     const double c11 =  0.0009461589032;
 322     const double c12 =  0.0013200243174;
 323     const double c13 = -0.00142906;
 324     const double c14 =  0.00048204;
 325     double series = c8 + y*(c9 + y*(c10 + y*(c11 + y*(c12 + y*(c13 + c14*y)))));
 326     series = y*(1.0 + y*(1.0 + y*(c3 + y*(c4 + y*(c5 + y*(c6 + y*(c7 + y*series)))))));
 327     result->val = -2.0 * series;
 328     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 329     return GSL_SUCCESS;
 330   }
 331   /*
 332   don't like use of log1p(); added above series stuff for small x instead, should be ok [GJ]
 333   else if (fabs(x) < 1.0) {
 334     gsl_sf_result result_erf;
 335     gsl_sf_erf_e(x, &result_erf);
 336     result->val  = log1p(-result_erf.val);
 337     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 338     return GSL_SUCCESS;
 339   }
 340   */
 341   else if(x > 8.0) {
 342     result->val = log_erfc8(x);
 343     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 344     return GSL_SUCCESS;
 345   }
 346   else {
 347     gsl_sf_result result_erfc;
 348     gsl_sf_erfc_e(x, &result_erfc);
 349     result->val  = log(result_erfc.val);
 350     result->err  = fabs(result_erfc.err / result_erfc.val);
 351     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 352     return GSL_SUCCESS;
 353   }
 354 }
 355
 356
 357 int gsl_sf_erf_e(double x, gsl_sf_result * result)
 358 {
 359   /* CHECK_POINTER(result) */
 360
 361   if(fabs(x) < 1.0) {
 362     return erfseries(x, result);
 363   }
 364   else {
 365     gsl_sf_result result_erfc;
 366     gsl_sf_erfc_e(x, &result_erfc);
 367     result->val  = 1.0 - result_erfc.val;
 368     result->err  = result_erfc.err;
 369     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 370     return GSL_SUCCESS;
 371   }
 372 }
 373
 374
 375 int gsl_sf_erf_Z_e(double x, gsl_sf_result * result)
 376 {
 377   /* CHECK_POINTER(result) */
 378
 379   {
 380     const double ex2 = exp(-x*x/2.0);
 381     result->val  = ex2 / (M_SQRT2 * M_SQRTPI);
 382     result->err  = fabs(x * result->val) * GSL_DBL_EPSILON;
 383     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 384     CHECK_UNDERFLOW(result);
 385     return GSL_SUCCESS;
 386   }
 387 }
 388
 389
 390 int gsl_sf_erf_Q_e(double x, gsl_sf_result * result)
 391 {
 392   /* CHECK_POINTER(result) */
 393
 394   {
 395     gsl_sf_result result_erfc;
 396     int stat = gsl_sf_erfc_e(x/M_SQRT2, &result_erfc);
 397     result->val  = 0.5 * result_erfc.val;
 398     result->err  = 0.5 * result_erfc.err;
 399     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 400     return stat;
 401   }
 402 }
 403
 404
 405 int gsl_sf_hazard_e(double x, gsl_sf_result * result)
 406 {
 407   if(x < 25.0)
 408   {
 409     gsl_sf_result result_ln_erfc;
 410     const int stat_l = gsl_sf_log_erfc_e(x/M_SQRT2, &result_ln_erfc);
 411     const double lnc = -0.22579135264472743236; /* ln(sqrt(2/pi)) */
 412     const double arg = lnc - 0.5*x*x - result_ln_erfc.val;
 413     const int stat_e = gsl_sf_exp_e(arg, result);
 414     result->err += 3.0 * (1.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
 415     result->err += fabs(result_ln_erfc.err * result->val);
 416     return GSL_ERROR_SELECT_2(stat_l, stat_e);
 417   }
 418   else
 419   {
 420     const double ix2 = 1.0/(x*x);
 421     const double corrB = 1.0 - 9.0*ix2 * (1.0 - 11.0*ix2);
 422     const double corrM = 1.0 - 5.0*ix2 * (1.0 - 7.0*ix2 * corrB);
 423     const double corrT = 1.0 - ix2 * (1.0 - 3.0*ix2*corrM);
 424     result->val = x / corrT;
 425     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 426     return GSL_SUCCESS;
 427   }
 428 }
 429
 430
 431
 432 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 433
 434 #include "eval.h"
 435
 436 double gsl_sf_erfc(double x)
 437 {
 438   EVAL_RESULT(gsl_sf_erfc_e(x, &result));
 439 }
 440
 441 double gsl_sf_log_erfc(double x)
 442 {
 443   EVAL_RESULT(gsl_sf_log_erfc_e(x, &result));
 444 }
 445
 446 double gsl_sf_erf(double x)
 447 {
 448   EVAL_RESULT(gsl_sf_erf_e(x, &result));
 449 }
 450
 451 double gsl_sf_erf_Z(double x)
 452 {
 453   EVAL_RESULT(gsl_sf_erf_Z_e(x, &result));
 454 }
 455
 456 double gsl_sf_erf_Q(double x)
 457 {
 458   EVAL_RESULT(gsl_sf_erf_Q_e(x, &result));
 459 }
 460
 461 double gsl_sf_hazard(double x)
 462 {
 463   EVAL_RESULT(gsl_sf_hazard_e(x, &result));
 464 }
 465