1 /* specfunc/bessel_I0.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
4 *
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.
9 *
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.
14 *
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.
18 */
19
20 /* Author: G. Jungman */
21
22 #include <config.h>
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_sf_bessel.h>
26
27 #include "error.h"
28
29 #include "chebyshev.h"
30 #include "cheb_eval.c"
31
32 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
33
34
35 /* based on SLATEC besi0 */
36
37 /* chebyshev expansions
38
39 series for bi0 on the interval 0. to 9.00000d+00
40 with weighted error 2.46e-18
41 log weighted error 17.61
42 significant figures required 17.90
43 decimal places required 18.15
44
45 series for ai0 on the interval 1.25000d-01 to 3.33333d-01
46 with weighted error 7.87e-17
47 log weighted error 16.10
48 significant figures required 14.69
49 decimal places required 16.76
50
51
52 series for ai02 on the interval 0. to 1.25000d-01
53 with weighted error 3.79e-17
54 log weighted error 16.42
55 significant figures required 14.86
56 decimal places required 17.09
57 */
58
59 static double bi0_data[12] = {
60 -.07660547252839144951,
61 1.92733795399380827000,
62 .22826445869203013390,
63 .01304891466707290428,
64 .00043442709008164874,
65 .00000942265768600193,
66 .00000014340062895106,
67 .00000000161384906966,
68 .00000000001396650044,
69 .00000000000009579451,
70 .00000000000000053339,
71 .00000000000000000245
72 };
73 static cheb_series bi0_cs = {
74 bi0_data,
75 11,
76 -1, 1,
77 11
78 };
79
80 static double ai0_data[21] = {
81 .07575994494023796,
82 .00759138081082334,
83 .00041531313389237,
84 .00001070076463439,
85 -.00000790117997921,
86 -.00000078261435014,
87 .00000027838499429,
88 .00000000825247260,
89 -.00000001204463945,
90 .00000000155964859,
91 .00000000022925563,
92 -.00000000011916228,
93 .00000000001757854,
94 .00000000000112822,
95 -.00000000000114684,
96 .00000000000027155,
97 -.00000000000002415,
98 -.00000000000000608,
99 .00000000000000314,
100 -.00000000000000071,
101 .00000000000000007
102 };
103 static cheb_series ai0_cs = {
104 ai0_data,
105 20,
106 -1, 1,
107 13
108 };
109
110 static double ai02_data[22] = {
111 .05449041101410882,
112 .00336911647825569,
113 .00006889758346918,
114 .00000289137052082,
115 .00000020489185893,
116 .00000002266668991,
117 .00000000339623203,
118 .00000000049406022,
119 .00000000001188914,
120 -.00000000003149915,
121 -.00000000001321580,
122 -.00000000000179419,
123 .00000000000071801,
124 .00000000000038529,
125 .00000000000001539,
126 -.00000000000004151,
127 -.00000000000000954,
128 .00000000000000382,
129 .00000000000000176,
130 -.00000000000000034,
131 -.00000000000000027,
132 .00000000000000003
133 };
134 static cheb_series ai02_cs = {
135 ai02_data,
136 21,
137 -1, 1,
138 11
139 };
140
141
142 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
143
gsl_sf_bessel_I0_scaled_e(const double x,gsl_sf_result * result)144 int gsl_sf_bessel_I0_scaled_e(const double x, gsl_sf_result * result)
145 {
146 double y = fabs(x);
147
148 /* CHECK_POINTER(result) */
149
150 if(y < 2.0 * GSL_SQRT_DBL_EPSILON) {
151 result->val = 1.0 - y;
152 result->err = 0.5*y*y;
153 return GSL_SUCCESS;
154 }
155 else if(y <= 3.0) {
156 const double ey = exp(-y);
157 gsl_sf_result c;
158 cheb_eval_e(&bi0_cs, y*y/4.5-1.0, &c);
159 result->val = ey * (2.75 + c.val);
160 result->err = GSL_DBL_EPSILON * fabs(result->val) + ey * c.err;
161 return GSL_SUCCESS;
162 }
163 else if(y <= 8.0) {
164 const double sy = sqrt(y);
165 gsl_sf_result c;
166 cheb_eval_e(&ai0_cs, (48.0/y-11.0)/5.0, &c);
167 result->val = (0.375 + c.val) / sy;
168 result->err = 2.0 * GSL_DBL_EPSILON * (0.375 + fabs(c.val)) / sy;
169 result->err += c.err / sy;
170 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
171 return GSL_SUCCESS;
172 }
173 else {
174 const double sy = sqrt(y);
175 gsl_sf_result c;
176 cheb_eval_e(&ai02_cs, 16.0/y-1.0, &c);
177 result->val = (0.375 + c.val) / sy;
178 result->err = 2.0 * GSL_DBL_EPSILON * (0.375 + fabs(c.val)) / sy;
179 result->err += c.err / sy;
180 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
181 return GSL_SUCCESS;
182 }
183 }
184
185
gsl_sf_bessel_I0_e(const double x,gsl_sf_result * result)186 int gsl_sf_bessel_I0_e(const double x, gsl_sf_result * result)
187 {
188 double y = fabs(x);
189
190 /* CHECK_POINTER(result) */
191
192 if(y < 2.0 * GSL_SQRT_DBL_EPSILON) {
193 result->val = 1.0;
194 result->err = 0.5*y*y;
195 return GSL_SUCCESS;
196 }
197 else if(y <= 3.0) {
198 gsl_sf_result c;
199 cheb_eval_e(&bi0_cs, y*y/4.5-1.0, &c);
200 result->val = 2.75 + c.val;
201 result->err = GSL_DBL_EPSILON * (2.75 + fabs(c.val));
202 result->err += c.err;
203 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
204 return GSL_SUCCESS;
205 }
206 else if(y < GSL_LOG_DBL_MAX - 1.0) {
207 const double ey = exp(y);
208 gsl_sf_result b_scaled;
209 gsl_sf_bessel_I0_scaled_e(x, &b_scaled);
210 result->val = ey * b_scaled.val;
211 result->err = ey * b_scaled.err + y*GSL_DBL_EPSILON*fabs(result->val);
212 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
213 return GSL_SUCCESS;
214 }
215 else {
216 OVERFLOW_ERROR(result);
217 }
218 }
219
220 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
221
222 #include "eval.h"
223
gsl_sf_bessel_I0_scaled(const double x)224 double gsl_sf_bessel_I0_scaled(const double x)
225 {
226 EVAL_RESULT(gsl_sf_bessel_I0_scaled_e(x, &result); )
227 }
228
gsl_sf_bessel_I0(const double x)229 double gsl_sf_bessel_I0(const double x)
230 {
231 EVAL_RESULT(gsl_sf_bessel_I0_e(x, &result); )
232 }
233