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

   1 /* specfunc/expint.c
   2  *
   3  * Copyright (C) 2007 Brian Gough
   4  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002 Gerard Jungman
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19  */
  20
  21 /* Author: G. Jungman */
  22
  23 #include <config.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_errno.h>
  26 #include <gsl/gsl_sf_expint.h>
  27 #include <gsl/gsl_sf_gamma.h>
  28
  29 #include "error.h"
  30 #include "check.h"
  31
  32 #include "chebyshev.h"
  33 #include "cheb_eval.c"
  34
  35 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  36
  37 /*
  38  Chebyshev expansions: based on SLATEC e1.f, W. Fullerton
  39
  40  Series for AE11       on the interval -1.00000D-01 to  0.
  41                                         with weighted error   1.76E-17
  42                                          log weighted error  16.75
  43                                significant figures required  15.70
  44                                     decimal places required  17.55
  45
  46
  47  Series for AE12       on the interval -2.50000D-01 to -1.00000D-01
  48                                         with weighted error   5.83E-17
  49                                          log weighted error  16.23
  50                                significant figures required  15.76
  51                                     decimal places required  16.93
  52
  53
  54  Series for E11        on the interval -4.00000D+00 to -1.00000D+00
  55                                         with weighted error   1.08E-18
  56                                          log weighted error  17.97
  57                                significant figures required  19.02
  58                                     decimal places required  18.61
  59
  60
  61  Series for E12        on the interval -1.00000D+00 to  1.00000D+00
  62                                         with weighted error   3.15E-18
  63                                          log weighted error  17.50
  64                         approx significant figures required  15.8
  65                                     decimal places required  18.10
  66
  67
  68  Series for AE13       on the interval  2.50000D-01 to  1.00000D+00
  69                                         with weighted error   2.34E-17
  70                                          log weighted error  16.63
  71                                significant figures required  16.14
  72                                     decimal places required  17.33
  73
  74
  75  Series for AE14       on the interval  0.          to  2.50000D-01
  76                                         with weighted error   5.41E-17
  77                                          log weighted error  16.27
  78                                significant figures required  15.38
  79                                     decimal places required  16.97
  80 */
  81
  82 static double AE11_data[39] = {
  83    0.121503239716065790,
  84   -0.065088778513550150,
  85    0.004897651357459670,
  86   -0.000649237843027216,
  87    0.000093840434587471,
  88    0.000000420236380882,
  89   -0.000008113374735904,
  90    0.000002804247688663,
  91    0.000000056487164441,
  92   -0.000000344809174450,
  93    0.000000058209273578,
  94    0.000000038711426349,
  95   -0.000000012453235014,
  96   -0.000000005118504888,
  97    0.000000002148771527,
  98    0.000000000868459898,
  99   -0.000000000343650105,
 100   -0.000000000179796603,
 101    0.000000000047442060,
 102    0.000000000040423282,
 103   -0.000000000003543928,
 104   -0.000000000008853444,
 105   -0.000000000000960151,
 106    0.000000000001692921,
 107    0.000000000000607990,
 108   -0.000000000000224338,
 109   -0.000000000000200327,
 110   -0.000000000000006246,
 111    0.000000000000045571,
 112    0.000000000000016383,
 113   -0.000000000000005561,
 114   -0.000000000000006074,
 115   -0.000000000000000862,
 116    0.000000000000001223,
 117    0.000000000000000716,
 118   -0.000000000000000024,
 119   -0.000000000000000201,
 120   -0.000000000000000082,
 121    0.000000000000000017
 122 };
 123 static cheb_series AE11_cs = {
 124   AE11_data,
 125   38,
 126   -1, 1,
 127   20
 128 };
 129
 130 static double AE12_data[25] = {
 131    0.582417495134726740,
 132   -0.158348850905782750,
 133   -0.006764275590323141,
 134    0.005125843950185725,
 135    0.000435232492169391,
 136   -0.000143613366305483,
 137   -0.000041801320556301,
 138   -0.000002713395758640,
 139    0.000001151381913647,
 140    0.000000420650022012,
 141    0.000000066581901391,
 142    0.000000000662143777,
 143   -0.000000002844104870,
 144   -0.000000000940724197,
 145   -0.000000000177476602,
 146   -0.000000000015830222,
 147    0.000000000002905732,
 148    0.000000000001769356,
 149    0.000000000000492735,
 150    0.000000000000093709,
 151    0.000000000000010707,
 152   -0.000000000000000537,
 153   -0.000000000000000716,
 154   -0.000000000000000244,
 155   -0.000000000000000058
 156 };
 157 static cheb_series AE12_cs = {
 158   AE12_data,
 159   24,
 160   -1, 1,
 161   15
 162 };
 163
 164 static double E11_data[19] = {
 165   -16.11346165557149402600,
 166     7.79407277874268027690,
 167    -1.95540581886314195070,
 168     0.37337293866277945612,
 169    -0.05692503191092901938,
 170     0.00721107776966009185,
 171    -0.00078104901449841593,
 172     0.00007388093356262168,
 173    -0.00000620286187580820,
 174     0.00000046816002303176,
 175    -0.00000003209288853329,
 176     0.00000000201519974874,
 177    -0.00000000011673686816,
 178     0.00000000000627627066,
 179    -0.00000000000031481541,
 180     0.00000000000001479904,
 181    -0.00000000000000065457,
 182     0.00000000000000002733,
 183    -0.00000000000000000108
 184 };
 185 static cheb_series E11_cs = {
 186   E11_data,
 187   18,
 188   -1, 1,
 189   13
 190 };
 191
 192 static double E12_data[16] = {
 193   -0.03739021479220279500,
 194    0.04272398606220957700,
 195   -0.13031820798497005440,
 196    0.01441912402469889073,
 197   -0.00134617078051068022,
 198    0.00010731029253063780,
 199   -0.00000742999951611943,
 200    0.00000045377325690753,
 201   -0.00000002476417211390,
 202    0.00000000122076581374,
 203   -0.00000000005485141480,
 204    0.00000000000226362142,
 205   -0.00000000000008635897,
 206    0.00000000000000306291,
 207   -0.00000000000000010148,
 208    0.00000000000000000315
 209 };
 210 static cheb_series E12_cs = {
 211   E12_data,
 212   15,
 213   -1, 1,
 214   10
 215 };
 216
 217 static double AE13_data[25] = {
 218   -0.605773246640603460,
 219   -0.112535243483660900,
 220    0.013432266247902779,
 221   -0.001926845187381145,
 222    0.000309118337720603,
 223   -0.000053564132129618,
 224    0.000009827812880247,
 225   -0.000001885368984916,
 226    0.000000374943193568,
 227   -0.000000076823455870,
 228    0.000000016143270567,
 229   -0.000000003466802211,
 230    0.000000000758754209,
 231   -0.000000000168864333,
 232    0.000000000038145706,
 233   -0.000000000008733026,
 234    0.000000000002023672,
 235   -0.000000000000474132,
 236    0.000000000000112211,
 237   -0.000000000000026804,
 238    0.000000000000006457,
 239   -0.000000000000001568,
 240    0.000000000000000383,
 241   -0.000000000000000094,
 242    0.000000000000000023
 243 };
 244 static cheb_series AE13_cs = {
 245   AE13_data,
 246   24,
 247   -1, 1,
 248   15
 249 };
 250
 251 static double AE14_data[26] = {
 252   -0.18929180007530170,
 253   -0.08648117855259871,
 254    0.00722410154374659,
 255   -0.00080975594575573,
 256    0.00010999134432661,
 257   -0.00001717332998937,
 258    0.00000298562751447,
 259   -0.00000056596491457,
 260    0.00000011526808397,
 261   -0.00000002495030440,
 262    0.00000000569232420,
 263   -0.00000000135995766,
 264    0.00000000033846628,
 265   -0.00000000008737853,
 266    0.00000000002331588,
 267   -0.00000000000641148,
 268    0.00000000000181224,
 269   -0.00000000000052538,
 270    0.00000000000015592,
 271   -0.00000000000004729,
 272    0.00000000000001463,
 273   -0.00000000000000461,
 274    0.00000000000000148,
 275   -0.00000000000000048,
 276    0.00000000000000016,
 277   -0.00000000000000005
 278 };
 279 static cheb_series AE14_cs = {
 280   AE14_data,
 281   25,
 282   -1, 1,
 283   13
 284 };
 285
 286
 287
 288 /* implementation for E1, allowing for scaling by exp(x) */
 289 static
 290 int expint_E1_impl(const double x, gsl_sf_result * result, const int scale)
 291 {
 292   const double xmaxt = -GSL_LOG_DBL_MIN;      /* XMAXT = -LOG (R1MACH(1)) */
 293   const double xmax  = xmaxt - log(xmaxt);    /* XMAX = XMAXT - LOG(XMAXT) */
 294
 295   /* CHECK_POINTER(result) */
 296
 297   if(x < -xmax && !scale) {
 298       OVERFLOW_ERROR(result);
 299   }
 300   else if(x <= -10.0) {
 301     const double s = 1.0/x * ( scale ? 1.0 : exp(-x) );
 302     gsl_sf_result result_c;
 303     cheb_eval_e(&AE11_cs, 20.0/x+1.0, &result_c);
 304     result->val  = s * (1.0 + result_c.val);
 305     result->err  = s * result_c.err;
 306     result->err += 2.0 * GSL_DBL_EPSILON * (fabs(x) + 1.0) * fabs(result->val);
 307     return GSL_SUCCESS;
 308   }
 309   else if(x <= -4.0) {
 310     const double s = 1.0/x * ( scale ? 1.0 : exp(-x) );
 311     gsl_sf_result result_c;
 312     cheb_eval_e(&AE12_cs, (40.0/x+7.0)/3.0, &result_c);
 313     result->val  = s * (1.0 + result_c.val);
 314     result->err  = s * result_c.err;
 315     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 316     return GSL_SUCCESS;
 317   }
 318   else if(x <= -1.0) {
 319     const double ln_term = -log(fabs(x));
 320     const double scale_factor = ( scale ? exp(x) : 1.0 );
 321     gsl_sf_result result_c;
 322     cheb_eval_e(&E11_cs, (2.0*x+5.0)/3.0, &result_c);
 323     result->val  = scale_factor * (ln_term + result_c.val);
 324     result->err  = scale_factor * (result_c.err + GSL_DBL_EPSILON * fabs(ln_term));
 325     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 326     return GSL_SUCCESS;
 327   }
 328   else if(x == 0.0) {
 329     DOMAIN_ERROR(result);
 330   }
 331   else if(x <= 1.0) {
 332     const double ln_term = -log(fabs(x));
 333     const double scale_factor = ( scale ? exp(x) : 1.0 );
 334     gsl_sf_result result_c;
 335     cheb_eval_e(&E12_cs, x, &result_c);
 336     result->val  = scale_factor * (ln_term - 0.6875 + x + result_c.val);
 337     result->err  = scale_factor * (result_c.err + GSL_DBL_EPSILON * fabs(ln_term));
 338     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 339     return GSL_SUCCESS;
 340   }
 341   else if(x <= 4.0) {
 342     const double s = 1.0/x * ( scale ? 1.0 : exp(-x) );
 343     gsl_sf_result result_c;
 344     cheb_eval_e(&AE13_cs, (8.0/x-5.0)/3.0, &result_c);
 345     result->val  = s * (1.0 + result_c.val);
 346     result->err  = s * result_c.err;
 347     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 348     return GSL_SUCCESS;
 349   }
 350   else if(x <= xmax || scale) {
 351     const double s = 1.0/x * ( scale ? 1.0 : exp(-x) );
 352     gsl_sf_result result_c;
 353     cheb_eval_e(&AE14_cs, 8.0/x-1.0, &result_c);
 354     result->val  = s * (1.0 +  result_c.val);
 355     result->err  = s * (GSL_DBL_EPSILON + result_c.err);
 356     result->err += 2.0 * (x + 1.0) * GSL_DBL_EPSILON * fabs(result->val);
 357     if(result->val == 0.0)
 358       UNDERFLOW_ERROR(result);
 359     else
 360       return GSL_SUCCESS;
 361   }
 362   else {
 363     UNDERFLOW_ERROR(result);
 364   }
 365 }
 366
 367
 368 static
 369 int expint_E2_impl(const double x, gsl_sf_result * result, const int scale)
 370 {
 371   const double xmaxt = -GSL_LOG_DBL_MIN;
 372   const double xmax  = xmaxt - log(xmaxt);
 373
 374   /* CHECK_POINTER(result) */
 375
 376   if(x < -xmax && !scale) {
 377     OVERFLOW_ERROR(result);
 378   }
 379   else if (x == 0.0) {
 380     result->val = (scale ? 1.0 : 1.0);
 381     result->err = 0.0;
 382     return GSL_SUCCESS;
 383   } else if(x < 100.0) {
 384     const double ex = ( scale ? 1.0 : exp(-x) );
 385     gsl_sf_result result_E1;
 386     int stat_E1 = expint_E1_impl(x, &result_E1, scale);
 387     result->val  = ex - x*result_E1.val;
 388     result->err  = GSL_DBL_EPSILON*ex + fabs(x) * result_E1.err;
 389     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 390     return stat_E1;
 391   }
 392   else if(x < xmax || scale) {
 393     const double s = ( scale ? 1.0 : exp(-x) );
 394     const double c1  = -2.0;
 395     const double c2  =  6.0;
 396     const double c3  = -24.0;
 397     const double c4  =  120.0;
 398     const double c5  = -720.0;
 399     const double c6  =  5040.0;
 400     const double c7  = -40320.0;
 401     const double c8  =  362880.0;
 402     const double c9  = -3628800.0;
 403     const double c10 =  39916800.0;
 404     const double c11 = -479001600.0;
 405     const double c12 =  6227020800.0;
 406     const double c13 = -87178291200.0;
 407     const double y = 1.0/x;
 408     const double sum6 = c6+y*(c7+y*(c8+y*(c9+y*(c10+y*(c11+y*(c12+y*c13))))));
 409     const double sum  = y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*sum6)))));
 410     result->val = s * (1.0 + sum)/x;
 411     result->err = 2.0 * (x + 1.0) * GSL_DBL_EPSILON * result->val;
 412     if(result->val == 0.0)
 413       UNDERFLOW_ERROR(result);
 414     else
 415       return GSL_SUCCESS;
 416   }
 417   else {
 418     UNDERFLOW_ERROR(result);
 419   }
 420 }
 421
 422 static
 423 int expint_En_impl(const int n, const double x, gsl_sf_result * result, const int scale)
 424 {
 425   if (n < 0) {
 426     DOMAIN_ERROR(result);
 427   } else if (n == 0) {
 428     if (x == 0) {
 429       DOMAIN_ERROR(result);
 430     } else {
 431       result->val = (scale ? 1.0 : exp(-x)) / x;
 432       result->err = 2 * GSL_DBL_EPSILON * fabs(result->val);
 433       CHECK_UNDERFLOW(result);
 434       return GSL_SUCCESS;
 435     }
 436   } else if (n == 1) {
 437     return expint_E1_impl(x, result, scale);
 438   } else if (n == 2) {
 439     return expint_E2_impl(x, result, scale);
 440   } else {
 441     if(x < 0) {
 442       DOMAIN_ERROR(result);
 443     }
 444     if (x == 0) {
 445       result->val = (scale ? exp(x) : 1 ) * (1/(n-1.0));
 446       result->err = 2 * GSL_DBL_EPSILON * fabs(result->val);
 447       CHECK_UNDERFLOW(result);
 448       return GSL_SUCCESS;
 449     } else {
 450       gsl_sf_result result_g;
 451       double prefactor = pow(x, n-1);
 452       int status = gsl_sf_gamma_inc_e (1-n, x, &result_g);
 453       double scale_factor = ( scale ? exp(x) : 1.0 );
 454       result->val = scale_factor * prefactor * result_g.val;
 455       result->err = 2 * GSL_DBL_EPSILON * fabs(result->val);
 456       result->err += 2 * fabs(scale_factor * prefactor * result_g.err);
 457       if (status == GSL_SUCCESS) CHECK_UNDERFLOW(result);
 458       return status;
 459     }
 460   }
 461 }
 462
 463 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 464
 465
 466 int gsl_sf_expint_E1_e(const double x, gsl_sf_result * result)
 467 {
 468   return expint_E1_impl(x, result, 0);
 469 }
 470
 471
 472 int gsl_sf_expint_E1_scaled_e(const double x, gsl_sf_result * result)
 473 {
 474   return expint_E1_impl(x, result, 1);
 475 }
 476
 477
 478 int gsl_sf_expint_E2_e(const double x, gsl_sf_result * result)
 479 {
 480   return expint_E2_impl(x, result, 0);
 481 }
 482
 483
 484 int gsl_sf_expint_E2_scaled_e(const double x, gsl_sf_result * result)
 485 {
 486   return expint_E2_impl(x, result, 1);
 487 }
 488
 489 int gsl_sf_expint_En_e(const int n, const double x, gsl_sf_result * result)
 490 {
 491   return expint_En_impl(n, x, result, 0);
 492 }
 493
 494
 495 int gsl_sf_expint_En_scaled_e(const int n, const double x, gsl_sf_result * result)
 496 {
 497   return expint_En_impl(n, x, result, 1);
 498 }
 499
 500
 501 int gsl_sf_expint_Ei_e(const double x, gsl_sf_result * result)
 502 {
 503   /* CHECK_POINTER(result) */
 504
 505   {
 506     int status = gsl_sf_expint_E1_e(-x, result);
 507     result->val = -result->val;
 508     return status;
 509   }
 510 }
 511
 512
 513 int gsl_sf_expint_Ei_scaled_e(const double x, gsl_sf_result * result)
 514 {
 515   /* CHECK_POINTER(result) */
 516
 517   {
 518     int status = gsl_sf_expint_E1_scaled_e(-x, result);
 519     result->val = -result->val;
 520     return status;
 521   }
 522 }
 523
 524
 525 #if 0
 526 static double recurse_En(int n, double x, double E1)
 527 {
 528   int i;
 529   double En = E1;
 530   double ex = exp(-x);
 531   for(i=2; i<=n; i++) {
 532     En = 1./(i-1) * (ex - x * En);
 533   }
 534   return En;
 535 }
 536 #endif
 537
 538
 539 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 540
 541 #include "eval.h"
 542
 543 double gsl_sf_expint_E1(const double x)
 544 {
 545   EVAL_RESULT(gsl_sf_expint_E1_e(x, &result));
 546 }
 547
 548 double gsl_sf_expint_E1_scaled(const double x)
 549 {
 550   EVAL_RESULT(gsl_sf_expint_E1_scaled_e(x, &result));
 551 }
 552
 553 double gsl_sf_expint_E2(const double x)
 554 {
 555   EVAL_RESULT(gsl_sf_expint_E2_e(x, &result));
 556 }
 557
 558 double gsl_sf_expint_E2_scaled(const double x)
 559 {
 560   EVAL_RESULT(gsl_sf_expint_E2_scaled_e(x, &result));
 561 }
 562
 563 double gsl_sf_expint_En(const int n, const double x)
 564 {
 565   EVAL_RESULT(gsl_sf_expint_En_e(n, x, &result));
 566 }
 567
 568 double gsl_sf_expint_En_scaled(const int n, const double x)
 569 {
 570   EVAL_RESULT(gsl_sf_expint_En_scaled_e(n, x, &result));
 571 }
 572
 573 double gsl_sf_expint_Ei(const double x)
 574 {
 575   EVAL_RESULT(gsl_sf_expint_Ei_e(x, &result));
 576 }
 577
 578 double gsl_sf_expint_Ei_scaled(const double x)
 579 {
 580   EVAL_RESULT(gsl_sf_expint_Ei_scaled_e(x, &result));
 581 }