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

   1 /* specfunc/debye.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 /* augmented to n=5 and 6 2005-11-08 by R. J. Mathar, http://www.strw.leidenuniv.nl/~mathar */
  22
  23 #include <config.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_errno.h>
  26 #include <gsl/gsl_sf_debye.h>
  27
  28 #include "error.h"
  29 #include "check.h"
  30
  31 #include "chebyshev.h"
  32 #include "cheb_eval.c"
  33
  34 static double adeb1_data[17] = {
  35    2.4006597190381410194,
  36    0.1937213042189360089,
  37   -0.62329124554895770e-02,
  38    0.3511174770206480e-03,
  39   -0.228222466701231e-04,
  40    0.15805467875030e-05,
  41   -0.1135378197072e-06,
  42    0.83583361188e-08,
  43   -0.6264424787e-09,
  44    0.476033489e-10,
  45   -0.36574154e-11,
  46    0.2835431e-12,
  47   -0.221473e-13,
  48    0.17409e-14,
  49   -0.1376e-15,
  50    0.109e-16,
  51   -0.9e-18
  52 };
  53 static cheb_series adeb1_cs = {
  54   adeb1_data,
  55   16,
  56   -1.0, 1.0,
  57   9
  58 };
  59
  60 static double adeb2_data[18] = {
  61    2.5943810232570770282,
  62    0.2863357204530719834,
  63   -0.102062656158046713e-01,
  64    0.6049109775346844e-03,
  65   -0.405257658950210e-04,
  66    0.28633826328811e-05,
  67   -0.2086394303065e-06,
  68    0.155237875826e-07,
  69   -0.11731280087e-08,
  70    0.897358589e-10,
  71   -0.69317614e-11,
  72    0.5398057e-12,
  73   -0.423241e-13,
  74    0.33378e-14,
  75   -0.2645e-15,
  76    0.211e-16,
  77   -0.17e-17,
  78    0.1e-18
  79 };
  80 static cheb_series adeb2_cs = {
  81   adeb2_data,
  82   17,
  83   -1.0, 1.0,
  84   10
  85 };
  86
  87 static double adeb3_data[17] = {
  88    2.707737068327440945,
  89    0.340068135211091751,
  90   -0.12945150184440869e-01,
  91    0.7963755380173816e-03,
  92   -0.546360009590824e-04,
  93    0.39243019598805e-05,
  94   -0.2894032823539e-06,
  95    0.217317613962e-07,
  96   -0.16542099950e-08,
  97    0.1272796189e-09,
  98   -0.987963460e-11,
  99    0.7725074e-12,
 100   -0.607797e-13,
 101    0.48076e-14,
 102   -0.3820e-15,
 103    0.305e-16,
 104   -0.24e-17
 105 };
 106 static cheb_series adeb3_cs = {
 107   adeb3_data,
 108   16,
 109   -1.0, 1.0,
 110   10
 111 };
 112
 113 static double adeb4_data[17] = {
 114    2.781869415020523460,
 115    0.374976783526892863,
 116   -0.14940907399031583e-01,
 117    0.945679811437042e-03,
 118   -0.66132916138933e-04,
 119    0.4815632982144e-05,
 120   -0.3588083958759e-06,
 121    0.271601187416e-07,
 122   -0.20807099122e-08,
 123    0.1609383869e-09,
 124   -0.125470979e-10,
 125    0.9847265e-12,
 126   -0.777237e-13,
 127    0.61648e-14,
 128   -0.4911e-15,
 129    0.393e-16,
 130   -0.32e-17
 131 };
 132 static cheb_series adeb4_cs = {
 133   adeb4_data,
 134   16,
 135   -1.0, 1.0,
 136   10
 137 };
 138
 139 static double adeb5_data[17] = {
 140    2.8340269546834530149,
 141    0.3994098857106266445,
 142   -0.164566764773099646e-1,
 143    0.10652138340664541e-2,
 144   -0.756730374875418e-4,
 145    0.55745985240273e-5,
 146   -0.4190692330918e-6,
 147    0.319456143678e-7,
 148   -0.24613318171e-8,
 149    0.1912801633e-9,
 150   -0.149720049e-10,
 151    0.11790312e-11,
 152   -0.933329e-13,
 153    0.74218e-14,
 154   -0.5925e-15,
 155    0.475e-16,
 156   -0.39e-17
 157 };
 158 static cheb_series adeb5_cs = {
 159   adeb5_data,
 160   16,
 161   -1.0, 1.0,
 162   10
 163 };
 164
 165 static double adeb6_data[17] = {
 166   2.8726727134130122113,
 167   0.4174375352339027746,
 168  -0.176453849354067873e-1,
 169   0.11629852733494556e-2,
 170  -0.837118027357117e-4,
 171   0.62283611596189e-5,
 172  -0.4718644465636e-6,
 173   0.361950397806e-7,
 174  -0.28030368010e-8,
 175   0.2187681983e-9,
 176  -0.171857387e-10,
 177   0.13575809e-11,
 178  -0.1077580e-12,
 179   0.85893e-14,
 180  -0.6872e-15,
 181   0.552e-16,
 182  -0.44e-17
 183 };
 184 static cheb_series adeb6_cs = {
 185   adeb6_data,
 186   16,
 187   -1.0, 1.0,
 188   10
 189 };
 190
 191
 192 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 193
 194 int gsl_sf_debye_1_e(const double x, gsl_sf_result * result)
 195 {
 196   const double val_infinity = 1.64493406684822644;
 197   const double xcut = -GSL_LOG_DBL_MIN;
 198
 199   /* CHECK_POINTER(result) */
 200
 201   if(x < 0.0) {
 202     DOMAIN_ERROR(result);
 203   }
 204   else if(x < 2.0*GSL_SQRT_DBL_EPSILON) {
 205     result->val = 1.0 - 0.25*x + x*x/36.0;
 206     result->err = GSL_DBL_EPSILON * fabs(result->val);
 207     return GSL_SUCCESS;
 208   }
 209   else if(x <= 4.0) {
 210     const double t = x*x/8.0 - 1.0;
 211     gsl_sf_result c;
 212     cheb_eval_e(&adeb1_cs, t, &c);
 213     result->val = c.val - 0.25 * x;
 214     result->err = c.err + 0.25 * x * GSL_DBL_EPSILON;
 215     return GSL_SUCCESS;
 216   }
 217   else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
 218     const int nexp = floor(xcut/x);
 219     const double ex  = exp(-x);
 220     double sum = 0.0;
 221     double xk  = nexp * x;
 222     double rk  = nexp;
 223     int i;
 224     for(i=nexp; i>=1; i--) {
 225       sum *= ex;
 226       sum += (1.0 + 1.0/xk)/rk;
 227       rk -= 1.0;
 228       xk -= x;
 229     }
 230     result->val = val_infinity/x - sum*ex;
 231     result->err = GSL_DBL_EPSILON * fabs(result->val);
 232     return GSL_SUCCESS;
 233   }
 234   else if(x < xcut) {
 235     result->val = (val_infinity - exp(-x)*(x+1.0)) / x;
 236     result->err = GSL_DBL_EPSILON * fabs(result->val);
 237     return GSL_SUCCESS;
 238   }
 239   else {
 240     result->val = val_infinity/x;
 241     result->err = GSL_DBL_EPSILON * fabs(result->val);
 242     return GSL_SUCCESS;
 243   }
 244 }
 245
 246
 247 int gsl_sf_debye_2_e(const double x, gsl_sf_result * result)
 248 {
 249   const double val_infinity = 4.80822761263837714;
 250   const double xcut = -GSL_LOG_DBL_MIN;
 251
 252   /* CHECK_POINTER(result) */
 253
 254   if(x < 0.0) {
 255     DOMAIN_ERROR(result);
 256   }
 257   else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
 258     result->val = 1.0 - x/3.0 + x*x/24.0;
 259     result->err = GSL_DBL_EPSILON * result->val;
 260     return GSL_SUCCESS;
 261   }
 262   else if(x <= 4.0) {
 263     const double t = x*x/8.0 - 1.0;
 264     gsl_sf_result c;
 265     cheb_eval_e(&adeb2_cs, t, &c);
 266     result->val = c.val - x/3.0;
 267     result->err = c.err + GSL_DBL_EPSILON * x/3.0;
 268     return GSL_SUCCESS;
 269   }
 270   else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
 271     const int nexp = floor(xcut/x);
 272     const double ex  = exp(-x);
 273     double xk  = nexp * x;
 274     double rk  = nexp;
 275     double sum = 0.0;
 276     int i;
 277     for(i=nexp; i>=1; i--) {
 278       sum *= ex;
 279       sum += (1.0 + 2.0/xk + 2.0/(xk*xk)) / rk;
 280       rk -= 1.0;
 281       xk -= x;
 282     }
 283     result->val = val_infinity/(x*x) - 2.0 * sum * ex;
 284     result->err = GSL_DBL_EPSILON * fabs(result->val);
 285     return GSL_SUCCESS;
 286   }
 287   else if(x < xcut) {
 288     const double x2  = x*x;
 289     const double sum = 2.0 + 2.0*x + x2;
 290     result->val = (val_infinity - 2.0 * sum * exp(-x)) / x2;
 291     result->err = GSL_DBL_EPSILON * fabs(result->val);
 292     return GSL_SUCCESS;
 293   }
 294   else {
 295     result->val = (val_infinity/x)/x;
 296     result->err = GSL_DBL_EPSILON * result->val;
 297     CHECK_UNDERFLOW(result);
 298     return GSL_SUCCESS;
 299   }
 300 }
 301
 302
 303 int gsl_sf_debye_3_e(const double x, gsl_sf_result * result)
 304 {
 305   const double val_infinity = 19.4818182068004875;
 306   const double xcut = -GSL_LOG_DBL_MIN;
 307
 308   /* CHECK_POINTER(result) */
 309
 310   if(x < 0.0) {
 311     DOMAIN_ERROR(result);
 312   }
 313   else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
 314     result->val = 1.0 - 3.0*x/8.0 + x*x/20.0;
 315     result->err = GSL_DBL_EPSILON * result->val;
 316     return GSL_SUCCESS;
 317   }
 318   else if(x <= 4.0) {
 319     const double t = x*x/8.0 - 1.0;
 320     gsl_sf_result c;
 321     cheb_eval_e(&adeb3_cs, t, &c);
 322     result->val = c.val - 0.375*x;
 323     result->err = c.err + GSL_DBL_EPSILON * 0.375*x;
 324     return GSL_SUCCESS;
 325   }
 326   else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
 327     const int nexp = floor(xcut/x);
 328     const double ex  = exp(-x);
 329     double xk  = nexp * x;
 330     double rk  = nexp;
 331     double sum = 0.0;
 332     int i;
 333     for(i=nexp; i>=1; i--) {
 334       double xk_inv = 1.0/xk;
 335       sum *= ex;
 336       sum += (((6.0*xk_inv + 6.0)*xk_inv + 3.0)*xk_inv + 1.0) / rk;
 337       rk -= 1.0;
 338       xk -= x;
 339     }
 340     result->val = val_infinity/(x*x*x) - 3.0 * sum * ex;
 341     result->err = GSL_DBL_EPSILON * result->val;
 342     return GSL_SUCCESS;
 343   }
 344   else if(x < xcut) {
 345     const double x3 = x*x*x;
 346     const double sum = 6.0 + 6.0*x + 3.0*x*x + x3;
 347     result->val = (val_infinity - 3.0 * sum * exp(-x)) / x3;
 348     result->err = GSL_DBL_EPSILON * result->val;
 349     return GSL_SUCCESS;
 350   }
 351   else {
 352     result->val = ((val_infinity/x)/x)/x;
 353     result->err = GSL_DBL_EPSILON * result->val;
 354     CHECK_UNDERFLOW(result);
 355     return GSL_SUCCESS;
 356   }
 357 }
 358
 359
 360 int gsl_sf_debye_4_e(const double x, gsl_sf_result * result)
 361 {
 362   const double val_infinity = 99.5450644937635129;
 363   const double xcut = -GSL_LOG_DBL_MIN;
 364
 365   /* CHECK_POINTER(result) */
 366
 367   if(x < 0.0) {
 368     DOMAIN_ERROR(result);
 369   }
 370   else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
 371     result->val = 1.0 - 2.0*x/5.0 + x*x/18.0;
 372     result->err = GSL_DBL_EPSILON * result->val;
 373     return GSL_SUCCESS;
 374   }
 375   else if(x <= 4.0) {
 376     const double t = x*x/8.0 - 1.0;
 377     gsl_sf_result c;
 378     cheb_eval_e(&adeb4_cs, t, &c);
 379     result->val = c.val - 2.0*x/5.0;
 380     result->err = c.err + GSL_DBL_EPSILON * 2.0*x/5.0;
 381     return GSL_SUCCESS;
 382   }
 383   else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
 384     const int nexp = floor(xcut/x);
 385     const double ex  = exp(-x);
 386     double xk  = nexp * x;
 387     double rk  = nexp;
 388     double sum = 0.0;
 389     int i;
 390     for(i=nexp; i>=1; i--) {
 391       double xk_inv = 1.0/xk;
 392       sum *= ex;
 393       sum += ((((24.0*xk_inv + 24.0)*xk_inv + 12.0)*xk_inv + 4.0)*xk_inv + 1.0) / rk;
 394       rk -= 1.0;
 395       xk -= x;
 396     }
 397     result->val = val_infinity/(x*x*x*x) - 4.0 * sum * ex;
 398     result->err = GSL_DBL_EPSILON * result->val;
 399     return GSL_SUCCESS;
 400   }
 401   else if(x < xcut) {
 402     const double x2 = x*x;
 403     const double x4 = x2*x2;
 404     const double sum = 24.0 + 24.0*x + 12.0*x2 + 4.0*x2*x + x4;
 405     result->val = (val_infinity - 4.0 * sum * exp(-x)) / x4;
 406     result->err = GSL_DBL_EPSILON * result->val;
 407     return GSL_SUCCESS;
 408   }
 409   else {
 410     result->val = (((val_infinity/x)/x)/x)/x;
 411     result->err = GSL_DBL_EPSILON * result->val;
 412     CHECK_UNDERFLOW(result);
 413     return GSL_SUCCESS;
 414   }
 415 }
 416
 417 int gsl_sf_debye_5_e(const double x, gsl_sf_result * result)
 418 {
 419   const double val_infinity = 610.405837190669483828710757875 ;
 420   const double xcut = -GSL_LOG_DBL_MIN;
 421
 422   /* CHECK_POINTER(result) */
 423
 424   if(x < 0.0) {
 425     DOMAIN_ERROR(result);
 426   }
 427   else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
 428     result->val = 1.0 - 5.0*x/12.0 + 5.0*x*x/84.0;
 429     result->err = GSL_DBL_EPSILON * result->val;
 430     return GSL_SUCCESS;
 431   }
 432   else if(x <= 4.0) {
 433     const double t = x*x/8.0 - 1.0;
 434     gsl_sf_result c;
 435     cheb_eval_e(&adeb5_cs, t, &c);
 436     result->val = c.val - 5.0*x/12.0;
 437     result->err = c.err + GSL_DBL_EPSILON * 5.0*x/12.0;
 438     return GSL_SUCCESS;
 439   }
 440   else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
 441     const int nexp = floor(xcut/x);
 442     const double ex  = exp(-x);
 443     double xk  = nexp * x;
 444     double rk  = nexp;
 445     double sum = 0.0;
 446     int i;
 447     for(i=nexp; i>=1; i--) {
 448       double xk_inv = 1.0/xk;
 449       sum *= ex;
 450       sum += (((((120.0*xk_inv + 120.0)*xk_inv + 60.0)*xk_inv + 20.0)*xk_inv + 5.0)*xk_inv+ 1.0) / rk;
 451       rk -= 1.0;
 452       xk -= x;
 453     }
 454     result->val = val_infinity/(x*x*x*x*x) - 5.0 * sum * ex;
 455     result->err = GSL_DBL_EPSILON * result->val;
 456     return GSL_SUCCESS;
 457   }
 458   else if(x < xcut) {
 459     const double x2 = x*x;
 460     const double x4 = x2*x2;
 461     const double x5 = x4*x;
 462     const double sum = 120.0 + 120.0*x + 60.0*x2 + 20.0*x2*x + 5.0*x4 + x5;
 463     result->val = (val_infinity - 5.0 * sum * exp(-x)) / x5;
 464     result->err = GSL_DBL_EPSILON * result->val;
 465     return GSL_SUCCESS;
 466   }
 467   else {
 468     result->val = ((((val_infinity/x)/x)/x)/x)/x;
 469     result->err = GSL_DBL_EPSILON * result->val;
 470     CHECK_UNDERFLOW(result);
 471     return GSL_SUCCESS;
 472   }
 473 }
 474
 475 int gsl_sf_debye_6_e(const double x, gsl_sf_result * result)
 476 {
 477   const double val_infinity = 4356.06887828990661194792541535 ;
 478   const double xcut = -GSL_LOG_DBL_MIN;
 479
 480   /* CHECK_POINTER(result) */
 481
 482   if(x < 0.0) {
 483     DOMAIN_ERROR(result);
 484   }
 485   else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
 486     result->val = 1.0 - 3.0*x/7.0 + x*x/16.0;
 487     result->err = GSL_DBL_EPSILON * result->val;
 488     return GSL_SUCCESS;
 489   }
 490   else if(x <= 4.0) {
 491     const double t = x*x/8.0 - 1.0;
 492     gsl_sf_result c;
 493     cheb_eval_e(&adeb6_cs, t, &c);
 494     result->val = c.val - 3.0*x/7.0;
 495     result->err = c.err + GSL_DBL_EPSILON * 3.0*x/7.0;
 496     return GSL_SUCCESS;
 497   }
 498   else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
 499     const int nexp = floor(xcut/x);
 500     const double ex  = exp(-x);
 501     double xk  = nexp * x;
 502     double rk  = nexp;
 503     double sum = 0.0;
 504     int i;
 505     for(i=nexp; i>=1; i--) {
 506       double xk_inv = 1.0/xk;
 507       sum *= ex;
 508       sum += ((((((720.0*xk_inv + 720.0)*xk_inv + 360.0)*xk_inv + 120.0)*xk_inv + 30.0)*xk_inv+ 6.0)*xk_inv+ 1.0) / rk;
 509       rk -= 1.0;
 510       xk -= x;
 511     }
 512     result->val = val_infinity/(x*x*x*x*x*x) - 6.0 * sum * ex;
 513     result->err = GSL_DBL_EPSILON * result->val;
 514     return GSL_SUCCESS;
 515   }
 516   else if(x < xcut) {
 517     const double x2 = x*x;
 518     const double x4 = x2*x2;
 519     const double x6 = x4*x2;
 520     const double sum = 720.0 + 720.0*x + 360.0*x2 + 120.0*x2*x + 30.0*x4 + 6.0*x4*x +x6 ;
 521     result->val = (val_infinity - 6.0 * sum * exp(-x)) / x6;
 522     result->err = GSL_DBL_EPSILON * result->val;
 523     return GSL_SUCCESS;
 524   }
 525   else {
 526     result->val = (((((val_infinity/x)/x)/x)/x)/x)/x ;
 527     result->err = GSL_DBL_EPSILON * result->val;
 528     CHECK_UNDERFLOW(result);
 529     return GSL_SUCCESS;
 530   }
 531 }
 532
 533
 534 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 535
 536 #include "eval.h"
 537
 538 double gsl_sf_debye_1(const double x)
 539 {
 540   EVAL_RESULT(gsl_sf_debye_1_e(x, &result));
 541 }
 542
 543 double gsl_sf_debye_2(const double x)
 544 {
 545   EVAL_RESULT(gsl_sf_debye_2_e(x, &result));
 546 }
 547
 548 double gsl_sf_debye_3(const double x)
 549 {
 550   EVAL_RESULT(gsl_sf_debye_3_e(x, &result));
 551 }
 552
 553 double gsl_sf_debye_4(const double x)
 554 {
 555   EVAL_RESULT(gsl_sf_debye_4_e(x, &result));
 556 }
 557
 558 double gsl_sf_debye_5(const double x)
 559 {
 560   EVAL_RESULT(gsl_sf_debye_5_e(x, &result));
 561 }
 562
 563 double gsl_sf_debye_6(const double x)
 564 {
 565   EVAL_RESULT(gsl_sf_debye_6_e(x, &result));
 566 }