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

   1 /* specfunc/bessel_I0.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_bessel.h>
  26
  27 #include "error.h"
  28
  29 #include "chebyshev.h"
  30 #include "cheb_eval.c"
  31
  32 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  33
  34
  35 /* based on SLATEC besi0 */
  36
  37 /* chebyshev expansions
  38
  39  series for bi0        on the interval  0.          to  9.00000d+00
  40                                         with weighted error   2.46e-18
  41                                          log weighted error  17.61
  42                                significant figures required  17.90
  43                                     decimal places required  18.15
  44
  45  series for ai0        on the interval  1.25000d-01 to  3.33333d-01
  46                                         with weighted error   7.87e-17
  47                                          log weighted error  16.10
  48                                significant figures required  14.69
  49                                     decimal places required  16.76
  50
  51
  52  series for ai02       on the interval  0.          to  1.25000d-01
  53                                         with weighted error   3.79e-17
  54                                          log weighted error  16.42
  55                                significant figures required  14.86
  56                                     decimal places required  17.09
  57 */
  58
  59 static double bi0_data[12] = {
  60   -.07660547252839144951,
  61   1.92733795399380827000,
  62    .22826445869203013390,
  63    .01304891466707290428,
  64    .00043442709008164874,
  65    .00000942265768600193,
  66    .00000014340062895106,
  67    .00000000161384906966,
  68    .00000000001396650044,
  69    .00000000000009579451,
  70    .00000000000000053339,
  71    .00000000000000000245
  72 };
  73 static cheb_series bi0_cs = {
  74   bi0_data,
  75   11,
  76   -1, 1,
  77   11
  78 };
  79
  80 static double ai0_data[21] = {
  81    .07575994494023796,
  82    .00759138081082334,
  83    .00041531313389237,
  84    .00001070076463439,
  85   -.00000790117997921,
  86   -.00000078261435014,
  87    .00000027838499429,
  88    .00000000825247260,
  89   -.00000001204463945,
  90    .00000000155964859,
  91    .00000000022925563,
  92   -.00000000011916228,
  93    .00000000001757854,
  94    .00000000000112822,
  95   -.00000000000114684,
  96    .00000000000027155,
  97   -.00000000000002415,
  98   -.00000000000000608,
  99    .00000000000000314,
 100   -.00000000000000071,
 101    .00000000000000007
 102 };
 103 static cheb_series ai0_cs = {
 104   ai0_data,
 105   20,
 106   -1, 1,
 107   13
 108 };
 109
 110 static double ai02_data[22] = {
 111    .05449041101410882,
 112    .00336911647825569,
 113    .00006889758346918,
 114    .00000289137052082,
 115    .00000020489185893,
 116    .00000002266668991,
 117    .00000000339623203,
 118    .00000000049406022,
 119    .00000000001188914,
 120   -.00000000003149915,
 121   -.00000000001321580,
 122   -.00000000000179419,
 123    .00000000000071801,
 124    .00000000000038529,
 125    .00000000000001539,
 126   -.00000000000004151,
 127   -.00000000000000954,
 128    .00000000000000382,
 129    .00000000000000176,
 130   -.00000000000000034,
 131   -.00000000000000027,
 132    .00000000000000003
 133 };
 134 static cheb_series ai02_cs = {
 135   ai02_data,
 136   21,
 137   -1, 1,
 138   11
 139 };
 140
 141
 142 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 143
 144 int gsl_sf_bessel_I0_scaled_e(const double x, gsl_sf_result * result)
 145 {
 146   double y = fabs(x);
 147
 148   /* CHECK_POINTER(result) */
 149
 150   if(y < 2.0 * GSL_SQRT_DBL_EPSILON) {
 151     result->val = 1.0 - y;
 152     result->err = 0.5*y*y;
 153     return GSL_SUCCESS;
 154   }
 155   else if(y <= 3.0) {
 156     const double ey = exp(-y);
 157     gsl_sf_result c;
 158     cheb_eval_e(&bi0_cs, y*y/4.5-1.0, &c);
 159     result->val = ey * (2.75 + c.val);
 160     result->err = GSL_DBL_EPSILON * fabs(result->val) + ey * c.err;
 161     return GSL_SUCCESS;
 162   }
 163   else if(y <= 8.0) {
 164     const double sy = sqrt(y);
 165     gsl_sf_result c;
 166     cheb_eval_e(&ai0_cs, (48.0/y-11.0)/5.0, &c);
 167     result->val  = (0.375 + c.val) / sy;
 168     result->err  = 2.0 * GSL_DBL_EPSILON * (0.375 + fabs(c.val)) / sy;
 169     result->err += c.err / sy;
 170     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 171     return GSL_SUCCESS;
 172   }
 173   else {
 174     const double sy = sqrt(y);
 175     gsl_sf_result c;
 176     cheb_eval_e(&ai02_cs, 16.0/y-1.0, &c);
 177     result->val = (0.375 + c.val) / sy;
 178     result->err  = 2.0 * GSL_DBL_EPSILON * (0.375 + fabs(c.val)) / sy;
 179     result->err += c.err / sy;
 180     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 181     return GSL_SUCCESS;
 182   }
 183 }
 184
 185
 186 int gsl_sf_bessel_I0_e(const double x, gsl_sf_result * result)
 187 {
 188   double y = fabs(x);
 189
 190   /* CHECK_POINTER(result) */
 191
 192   if(y < 2.0 * GSL_SQRT_DBL_EPSILON) {
 193     result->val = 1.0;
 194     result->err = 0.5*y*y;
 195     return GSL_SUCCESS;
 196   }
 197   else if(y <= 3.0) {
 198     gsl_sf_result c;
 199     cheb_eval_e(&bi0_cs, y*y/4.5-1.0, &c);
 200     result->val  = 2.75 + c.val;
 201     result->err  = GSL_DBL_EPSILON * (2.75 + fabs(c.val));
 202     result->err += c.err;
 203     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 204     return GSL_SUCCESS;
 205   }
 206   else if(y < GSL_LOG_DBL_MAX - 1.0) {
 207     const double ey = exp(y);
 208     gsl_sf_result b_scaled;
 209     gsl_sf_bessel_I0_scaled_e(x, &b_scaled);
 210     result->val  = ey * b_scaled.val;
 211     result->err  = ey * b_scaled.err + y*GSL_DBL_EPSILON*fabs(result->val);
 212     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 213     return GSL_SUCCESS;
 214   }
 215   else {
 216     OVERFLOW_ERROR(result);
 217   }
 218 }
 219
 220 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 221
 222 #include "eval.h"
 223
 224 double gsl_sf_bessel_I0_scaled(const double x)
 225 {
 226   EVAL_RESULT(gsl_sf_bessel_I0_scaled_e(x, &result); )
 227 }
 228
 229 double gsl_sf_bessel_I0(const double x)
 230 {
 231   EVAL_RESULT(gsl_sf_bessel_I0_e(x, &result); )
 232 }