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

   1 /* specfunc/bessel_In.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 "bessel.h"
  30
  31 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  32
  33
  34 int
  35 gsl_sf_bessel_In_scaled_e(int n, const double x, gsl_sf_result * result)
  36 {
  37   const double ax = fabs(x);
  38
  39   n = abs(n);  /* I(-n, z) = I(n, z) */
  40
  41   /* CHECK_POINTER(result) */
  42
  43   if(n == 0) {
  44     return gsl_sf_bessel_I0_scaled_e(x, result);
  45   }
  46   else if(n == 1) {
  47     return gsl_sf_bessel_I1_scaled_e(x, result);
  48   }
  49   else if(x == 0.0) {
  50     result->val = 0.0;
  51     result->err = 0.0;
  52     return GSL_SUCCESS;
  53   }
  54   else if(x*x < 10.0*(n+1.0)/M_E) {
  55     gsl_sf_result t;
  56     double ex   = exp(-ax);
  57     int stat_In = gsl_sf_bessel_IJ_taylor_e((double)n, ax, 1, 50, GSL_DBL_EPSILON, &t);
  58     result->val  = t.val * ex;
  59     result->err  = t.err * ex;
  60     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  61     if(x < 0.0 && GSL_IS_ODD(n)) result->val = -result->val;
  62     return stat_In;
  63   }
  64   else if(n < 150 && ax < 1e7) {
  65     gsl_sf_result I0_scaled;
  66     int stat_I0 = gsl_sf_bessel_I0_scaled_e(ax, &I0_scaled);
  67     double rat;
  68     int stat_CF1 = gsl_sf_bessel_I_CF1_ser((double)n, ax, &rat);
  69     double Ikp1 = rat * GSL_SQRT_DBL_MIN;
  70     double Ik   = GSL_SQRT_DBL_MIN;
  71     double Ikm1;
  72     int k;
  73     for(k=n; k >= 1; k--) {
  74       Ikm1 = Ikp1 + 2.0*k/ax * Ik;
  75       Ikp1 = Ik;
  76       Ik   = Ikm1;
  77     }
  78     result->val  = I0_scaled.val * (GSL_SQRT_DBL_MIN / Ik);
  79     result->err  = I0_scaled.err * (GSL_SQRT_DBL_MIN / Ik);
  80     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  81     if(x < 0.0 && GSL_IS_ODD(n)) result->val = -result->val;
  82     return GSL_ERROR_SELECT_2(stat_I0, stat_CF1);
  83   }
  84   else if( GSL_MIN( 0.29/(n*n), 0.5/(n*n + x*x) ) < 0.5*GSL_ROOT3_DBL_EPSILON) {
  85     int stat_as = gsl_sf_bessel_Inu_scaled_asymp_unif_e((double)n, ax, result);
  86     if(x < 0.0 && GSL_IS_ODD(n)) result->val = -result->val;
  87     return stat_as;
  88   }
  89   else {
  90     const int nhi = 2 + (int) (1.2 / GSL_ROOT6_DBL_EPSILON);
  91     gsl_sf_result r_Ikp1;
  92     gsl_sf_result r_Ik;
  93     int stat_a1 = gsl_sf_bessel_Inu_scaled_asymp_unif_e(nhi+1.0,     ax, &r_Ikp1);
  94     int stat_a2 = gsl_sf_bessel_Inu_scaled_asymp_unif_e((double)nhi, ax, &r_Ik);
  95     double Ikp1 = r_Ikp1.val;
  96     double Ik   = r_Ik.val;
  97     double Ikm1;
  98     int k;
  99     for(k=nhi; k > n; k--) {
 100       Ikm1 = Ikp1 + 2.0*k/ax * Ik;
 101       Ikp1 = Ik;
 102       Ik   = Ikm1;
 103     }
 104     result->val = Ik;
 105     result->err = Ik * (r_Ikp1.err/r_Ikp1.val + r_Ik.err/r_Ik.val);
 106     if(x < 0.0 && GSL_IS_ODD(n)) result->val = -result->val;
 107     return GSL_ERROR_SELECT_2(stat_a1, stat_a2);
 108   }
 109 }
 110
 111
 112 int
 113 gsl_sf_bessel_In_scaled_array(const int nmin, const int nmax, const double x, double * result_array)
 114 {
 115   /* CHECK_POINTER(result_array) */
 116
 117   if(nmax < nmin || nmin < 0) {
 118     int j;
 119     for(j=0; j<=nmax-nmin; j++) result_array[j] = 0.0;
 120     GSL_ERROR ("domain error", GSL_EDOM);
 121   }
 122   else if(x == 0.0) {
 123     int j;
 124     for(j=0; j<=nmax-nmin; j++) result_array[j] = 0.0;
 125     if(nmin == 0) result_array[0] = 1.0;
 126     return GSL_SUCCESS;
 127   }
 128   else if(nmax == 0) {
 129     gsl_sf_result I0_scaled;
 130     int stat = gsl_sf_bessel_I0_scaled_e(x, &I0_scaled);
 131     result_array[0] = I0_scaled.val;
 132     return stat;
 133   }
 134   else {
 135     const double ax = fabs(x);
 136     const double two_over_x = 2.0/ax;
 137
 138     /* starting values */
 139     gsl_sf_result r_Inp1;
 140     gsl_sf_result r_In;
 141     int stat_0 = gsl_sf_bessel_In_scaled_e(nmax+1, ax, &r_Inp1);
 142     int stat_1 = gsl_sf_bessel_In_scaled_e(nmax,   ax, &r_In);
 143     double Inp1 = r_Inp1.val;
 144     double In   = r_In.val;
 145     double Inm1;
 146     int n;
 147
 148     for(n=nmax; n>=nmin; n--) {
 149       result_array[n-nmin] = In;
 150       Inm1 = Inp1 + n * two_over_x * In;
 151       Inp1 = In;
 152       In   = Inm1;
 153     }
 154
 155     /* deal with signs */
 156     if(x < 0.0) {
 157       for(n=nmin; n<=nmax; n++) {
 158         if(GSL_IS_ODD(n)) result_array[n-nmin] = -result_array[n-nmin];
 159       }
 160     }
 161
 162     return GSL_ERROR_SELECT_2(stat_0, stat_1);
 163   }
 164 }
 165
 166
 167 int
 168 gsl_sf_bessel_In_e(const int n_in, const double x, gsl_sf_result * result)
 169 {
 170   const double ax = fabs(x);
 171   const int n = abs(n_in);  /* I(-n, z) = I(n, z) */
 172   gsl_sf_result In_scaled;
 173   const int stat_In_scaled = gsl_sf_bessel_In_scaled_e(n, ax, &In_scaled);
 174
 175   /* In_scaled is always less than 1,
 176    * so this overflow check is conservative.
 177    */
 178   if(ax > GSL_LOG_DBL_MAX - 1.0) {
 179     OVERFLOW_ERROR(result);
 180   }
 181   else {
 182     const double ex = exp(ax);
 183     result->val  = ex * In_scaled.val;
 184     result->err  = ex * In_scaled.err;
 185     result->err += ax * GSL_DBL_EPSILON * fabs(result->val);
 186     if(x < 0.0 && GSL_IS_ODD(n)) result->val = -result->val;
 187     return stat_In_scaled;
 188   }
 189 }
 190
 191
 192 int
 193 gsl_sf_bessel_In_array(const int nmin, const int nmax, const double x, double * result_array)
 194 {
 195   double ax = fabs(x);
 196
 197   /* CHECK_POINTER(result_array) */
 198
 199   if(ax > GSL_LOG_DBL_MAX - 1.0) {
 200     int j;
 201     for(j=0; j<=nmax-nmin; j++) result_array[j] = 0.0; /* FIXME: should be Inf */
 202     GSL_ERROR ("overflow", GSL_EOVRFLW);
 203   }
 204   else {
 205     int j;
 206     double eax = exp(ax);
 207     int status = gsl_sf_bessel_In_scaled_array(nmin, nmax, x, result_array);
 208     for(j=0; j<=nmax-nmin; j++) result_array[j] *= eax;
 209     return status;
 210   }
 211 }
 212
 213 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 214
 215 #include "eval.h"
 216
 217 double gsl_sf_bessel_In_scaled(const int n, const double x)
 218 {
 219   EVAL_RESULT(gsl_sf_bessel_In_scaled_e(n, x, &result));
 220 }
 221
 222 double gsl_sf_bessel_In(const int n, const double x)
 223 {
 224   EVAL_RESULT(gsl_sf_bessel_In_e(n, x, &result));
 225 }