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

   1 /* specfunc/transport.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_transport.h>
  26
  27 #include "error.h"
  28 #include "check.h"
  29
  30 #include "chebyshev.h"
  31 #include "cheb_eval.c"
  32
  33 static double transport2_data[18] = {
  34    1.671760446434538503,
  35   -0.147735359946794490,
  36    0.148213819946936338e-01,
  37   -0.14195330326305613e-02,
  38    0.1306541324415708e-03,
  39   -0.117155795867579e-04,
  40    0.10333498445756e-05,
  41   -0.901911304223e-07,
  42    0.78177169833e-08,
  43   -0.6744565684e-09,
  44    0.579946394e-10,
  45   -0.49747619e-11,
  46    0.425961e-12,
  47   -0.36422e-13,
  48    0.3111e-14,
  49   -0.265e-15,
  50    0.23e-16,
  51   -0.19e-17
  52 };
  53 static cheb_series transport2_cs = {
  54   transport2_data,
  55   17,
  56   -1, 1,
  57   9
  58 };
  59
  60 static double transport3_data[18] = {
  61    0.762012543243872007,
  62   -0.105674387705058533,
  63    0.119778084819657810e-01,
  64   -0.12144015203698307e-02,
  65    0.1155099769392855e-03,
  66   -0.105815992124423e-04,
  67    0.9474663385302e-06,
  68   -0.836221212858e-07,
  69    0.73109099278e-08,
  70   -0.6350594779e-09,
  71    0.549118282e-10,
  72   -0.47321395e-11,
  73    0.4067695e-12,
  74   -0.348971e-13,
  75    0.29892e-14,
  76   -0.256e-15,
  77    0.219e-16,
  78   -0.19e-17
  79 };
  80 static cheb_series transport3_cs = {
  81   transport3_data,
  82   17,
  83   -1, 1,
  84   9
  85 };
  86
  87
  88 static double transport4_data[18] = {
  89   0.4807570994615110579,
  90  -0.8175378810321083956e-01,
  91   0.1002700665975162973e-01,
  92  -0.10599339359820151e-02,
  93   0.1034506245030405e-03,
  94  -0.96442705485899e-05,
  95   0.8745544408515e-06,
  96  -0.779321207981e-07,
  97   0.68649886141e-08,
  98  -0.5999571076e-09,
  99   0.521366241e-10,
 100  -0.45118382e-11,
 101   0.3892159e-12,
 102  -0.334936e-13,
 103   0.28767e-14,
 104  -0.2467e-15,
 105   0.211e-16,
 106  -0.18e-17
 107 };
 108 static cheb_series transport4_cs = {
 109   transport4_data,
 110   17,
 111   -1, 1,
 112   9
 113 };
 114
 115
 116 static double transport5_data[18] = {
 117    0.347777777133910789,
 118   -0.66456988976050428e-01,
 119    0.8611072656883309e-02,
 120   -0.9396682223755538e-03,
 121    0.936324806081513e-04,
 122   -0.88571319340833e-05,
 123    0.811914989145e-06,
 124   -0.72957654233e-07,
 125    0.646971455e-08,
 126   -0.568490283e-09,
 127    0.49625598e-10,
 128   -0.4310940e-11,
 129    0.373100e-12,
 130   -0.32198e-13,
 131    0.2772e-14,
 132   -0.238e-15,
 133    0.21e-16,
 134   -0.18e-17
 135 };
 136 static cheb_series transport5_cs = {
 137   transport5_data,
 138   17,
 139   -1, 1,
 140   9
 141 };
 142
 143
 144 static
 145 double
 146 transport_sumexp(const int numexp, const int order, const double t, double x)
 147 {
 148   double rk = (double)numexp;
 149   double sumexp = 0.0;
 150   int k;
 151   for(k=1; k<=numexp; k++) {
 152     double sum2 = 1.0;
 153     double xk  = 1.0/(rk*x);
 154     double xk1 = 1.0;
 155     int j;
 156     for(j=1; j<=order; j++) {
 157       sum2 = sum2*xk1*xk + 1.0;
 158       xk1 += 1.0;
 159     }
 160     sumexp *= t;
 161     sumexp += sum2;
 162     rk -= 1.0;
 163   }
 164   return sumexp;
 165 }
 166
 167
 168 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 169
 170 int
 171 gsl_sf_transport_2_e(const double x, gsl_sf_result * result)
 172 {
 173   const double val_infinity = 3.289868133696452873;
 174
 175   /* CHECK_POINTER(result) */
 176
 177   if(x < 0.0) {
 178     DOMAIN_ERROR(result);
 179   }
 180   else if(x < 3.0*GSL_SQRT_DBL_EPSILON) {
 181     result->val = x;
 182     result->err = GSL_DBL_EPSILON*fabs(x) + x*x/2.0;
 183     return GSL_SUCCESS;
 184   }
 185   else if(x <= 4.0) {
 186     double t = (x*x/8.0 - 0.5) - 0.5;
 187     gsl_sf_result result_c;
 188     cheb_eval_e(&transport2_cs, t, &result_c);
 189     result->val  = x * result_c.val;
 190     result->err  = x * result_c.err;
 191     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 192     return GSL_SUCCESS;
 193   }
 194   else if(x < -GSL_LOG_DBL_EPSILON) {
 195     const int    numexp = (int)((-GSL_LOG_DBL_EPSILON)/x) + 1;
 196     const double sumexp = transport_sumexp(numexp, 2, exp(-x), x);
 197     const double t = 2.0 * log(x) - x + log(sumexp);
 198     if(t < GSL_LOG_DBL_EPSILON) {
 199       result->val  = val_infinity;
 200       result->err  = 2.0 * GSL_DBL_EPSILON * val_infinity;
 201     }
 202     else {
 203       const double et = exp(t);
 204       result->val = val_infinity - et;
 205       result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + fabs(t) * et);
 206     }
 207     return GSL_SUCCESS;
 208   }
 209   else if(x < 2.0/GSL_DBL_EPSILON) {
 210     const int    numexp = 1;
 211     const double sumexp = transport_sumexp(numexp, 2, 1.0, x);
 212     const double t = 2.0 * log(x) - x + log(sumexp);
 213     if(t < GSL_LOG_DBL_EPSILON) {
 214       result->val = val_infinity;
 215       result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
 216     }
 217     else {
 218       const double et = exp(t);
 219       result->val = val_infinity - et;
 220       result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
 221     }
 222     return GSL_SUCCESS;
 223   }
 224   else {
 225     const double t = 2.0 * log(x) - x;
 226     if(t < GSL_LOG_DBL_EPSILON) {
 227       result->val = val_infinity;
 228       result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
 229     }
 230     else {
 231       const double et = exp(t);
 232       result->val = val_infinity - et;
 233       result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
 234     }
 235     return GSL_SUCCESS;
 236   }
 237 }
 238
 239
 240 int
 241 gsl_sf_transport_3_e(const double x, gsl_sf_result * result)
 242 {
 243   const double val_infinity = 7.212341418957565712;
 244
 245   /* CHECK_POINTER(result) */
 246
 247   if(x < 0.0) {
 248     DOMAIN_ERROR(result);
 249   }
 250   else if(x == 0.0) {
 251     result->val = 0.0;
 252     result->err = 0.0;
 253     return GSL_SUCCESS;
 254   }
 255   else if(x < 3.0*GSL_SQRT_DBL_EPSILON) {
 256     result->val = 0.5*x*x;
 257     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
 258     CHECK_UNDERFLOW(result);
 259     return GSL_SUCCESS;
 260   }
 261   else if(x <= 4.0) {
 262     const double x2 = x*x;
 263     const double t = (x2/8.0 - 0.5) - 0.5;
 264     gsl_sf_result result_c;
 265     cheb_eval_e(&transport3_cs, t, &result_c);
 266     result->val  = x2 * result_c.val;
 267     result->err  = x2 * result_c.err;
 268     result->err += GSL_DBL_EPSILON * fabs(result->val);
 269     return GSL_SUCCESS;
 270   }
 271   else if(x < -GSL_LOG_DBL_EPSILON) {
 272     const int    numexp = (int)((-GSL_LOG_DBL_EPSILON)/x) + 1;
 273     const double sumexp = transport_sumexp(numexp, 3, exp(-x), x);
 274     const double t = 3.0 * log(x) - x + log(sumexp);
 275     if(t < GSL_LOG_DBL_EPSILON) {
 276       result->val = val_infinity;
 277       result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
 278     }
 279     else {
 280       const double et = exp(t);
 281       result->val = val_infinity - et;
 282       result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + fabs(t) * et);
 283     }
 284     return GSL_SUCCESS;
 285   }
 286   else if(x < 3.0/GSL_DBL_EPSILON) {
 287     const int    numexp = 1;
 288     const double sumexp = transport_sumexp(numexp, 3, 1.0, x);
 289     const double t = 3.0 * log(x) - x + log(sumexp);
 290     if(t < GSL_LOG_DBL_EPSILON) {
 291       result->val = val_infinity;
 292       result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
 293     }
 294     else {
 295       const double et = exp(t);
 296       result->val = val_infinity - et;
 297       result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
 298     }
 299     return GSL_SUCCESS;
 300   }
 301   else {
 302     const double t = 3.0 * log(x) - x;
 303     if(t < GSL_LOG_DBL_EPSILON) {
 304       result->val = val_infinity;
 305       result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
 306     }
 307     else {
 308       const double et = exp(t);
 309       result->val = val_infinity - et;
 310       result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
 311     }
 312     return GSL_SUCCESS;
 313   }
 314 }
 315
 316
 317 int
 318 gsl_sf_transport_4_e(const double x, gsl_sf_result * result)
 319 {
 320   const double val_infinity = 25.97575760906731660;
 321
 322   /* CHECK_POINTER(result) */
 323
 324   if(x < 0.0) {
 325     DOMAIN_ERROR(result);
 326   }
 327   else if(x == 0.0) {
 328     result->val = 0.0;
 329     result->err = 0.0;
 330     return GSL_SUCCESS;
 331   }
 332   else if(x < 3.0*GSL_SQRT_DBL_EPSILON) {
 333     result->val = x*x*x/3.0;
 334     result->err = 3.0 * GSL_DBL_EPSILON * result->val;
 335     CHECK_UNDERFLOW(result);
 336     return GSL_SUCCESS;
 337   }
 338   else if(x <= 4.0) {
 339     const double x2 = x*x;
 340     const double t = (x2/8.0 - 0.5) - 0.5;
 341     gsl_sf_result result_c;
 342     cheb_eval_e(&transport4_cs, t, &result_c);
 343     result->val  = x2*x * result_c.val;
 344     result->err  = x2*x * result_c.err;
 345     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 346     return GSL_SUCCESS;
 347   }
 348   else if(x < -GSL_LOG_DBL_EPSILON) {
 349     const int    numexp = (int)((-GSL_LOG_DBL_EPSILON)/x) + 1;
 350     const double sumexp = transport_sumexp(numexp, 4, exp(-x), x);
 351     const double t = 4.0 * log(x) - x + log(sumexp);
 352     if(t < GSL_LOG_DBL_EPSILON) {
 353       result->val = val_infinity;
 354       result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
 355     }
 356     else {
 357       const double et = exp(t);
 358       result->val = val_infinity - et;
 359       result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
 360     }
 361     return GSL_SUCCESS;
 362   }
 363   else if(x < 3.0/GSL_DBL_EPSILON) {
 364     const int    numexp = 1;
 365     const double sumexp = transport_sumexp(numexp, 4, 1.0, x);
 366     const double t = 4.0 * log(x) - x + log(sumexp);
 367     if(t < GSL_LOG_DBL_EPSILON) {
 368       result->val = val_infinity;
 369       result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
 370     }
 371     else {
 372       const double et = exp(t);
 373       result->val = val_infinity - et;
 374       result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
 375     }
 376     return GSL_SUCCESS;
 377   }
 378   else {
 379     const double t = 4.0 * log(x) - x;
 380     if(t < GSL_LOG_DBL_EPSILON) {
 381       result->val = val_infinity;
 382       result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
 383     }
 384     else {
 385       const double et = exp(t);
 386       result->val = val_infinity - et;
 387       result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
 388     }
 389     return GSL_SUCCESS;
 390   }
 391 }
 392
 393
 394 int
 395 gsl_sf_transport_5_e(const double x, gsl_sf_result * result)
 396 {
 397   const double val_infinity = 124.4313306172043912;
 398
 399   /* CHECK_POINTER(result) */
 400
 401   if(x < 0.0) {
 402     DOMAIN_ERROR(result);
 403   }
 404   else if(x == 0.0) {
 405     result->val = 0.0;
 406     result->err = 0.0;
 407     return GSL_SUCCESS;
 408   }
 409   else if(x < 3.0*GSL_SQRT_DBL_EPSILON) {
 410     result->val = x*x*x*x/4.0;
 411     result->err = 4.0 * GSL_DBL_EPSILON * result->val;
 412     CHECK_UNDERFLOW(result);
 413     return GSL_SUCCESS;
 414   }
 415   else if(x <= 4.0) {
 416     const double x2 = x*x;
 417     const double t = (x2/8.0 - 0.5) - 0.5;
 418     gsl_sf_result result_c;
 419     cheb_eval_e(&transport5_cs, t, &result_c);
 420     result->val  = x2*x2 * result_c.val;
 421     result->err  = x2*x2 * result_c.err;
 422     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 423     return GSL_SUCCESS;
 424   }
 425   else if(x < -GSL_LOG_DBL_EPSILON) {
 426     const int    numexp = (int)((-GSL_LOG_DBL_EPSILON)/x) + 1;
 427     const double sumexp = transport_sumexp(numexp, 5, exp(-x), x);
 428     const double t = 5.0 * log(x) - x + log(sumexp);
 429     if(t < GSL_LOG_DBL_EPSILON) {
 430       result->val = val_infinity;
 431       result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
 432     }
 433     else {
 434       const double et = exp(t);
 435       result->val = val_infinity - et;
 436       result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
 437     }
 438     return GSL_SUCCESS;
 439   }
 440   else if(x < 3.0/GSL_DBL_EPSILON) {
 441     const int    numexp = 1;
 442     const double sumexp = transport_sumexp(numexp, 5, 1.0, x);
 443     const double t = 5.0 * log(x) - x + log(sumexp);
 444     if(t < GSL_LOG_DBL_EPSILON) {
 445       result->val = val_infinity;
 446       result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
 447     }
 448     else {
 449       const double et = exp(t);
 450       result->val = val_infinity - et;
 451       result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
 452     }
 453     return GSL_SUCCESS;
 454   }
 455   else {
 456     const double t = 5.0 * log(x) - x;
 457     if(t < GSL_LOG_DBL_EPSILON) {
 458       result->val = val_infinity;
 459       result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
 460     }
 461     else {
 462       const double et = exp(t);
 463       result->val = val_infinity - et;
 464       result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
 465     }
 466     return GSL_SUCCESS;
 467   }
 468 }
 469
 470 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 471
 472 #include "eval.h"
 473
 474 double gsl_sf_transport_2(const double x)
 475 {
 476   EVAL_RESULT(gsl_sf_transport_2_e(x, &result));
 477 }
 478
 479 double gsl_sf_transport_3(const double x)
 480 {
 481   EVAL_RESULT(gsl_sf_transport_3_e(x, &result));
 482 }
 483
 484 double gsl_sf_transport_4(const double x)
 485 {
 486   EVAL_RESULT(gsl_sf_transport_4_e(x, &result));
 487 }
 488
 489 double gsl_sf_transport_5(const double x)
 490 {
 491   EVAL_RESULT(gsl_sf_transport_5_e(x, &result));
 492 }