diff options
Diffstat (limited to 'gsl-1.9/specfunc/bessel_K0.c')
| -rw-r--r-- | gsl-1.9/specfunc/bessel_K0.c | 214 |
1 files changed, 214 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/bessel_K0.c b/gsl-1.9/specfunc/bessel_K0.c new file mode 100644 index 0000000..ebbb0a1 --- /dev/null +++ b/gsl-1.9/specfunc/bessel_K0.c @@ -0,0 +1,214 @@ +/* specfunc/bessel_K0.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_exp.h> +#include <gsl/gsl_sf_bessel.h> + +#include "error.h" + +#include "chebyshev.h" +#include "cheb_eval.c" + +/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/ + +/* based on SLATEC bk0(), bk0e() */ + +/* chebyshev expansions + + series for bk0 on the interval 0. to 4.00000d+00 + with weighted error 3.57e-19 + log weighted error 18.45 + significant figures required 17.99 + decimal places required 18.97 + + series for ak0 on the interval 1.25000d-01 to 5.00000d-01 + with weighted error 5.34e-17 + log weighted error 16.27 + significant figures required 14.92 + decimal places required 16.89 + + series for ak02 on the interval 0. to 1.25000d-01 + with weighted error 2.34e-17 + log weighted error 16.63 + significant figures required 14.67 + decimal places required 17.20 +*/ + +static double bk0_data[11] = { + -0.03532739323390276872, + 0.3442898999246284869, + 0.03597993651536150163, + 0.00126461541144692592, + 0.00002286212103119451, + 0.00000025347910790261, + 0.00000000190451637722, + 0.00000000001034969525, + 0.00000000000004259816, + 0.00000000000000013744, + 0.00000000000000000035 +}; +static cheb_series bk0_cs = { + bk0_data, + 10, + -1, 1, + 10 +}; + +static double ak0_data[17] = { + -0.07643947903327941, + -0.02235652605699819, + 0.00077341811546938, + -0.00004281006688886, + 0.00000308170017386, + -0.00000026393672220, + 0.00000002563713036, + -0.00000000274270554, + 0.00000000031694296, + -0.00000000003902353, + 0.00000000000506804, + -0.00000000000068895, + 0.00000000000009744, + -0.00000000000001427, + 0.00000000000000215, + -0.00000000000000033, + 0.00000000000000005 +}; +static cheb_series ak0_cs = { + ak0_data, + 16, + -1, 1, + 10 +}; + +static double ak02_data[14] = { + -0.01201869826307592, + -0.00917485269102569, + 0.00014445509317750, + -0.00000401361417543, + 0.00000015678318108, + -0.00000000777011043, + 0.00000000046111825, + -0.00000000003158592, + 0.00000000000243501, + -0.00000000000020743, + 0.00000000000001925, + -0.00000000000000192, + 0.00000000000000020, + -0.00000000000000002 +}; +static cheb_series ak02_cs = { + ak02_data, + 13, + -1, 1, + 8 +}; + + +/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ + +int gsl_sf_bessel_K0_scaled_e(const double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(x <= 0.0) { + DOMAIN_ERROR(result); + } + else if(x <= 2.0) { + const double lx = log(x); + const double ex = exp(x); + int stat_I0; + gsl_sf_result I0; + gsl_sf_result c; + cheb_eval_e(&bk0_cs, 0.5*x*x-1.0, &c); + stat_I0 = gsl_sf_bessel_I0_e(x, &I0); + result->val = ex * ((-lx+M_LN2)*I0.val - 0.25 + c.val); + result->err = ex * ((M_LN2+fabs(lx))*I0.err + c.err); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return stat_I0; + } + else if(x <= 8.0) { + const double sx = sqrt(x); + gsl_sf_result c; + cheb_eval_e(&ak0_cs, (16.0/x-5.0)/3.0, &c); + result->val = (1.25 + c.val) / sx; + result->err = c.err / sx; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + const double sx = sqrt(x); + gsl_sf_result c; + cheb_eval_e(&ak02_cs, 16.0/x-1.0, &c); + result->val = (1.25 + c.val) / sx; + result->err = (c.err + GSL_DBL_EPSILON) / sx; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } +} + + +int gsl_sf_bessel_K0_e(const double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(x <= 0.0) { + DOMAIN_ERROR(result); + } + else if(x <= 2.0) { + const double lx = log(x); + int stat_I0; + gsl_sf_result I0; + gsl_sf_result c; + cheb_eval_e(&bk0_cs, 0.5*x*x-1.0, &c); + stat_I0 = gsl_sf_bessel_I0_e(x, &I0); + result->val = (-lx+M_LN2)*I0.val - 0.25 + c.val; + result->err = (fabs(lx) + M_LN2) * I0.err + c.err; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return stat_I0; + } + else { + gsl_sf_result K0_scaled; + int stat_K0 = gsl_sf_bessel_K0_scaled_e(x, &K0_scaled); + int stat_e = gsl_sf_exp_mult_err_e(-x, GSL_DBL_EPSILON*fabs(x), + K0_scaled.val, K0_scaled.err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_K0); + } +} + + +/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ + +#include "eval.h" + +double gsl_sf_bessel_K0_scaled(const double x) +{ + EVAL_RESULT(gsl_sf_bessel_K0_scaled_e(x, &result)); +} + +double gsl_sf_bessel_K0(const double x) +{ + EVAL_RESULT(gsl_sf_bessel_K0_e(x, &result)); +} + |