1 /* specfunc/bessel_I1.c
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
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.
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.
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.
20 /* Author: G. Jungman */
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_sf_bessel.h>
29 #include "chebyshev.h"
30 #include "cheb_eval.c"
32 #define ROOT_EIGHT (2.0*M_SQRT2)
35 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
37 /* based on SLATEC besi1(), besi1e() */
39 /* chebyshev expansions
41 series for bi1 on the interval 0. to 9.00000d+00
42 with weighted error 2.40e-17
43 log weighted error 16.62
44 significant figures required 16.23
45 decimal places required 17.14
47 series for ai1 on the interval 1.25000d-01 to 3.33333d-01
48 with weighted error 6.98e-17
49 log weighted error 16.16
50 significant figures required 14.53
51 decimal places required 16.82
53 series for ai12 on the interval 0. to 1.25000d-01
54 with weighted error 3.55e-17
55 log weighted error 16.45
56 significant figures required 14.69
57 decimal places required 17.12
60 static double bi1_data[11] = {
61 -0.001971713261099859,
73 static cheb_series bi1_cs = {
80 static double ai1_data[21] = {
103 static cheb_series ai1_cs = {
110 static double ai12_data[22] = {
112 -0.00976109749136147,
113 -0.00011058893876263,
114 -0.00000388256480887,
115 -0.00000025122362377,
116 -0.00000002631468847,
117 -0.00000000383538039,
118 -0.00000000055897433,
119 -0.00000000001897495,
123 -0.00000000000071985,
124 -0.00000000000040836,
125 -0.00000000000002101,
128 -0.00000000000000382,
129 -0.00000000000000186,
134 static cheb_series ai12_cs = {
142 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
144 int gsl_sf_bessel_I1_scaled_e(const double x, gsl_sf_result * result)
146 const double xmin = 2.0 * GSL_DBL_MIN;
147 const double x_small = ROOT_EIGHT * GSL_SQRT_DBL_EPSILON;
148 const double y = fabs(x);
150 /* CHECK_POINTER(result) */
158 UNDERFLOW_ERROR(result);
160 else if(y < x_small) {
166 const double ey = exp(-y);
168 cheb_eval_e(&bi1_cs, y*y/4.5-1.0, &c);
169 result->val = x * ey * (0.875 + c.val);
170 result->err = ey * c.err + y * GSL_DBL_EPSILON * fabs(result->val);
171 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
175 const double sy = sqrt(y);
179 cheb_eval_e(&ai1_cs, (48.0/y-11.0)/5.0, &c);
180 b = (0.375 + c.val) / sy;
181 s = (x > 0.0 ? 1.0 : -1.0);
183 result->err = c.err / sy;
184 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
188 const double sy = sqrt(y);
192 cheb_eval_e(&ai12_cs, 16.0/y-1.0, &c);
193 b = (0.375 + c.val) / sy;
194 s = (x > 0.0 ? 1.0 : -1.0);
196 result->err = c.err / sy;
197 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
203 int gsl_sf_bessel_I1_e(const double x, gsl_sf_result * result)
205 const double xmin = 2.0 * GSL_DBL_MIN;
206 const double x_small = ROOT_EIGHT * GSL_SQRT_DBL_EPSILON;
207 const double y = fabs(x);
209 /* CHECK_POINTER(result) */
217 UNDERFLOW_ERROR(result);
219 else if(y < x_small) {
226 cheb_eval_e(&bi1_cs, y*y/4.5-1.0, &c);
227 result->val = x * (0.875 + c.val);
228 result->err = y * c.err;
229 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
232 else if(y < GSL_LOG_DBL_MAX) {
233 const double ey = exp(y);
234 gsl_sf_result I1_scaled;
235 gsl_sf_bessel_I1_scaled_e(x, &I1_scaled);
236 result->val = ey * I1_scaled.val;
237 result->err = ey * I1_scaled.err + y * GSL_DBL_EPSILON * fabs(result->val);
238 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
242 OVERFLOW_ERROR(result);
246 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
250 double gsl_sf_bessel_I1_scaled(const double x)
252 EVAL_RESULT(gsl_sf_bessel_I1_scaled_e(x, &result));
255 double gsl_sf_bessel_I1(const double x)
257 EVAL_RESULT(gsl_sf_bessel_I1_e(x, &result));