diff options
Diffstat (limited to 'gsl-1.9/specfunc/bessel_I1.c')
-rw-r--r-- | gsl-1.9/specfunc/bessel_I1.c | 258 |
1 files changed, 258 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/bessel_I1.c b/gsl-1.9/specfunc/bessel_I1.c new file mode 100644 index 0000000..18b2c0e --- /dev/null +++ b/gsl-1.9/specfunc/bessel_I1.c @@ -0,0 +1,258 @@ +/* specfunc/bessel_I1.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" + +#define ROOT_EIGHT (2.0*M_SQRT2) + + +/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/ + +/* based on SLATEC besi1(), besi1e() */ + +/* chebyshev expansions + + series for bi1 on the interval 0. to 9.00000d+00 + with weighted error 2.40e-17 + log weighted error 16.62 + significant figures required 16.23 + decimal places required 17.14 + + series for ai1 on the interval 1.25000d-01 to 3.33333d-01 + with weighted error 6.98e-17 + log weighted error 16.16 + significant figures required 14.53 + decimal places required 16.82 + + series for ai12 on the interval 0. to 1.25000d-01 + with weighted error 3.55e-17 + log weighted error 16.45 + significant figures required 14.69 + decimal places required 17.12 +*/ + +static double bi1_data[11] = { + -0.001971713261099859, + 0.407348876675464810, + 0.034838994299959456, + 0.001545394556300123, + 0.000041888521098377, + 0.000000764902676483, + 0.000000010042493924, + 0.000000000099322077, + 0.000000000000766380, + 0.000000000000004741, + 0.000000000000000024 +}; +static cheb_series bi1_cs = { + bi1_data, + 10, + -1, 1, + 10 +}; + +static double ai1_data[21] = { + -0.02846744181881479, + -0.01922953231443221, + -0.00061151858579437, + -0.00002069971253350, + 0.00000858561914581, + 0.00000104949824671, + -0.00000029183389184, + -0.00000001559378146, + 0.00000001318012367, + -0.00000000144842341, + -0.00000000029085122, + 0.00000000012663889, + -0.00000000001664947, + -0.00000000000166665, + 0.00000000000124260, + -0.00000000000027315, + 0.00000000000002023, + 0.00000000000000730, + -0.00000000000000333, + 0.00000000000000071, + -0.00000000000000006 +}; +static cheb_series ai1_cs = { + ai1_data, + 20, + -1, 1, + 11 +}; + +static double ai12_data[22] = { + 0.02857623501828014, + -0.00976109749136147, + -0.00011058893876263, + -0.00000388256480887, + -0.00000025122362377, + -0.00000002631468847, + -0.00000000383538039, + -0.00000000055897433, + -0.00000000001897495, + 0.00000000003252602, + 0.00000000001412580, + 0.00000000000203564, + -0.00000000000071985, + -0.00000000000040836, + -0.00000000000002101, + 0.00000000000004273, + 0.00000000000001041, + -0.00000000000000382, + -0.00000000000000186, + 0.00000000000000033, + 0.00000000000000028, + -0.00000000000000003 +}; +static cheb_series ai12_cs = { + ai12_data, + 21, + -1, 1, + 9 +}; + + +/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ + +int gsl_sf_bessel_I1_scaled_e(const double x, gsl_sf_result * result) +{ + const double xmin = 2.0 * GSL_DBL_MIN; + const double x_small = ROOT_EIGHT * GSL_SQRT_DBL_EPSILON; + const double y = fabs(x); + + /* CHECK_POINTER(result) */ + + if(y == 0.0) { + result->val = 0.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(y < xmin) { + UNDERFLOW_ERROR(result); + } + else if(y < x_small) { + result->val = 0.5*x; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(y <= 3.0) { + const double ey = exp(-y); + gsl_sf_result c; + cheb_eval_e(&bi1_cs, y*y/4.5-1.0, &c); + result->val = x * ey * (0.875 + c.val); + result->err = ey * c.err + y * GSL_DBL_EPSILON * fabs(result->val); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(y <= 8.0) { + const double sy = sqrt(y); + gsl_sf_result c; + double b; + double s; + cheb_eval_e(&ai1_cs, (48.0/y-11.0)/5.0, &c); + b = (0.375 + c.val) / sy; + s = (x > 0.0 ? 1.0 : -1.0); + result->val = s * b; + 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; + double b; + double s; + cheb_eval_e(&ai12_cs, 16.0/y-1.0, &c); + b = (0.375 + c.val) / sy; + s = (x > 0.0 ? 1.0 : -1.0); + result->val = s * b; + result->err = c.err / sy; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } +} + + +int gsl_sf_bessel_I1_e(const double x, gsl_sf_result * result) +{ + const double xmin = 2.0 * GSL_DBL_MIN; + const double x_small = ROOT_EIGHT * GSL_SQRT_DBL_EPSILON; + const double y = fabs(x); + + /* CHECK_POINTER(result) */ + + if(y == 0.0) { + result->val = 0.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(y < xmin) { + UNDERFLOW_ERROR(result); + } + else if(y < x_small) { + result->val = 0.5*x; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(y <= 3.0) { + gsl_sf_result c; + cheb_eval_e(&bi1_cs, y*y/4.5-1.0, &c); + result->val = x * (0.875 + c.val); + result->err = y * c.err; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(y < GSL_LOG_DBL_MAX) { + const double ey = exp(y); + gsl_sf_result I1_scaled; + gsl_sf_bessel_I1_scaled_e(x, &I1_scaled); + result->val = ey * I1_scaled.val; + result->err = ey * I1_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_I1_scaled(const double x) +{ + EVAL_RESULT(gsl_sf_bessel_I1_scaled_e(x, &result)); +} + +double gsl_sf_bessel_I1(const double x) +{ + EVAL_RESULT(gsl_sf_bessel_I1_e(x, &result)); +} |