1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.numerics.special;
39
40 import static org.apache.commons.math3.util.FastMath.abs;
41 import static org.apache.commons.math3.util.FastMath.exp;
42 import static org.apache.commons.math3.util.FastMath.log;
43 import static org.apache.commons.math3.util.FastMath.sqrt;
44
45
46
47
48
49
50
51 public class ModifiedBessel {
52
53
54
55
56
57
58 private static final double[] A_i0 = {
59 -4.41534164647933937950E-18,
60 3.33079451882223809783E-17,
61 -2.43127984654795469359E-16,
62 1.71539128555513303061E-15,
63 -1.16853328779934516808E-14,
64 7.67618549860493561688E-14,
65 -4.85644678311192946090E-13,
66 2.95505266312963983461E-12,
67 -1.72682629144155570723E-11,
68 9.67580903537323691224E-11,
69 -5.18979560163526290666E-10,
70 2.65982372468238665035E-9,
71 -1.30002500998624804212E-8,
72 6.04699502254191894932E-8,
73 -2.67079385394061173391E-7,
74 1.11738753912010371815E-6,
75 -4.41673835845875056359E-6,
76 1.64484480707288970893E-5,
77 -5.75419501008210370398E-5,
78 1.88502885095841655729E-4,
79 -5.76375574538582365885E-4,
80 1.63947561694133579842E-3,
81 -4.32430999505057594430E-3,
82 1.05464603945949983183E-2,
83 -2.37374148058994688156E-2,
84 4.93052842396707084878E-2,
85 -9.49010970480476444210E-2,
86 1.71620901522208775349E-1,
87 -3.04682672343198398683E-1,
88 6.76795274409476084995E-1
89 };
90
91
92
93
94
95 private static final double[] B_i0 = {
96 -7.23318048787475395456E-18,
97 -4.83050448594418207126E-18,
98 4.46562142029675999901E-17,
99 3.46122286769746109310E-17,
100 -2.82762398051658348494E-16,
101 -3.42548561967721913462E-16,
102 1.77256013305652638360E-15,
103 3.81168066935262242075E-15,
104 -9.55484669882830764870E-15,
105 -4.15056934728722208663E-14,
106 1.54008621752140982691E-14,
107 3.85277838274214270114E-13,
108 7.18012445138366623367E-13,
109 -1.79417853150680611778E-12,
110 -1.32158118404477131188E-11,
111 -3.14991652796324136454E-11,
112 1.18891471078464383424E-11,
113 4.94060238822496958910E-10,
114 3.39623202570838634515E-9,
115 2.26666899049817806459E-8,
116 2.04891858946906374183E-7,
117 2.89137052083475648297E-6,
118 6.88975834691682398426E-5,
119 3.36911647825569408990E-3,
120 8.04490411014108831608E-1
121 };
122
123
124
125
126
127 private static final double[] A_i1 = {
128 2.77791411276104639959E-18,
129 -2.11142121435816608115E-17,
130 1.55363195773620046921E-16,
131 -1.10559694773538630805E-15,
132 7.60068429473540693410E-15,
133 -5.04218550472791168711E-14,
134 3.22379336594557470981E-13,
135 -1.98397439776494371520E-12,
136 1.17361862988909016308E-11,
137 -6.66348972350202774223E-11,
138 3.62559028155211703701E-10,
139 -1.88724975172282928790E-9,
140 9.38153738649577178388E-9,
141 -4.44505912879632808065E-8,
142 2.00329475355213526229E-7,
143 -8.56872026469545474066E-7,
144 3.47025130813767847674E-6,
145 -1.32731636560394358279E-5,
146 4.78156510755005422638E-5,
147 -1.61760815825896745588E-4,
148 5.12285956168575772895E-4,
149 -1.51357245063125314899E-3,
150 4.15642294431288815669E-3,
151 -1.05640848946261981558E-2,
152 2.47264490306265168283E-2,
153 -5.29459812080949914269E-2,
154 1.02643658689847095384E-1,
155 -1.76416518357834055153E-1,
156 2.52587186443633654823E-1
157 };
158
159
160
161
162
163 private static final double[] B_i1 = {
164 7.51729631084210481353E-18,
165 4.41434832307170791151E-18,
166 -4.65030536848935832153E-17,
167 -3.20952592199342395980E-17,
168 2.96262899764595013876E-16,
169 3.30820231092092828324E-16,
170 -1.88035477551078244854E-15,
171 -3.81440307243700780478E-15,
172 1.04202769841288027642E-14,
173 4.27244001671195135429E-14,
174 -2.10154184277266431302E-14,
175 -4.08355111109219731823E-13,
176 -7.19855177624590851209E-13,
177 2.03562854414708950722E-12,
178 1.41258074366137813316E-11,
179 3.25260358301548823856E-11,
180 -1.89749581235054123450E-11,
181 -5.58974346219658380687E-10,
182 -3.83538038596423702205E-9,
183 -2.63146884688951950684E-8,
184 -2.51223623787020892529E-7,
185 -3.88256480887769039346E-6,
186 -1.10588938762623716291E-4,
187 -9.76109749136146840777E-3,
188 7.78576235018280120474E-1
189 };
190
191
192
193
194
195
196
197
198
199
200 public static double i0(double x) {
201 double y;
202 if (x < 0) {
203 x = -x;
204 }
205 if (x <= 8.0) {
206 y = (x * 0.5) - 2.0;
207 return (eToThe(x) * evaluateChebyshev(y, A_i0, 30));
208 }
209 double ix = 1.0 / x;
210 return (eToThe(x) * evaluateChebyshev(32.0 * ix - 2.0, B_i0, 25) * sqrt(ix));
211 }
212
213
214
215
216
217
218
219
220
221
222 public static double i1(double x) {
223 double y, z;
224
225 z = abs(x);
226 if (z <= 8.0) {
227 y = (z * 0.5) - 2.0;
228 z = evaluateChebyshev(y, A_i1, 29) * z * eToThe(z);
229 } else {
230 double iz = 1.0 / z;
231 z = eToThe(z) * evaluateChebyshev(32.0 * iz - 2.0, B_i1, 25) * sqrt(iz);
232 }
233 if (x < 0.0) {
234 z = -z;
235 }
236 return (z);
237 }
238
239
240
241
242
243
244
245 public static double i1OverI0(double x) {
246 return (i1(x) / i0(x));
247 }
248
249
250
251
252
253
254
255 public static double lnI0(double x) {
256 return log(i0(x));
257 }
258
259
260
261
262
263
264
265
266
267
268
269
270
271 private static double evaluateChebyshev(double x, double[] coefficients, int N) {
272 double b0, b1, b2;
273
274 int p = 0;
275 int i;
276
277 b0 = coefficients[p++];
278 b1 = 0.0;
279 i = N - 1;
280
281 do {
282 b2 = b1;
283 b1 = b0;
284 b0 = x * b1 - b2 + coefficients[p++];
285 } while (--i > 0);
286
287 return (0.5 * (b0 - b2));
288 }
289
290
291
292
293
294
295
296
297 private static double eToThe(double x) {
298 double res = exp(x);
299 if (res == Double.POSITIVE_INFINITY) {
300 return Double.MAX_VALUE;
301 } else if (res == Double.NEGATIVE_INFINITY) {
302 return Double.MIN_VALUE;
303 }
304 return res;
305 }
306 }