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