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

   1 /* specfunc/bessel_i.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_pow_int.h>
  26 #include <gsl/gsl_sf_bessel.h>
  27
  28 #include "error.h"
  29
  30 #include "bessel.h"
  31
  32
  33 /* i_{l+1}/i_l
  34  */
  35 static
  36 int
  37 bessel_il_CF1(const int l, const double x, const double threshold, double * ratio)
  38 {
  39   const int kmax = 2000;
  40   double tk   = 1.0;
  41   double sum  = 1.0;
  42   double rhok = 0.0;
  43   int k;
  44
  45   for(k=1; k<=kmax; k++) {
  46     double ak = (x/(2.0*l+1.0+2.0*k)) * (x/(2.0*l+3.0+2.0*k));
  47     rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));
  48     tk  *= rhok;
  49     sum += tk;
  50     if(fabs(tk/sum) < threshold) break;
  51   }
  52
  53   *ratio = x/(2.0*l+3.0) * sum;
  54
  55   if(k == kmax)
  56     GSL_ERROR ("error", GSL_EMAXITER);
  57   else
  58     return GSL_SUCCESS;
  59 }
  60
  61
  62 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  63
  64 int gsl_sf_bessel_i0_scaled_e(const double x, gsl_sf_result * result)
  65 {
  66   double ax = fabs(x);
  67
  68   /* CHECK_POINTER(result) */
  69
  70   if(x == 0.0) {
  71     result->val = 1.0;
  72     result->err = 0.0;
  73     return GSL_SUCCESS;
  74   }
  75   else if(ax < 0.2) {
  76     const double eax = exp(-ax);
  77     const double y = ax*ax;
  78     const double c1 = 1.0/6.0;
  79     const double c2 = 1.0/120.0;
  80     const double c3 = 1.0/5040.0;
  81     const double c4 = 1.0/362880.0;
  82     const double c5 = 1.0/39916800.0;
  83     const double sum = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))));
  84     result->val = eax * sum;
  85     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  86   }
  87   else if(ax < -0.5*GSL_LOG_DBL_EPSILON) {
  88     result->val = (1.0 - exp(-2.0*ax))/(2.0*ax);
  89     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  90   }
  91   else {
  92     result->val = 1.0/(2.0*ax);
  93     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  94   }
  95   return GSL_SUCCESS;
  96 }
  97
  98
  99 int gsl_sf_bessel_i1_scaled_e(const double x, gsl_sf_result * result)
 100 {
 101   double ax = fabs(x);
 102
 103   /* CHECK_POINTER(result) */
 104
 105   if(x == 0.0) {
 106     result->val = 0.0;
 107     result->err = 0.0;
 108     return GSL_SUCCESS;
 109   }
 110   else if(ax < 3.0*GSL_DBL_MIN) {
 111     UNDERFLOW_ERROR(result);
 112   }
 113   else if(ax < 0.25) {
 114     const double eax = exp(-ax);
 115     const double y  = x*x;
 116     const double c1 = 1.0/10.0;
 117     const double c2 = 1.0/280.0;
 118     const double c3 = 1.0/15120.0;
 119     const double c4 = 1.0/1330560.0;
 120     const double c5 = 1.0/172972800.0;
 121     const double sum = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))));
 122     result->val = eax * x/3.0 * sum;
 123     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 124     return GSL_SUCCESS;
 125   }
 126   else {
 127     double ex = exp(-2.0*ax);
 128     result->val = 0.5 * (ax*(1.0+ex) - (1.0-ex)) / (ax*ax);
 129     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 130     if(x < 0.0) result->val = -result->val;
 131     return GSL_SUCCESS;
 132   }
 133 }
 134
 135
 136 int gsl_sf_bessel_i2_scaled_e(const double x, gsl_sf_result * result)
 137 {
 138   double ax = fabs(x);
 139
 140   /* CHECK_POINTER(result) */
 141
 142   if(x == 0.0) {
 143     result->val = 0.0;
 144     result->err = 0.0;
 145     return GSL_SUCCESS;
 146   }
 147   else if(ax < 4.0*GSL_SQRT_DBL_MIN) {
 148     UNDERFLOW_ERROR(result);
 149   }
 150   else if(ax < 0.25) {
 151     const double y = x*x;
 152     const double c1 = 1.0/14.0;
 153     const double c2 = 1.0/504.0;
 154     const double c3 = 1.0/33264.0;
 155     const double c4 = 1.0/3459456.0;
 156     const double c5 = 1.0/518918400.0;
 157     const double sum = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))));
 158     const double pre = exp(-ax) * x*x/15.0;
 159     result->val = pre * sum;
 160     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 161     return GSL_SUCCESS;
 162   }
 163   else {
 164     double ex = exp(-2.0*ax);
 165     double x2 = x*x;
 166     result->val = 0.5 * ((3.0+x2)*(1.0-ex) - 3.0*ax*(1.0+ex))/(ax*ax*ax);
 167     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 168     return GSL_SUCCESS;
 169   }
 170 }
 171
 172
 173 int gsl_sf_bessel_il_scaled_e(const int l, double x, gsl_sf_result * result)
 174 {
 175   double sgn = 1.0;
 176   double ax  = fabs(x);
 177
 178   if(x < 0.0) {
 179     /* i_l(-x) = (-1)^l i_l(x) */
 180     sgn = ( GSL_IS_ODD(l) ? -1.0 : 1.0 );
 181     x = -x;
 182   }
 183
 184   if(l < 0) {
 185     DOMAIN_ERROR(result);
 186   }
 187   else if(x == 0.0) {
 188     result->val = ( l == 0 ? 1.0 : 0.0 );
 189     result->err = 0.0;
 190     return GSL_SUCCESS;
 191   }
 192   else if(l == 0) {
 193     gsl_sf_result il;
 194     int stat_il = gsl_sf_bessel_i0_scaled_e(x, &il);
 195     result->val = sgn * il.val;
 196     result->err = il.err;
 197     return stat_il;
 198   }
 199   else if(l == 1) {
 200     gsl_sf_result il;
 201     int stat_il = gsl_sf_bessel_i1_scaled_e(x, &il);
 202     result->val = sgn * il.val;
 203     result->err = il.err;
 204     return stat_il;
 205   }
 206   else if(l == 2) {
 207     gsl_sf_result il;
 208     int stat_il = gsl_sf_bessel_i2_scaled_e(x, &il);
 209     result->val = sgn * il.val;
 210     result->err = il.err;
 211     return stat_il;
 212   }
 213   else if(x*x < 10.0*(l+1.5)/M_E) {
 214     gsl_sf_result b;
 215     int stat = gsl_sf_bessel_IJ_taylor_e(l+0.5, x, 1, 50, GSL_DBL_EPSILON, &b);
 216     double pre   = exp(-ax) * sqrt((0.5*M_PI)/x);
 217     result->val  = sgn * pre * b.val;
 218     result->err  = pre * b.err;
 219     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 220     return stat;
 221   }
 222   else if(l < 150) {
 223     gsl_sf_result i0_scaled;
 224     int stat_i0  = gsl_sf_bessel_i0_scaled_e(ax, &i0_scaled);
 225     double rat;
 226     int stat_CF1 = bessel_il_CF1(l, ax, GSL_DBL_EPSILON, &rat);
 227     double iellp1 = rat * GSL_SQRT_DBL_MIN;
 228     double iell   = GSL_SQRT_DBL_MIN;
 229     double iellm1;
 230     int ell;
 231     for(ell = l; ell >= 1; ell--) {
 232       iellm1 = iellp1 + (2*ell + 1)/x * iell;
 233       iellp1 = iell;
 234       iell   = iellm1;
 235     }
 236     result->val  = sgn * i0_scaled.val * (GSL_SQRT_DBL_MIN / iell);
 237     result->err  = i0_scaled.err * (GSL_SQRT_DBL_MIN / iell);
 238     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 239     return GSL_ERROR_SELECT_2(stat_i0, stat_CF1);
 240   }
 241   else if(GSL_MIN(0.29/(l*l+1.0), 0.5/(l*l+1.0+x*x)) < 0.5*GSL_ROOT3_DBL_EPSILON) {
 242     int status = gsl_sf_bessel_Inu_scaled_asymp_unif_e(l + 0.5, x, result);
 243     double pre = sqrt((0.5*M_PI)/x);
 244     result->val *= sgn * pre;
 245     result->err *= pre;
 246     return status;
 247   }
 248   else {
 249     /* recurse down from safe values */
 250     double rt_term = sqrt((0.5*M_PI)/x);
 251     const int LMAX = 2 + (int) (1.2 / GSL_ROOT6_DBL_EPSILON);
 252     gsl_sf_result r_iellp1;
 253     gsl_sf_result r_iell;
 254     int stat_a1 = gsl_sf_bessel_Inu_scaled_asymp_unif_e(LMAX + 1 + 0.5, x, &r_iellp1);
 255     int stat_a2 = gsl_sf_bessel_Inu_scaled_asymp_unif_e(LMAX     + 0.5, x, &r_iell);
 256     double iellp1 = r_iellp1.val;
 257     double iell   = r_iell.val;
 258     double iellm1 = 0.0;
 259     int ell;
 260     iellp1 *= rt_term;
 261     iell   *= rt_term;
 262     for(ell = LMAX; ell >= l+1; ell--) {
 263       iellm1 = iellp1 + (2*ell + 1)/x * iell;
 264       iellp1 = iell;
 265       iell   = iellm1;
 266     }
 267     result->val  = sgn * iellm1;
 268     result->err  = fabs(result->val)*(GSL_DBL_EPSILON + fabs(r_iellp1.err/r_iellp1.val) + fabs(r_iell.err/r_iell.val));
 269     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 270
 271     return GSL_ERROR_SELECT_2(stat_a1, stat_a2);
 272   }
 273 }
 274
 275
 276 int gsl_sf_bessel_il_scaled_array(const int lmax, const double x, double * result_array)
 277 {
 278   if(x == 0.0) {
 279     int ell;
 280     result_array[0] = 1.0;
 281     for (ell = lmax; ell >= 1; ell--) {
 282       result_array[ell] = 0.0;
 283     };
 284     return GSL_SUCCESS;
 285   } else {
 286     int ell;
 287     gsl_sf_result r_iellp1;
 288     gsl_sf_result r_iell;
 289     int stat_0 = gsl_sf_bessel_il_scaled_e(lmax+1, x, &r_iellp1);
 290     int stat_1 = gsl_sf_bessel_il_scaled_e(lmax,   x, &r_iell);
 291     double iellp1 = r_iellp1.val;
 292     double iell   = r_iell.val;
 293     double iellm1;
 294     result_array[lmax] = iell;
 295     for(ell = lmax; ell >= 1; ell--) {
 296       iellm1 = iellp1 + (2*ell + 1)/x * iell;
 297       iellp1 = iell;
 298       iell   = iellm1;
 299       result_array[ell-1] = iellm1;
 300     }
 301     return GSL_ERROR_SELECT_2(stat_0, stat_1);
 302   }
 303 }
 304
 305
 306 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 307
 308 #include "eval.h"
 309
 310 double gsl_sf_bessel_i0_scaled(const double x)
 311 {
 312   EVAL_RESULT(gsl_sf_bessel_i0_scaled_e(x, &result));
 313 }
 314
 315 double gsl_sf_bessel_i1_scaled(const double x)
 316 {
 317   EVAL_RESULT(gsl_sf_bessel_i1_scaled_e(x, &result));
 318 }
 319
 320 double gsl_sf_bessel_i2_scaled(const double x)
 321 {
 322   EVAL_RESULT(gsl_sf_bessel_i2_scaled_e(x, &result));
 323 }
 324
 325 double gsl_sf_bessel_il_scaled(const int l, const double x)
 326 {
 327   EVAL_RESULT(gsl_sf_bessel_il_scaled_e(l, x, &result));
 328 }
 329