diff options
Diffstat (limited to 'gsl-1.9/specfunc/bessel_I0.c')
| -rw-r--r-- | gsl-1.9/specfunc/bessel_I0.c | 232 |
1 files changed, 232 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/bessel_I0.c b/gsl-1.9/specfunc/bessel_I0.c new file mode 100644 index 0000000..564b093 --- /dev/null +++ b/gsl-1.9/specfunc/bessel_I0.c @@ -0,0 +1,232 @@ +/* specfunc/bessel_I0.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Author: G. Jungman */ + +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_sf_bessel.h> + +#include "error.h" + +#include "chebyshev.h" +#include "cheb_eval.c" + +/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/ + + +/* based on SLATEC besi0 */ + +/* chebyshev expansions + + series for bi0 on the interval 0. to 9.00000d+00 + with weighted error 2.46e-18 + log weighted error 17.61 + significant figures required 17.90 + decimal places required 18.15 + + series for ai0 on the interval 1.25000d-01 to 3.33333d-01 + with weighted error 7.87e-17 + log weighted error 16.10 + significant figures required 14.69 + decimal places required 16.76 + + + series for ai02 on the interval 0. to 1.25000d-01 + with weighted error 3.79e-17 + log weighted error 16.42 + significant figures required 14.86 + decimal places required 17.09 +*/ + +static double bi0_data[12] = { + -.07660547252839144951, + 1.92733795399380827000, + .22826445869203013390, + .01304891466707290428, + .00043442709008164874, + .00000942265768600193, + .00000014340062895106, + .00000000161384906966, + .00000000001396650044, + .00000000000009579451, + .00000000000000053339, + .00000000000000000245 +}; +static cheb_series bi0_cs = { + bi0_data, + 11, + -1, 1, + 11 +}; + +static double ai0_data[21] = { + .07575994494023796, + .00759138081082334, + .00041531313389237, + .00001070076463439, + -.00000790117997921, + -.00000078261435014, + .00000027838499429, + .00000000825247260, + -.00000001204463945, + .00000000155964859, + .00000000022925563, + -.00000000011916228, + .00000000001757854, + .00000000000112822, + -.00000000000114684, + .00000000000027155, + -.00000000000002415, + -.00000000000000608, + .00000000000000314, + -.00000000000000071, + .00000000000000007 +}; +static cheb_series ai0_cs = { + ai0_data, + 20, + -1, 1, + 13 +}; + +static double ai02_data[22] = { + .05449041101410882, + .00336911647825569, + .00006889758346918, + .00000289137052082, + .00000020489185893, + .00000002266668991, + .00000000339623203, + .00000000049406022, + .00000000001188914, + -.00000000003149915, + -.00000000001321580, + -.00000000000179419, + .00000000000071801, + .00000000000038529, + .00000000000001539, + -.00000000000004151, + -.00000000000000954, + .00000000000000382, + .00000000000000176, + -.00000000000000034, + -.00000000000000027, + .00000000000000003 +}; +static cheb_series ai02_cs = { + ai02_data, + 21, + -1, 1, + 11 +}; + + +/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ + +int gsl_sf_bessel_I0_scaled_e(const double x, gsl_sf_result * result) +{ + double y = fabs(x); + + /* CHECK_POINTER(result) */ + + if(y < 2.0 * GSL_SQRT_DBL_EPSILON) { + result->val = 1.0 - y; + result->err = 0.5*y*y; + return GSL_SUCCESS; + } + else if(y <= 3.0) { + const double ey = exp(-y); + gsl_sf_result c; + cheb_eval_e(&bi0_cs, y*y/4.5-1.0, &c); + result->val = ey * (2.75 + c.val); + result->err = GSL_DBL_EPSILON * fabs(result->val) + ey * c.err; + return GSL_SUCCESS; + } + else if(y <= 8.0) { + const double sy = sqrt(y); + gsl_sf_result c; + cheb_eval_e(&ai0_cs, (48.0/y-11.0)/5.0, &c); + result->val = (0.375 + c.val) / sy; + result->err = 2.0 * GSL_DBL_EPSILON * (0.375 + fabs(c.val)) / sy; + result->err += c.err / sy; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + const double sy = sqrt(y); + gsl_sf_result c; + cheb_eval_e(&ai02_cs, 16.0/y-1.0, &c); + result->val = (0.375 + c.val) / sy; + result->err = 2.0 * GSL_DBL_EPSILON * (0.375 + fabs(c.val)) / sy; + result->err += c.err / sy; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } +} + + +int gsl_sf_bessel_I0_e(const double x, gsl_sf_result * result) +{ + double y = fabs(x); + + /* CHECK_POINTER(result) */ + + if(y < 2.0 * GSL_SQRT_DBL_EPSILON) { + result->val = 1.0; + result->err = 0.5*y*y; + return GSL_SUCCESS; + } + else if(y <= 3.0) { + gsl_sf_result c; + cheb_eval_e(&bi0_cs, y*y/4.5-1.0, &c); + result->val = 2.75 + c.val; + result->err = GSL_DBL_EPSILON * (2.75 + fabs(c.val)); + result->err += c.err; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(y < GSL_LOG_DBL_MAX - 1.0) { + const double ey = exp(y); + gsl_sf_result b_scaled; + gsl_sf_bessel_I0_scaled_e(x, &b_scaled); + result->val = ey * b_scaled.val; + result->err = ey * b_scaled.err + y*GSL_DBL_EPSILON*fabs(result->val); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + OVERFLOW_ERROR(result); + } +} + +/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ + +#include "eval.h" + +double gsl_sf_bessel_I0_scaled(const double x) +{ + EVAL_RESULT(gsl_sf_bessel_I0_scaled_e(x, &result); ) +} + +double gsl_sf_bessel_I0(const double x) +{ + EVAL_RESULT(gsl_sf_bessel_I0_e(x, &result); ) +} |