Mercurial > projects > ldc

comparison tango/tango/math/Bessel.d @ 132:1700239cab2e trunk

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[svn r136] MAJOR UNSTABLE UPDATE!!! Initial commit after moving to Tango instead of Phobos. Lots of bugfixes... This build is not suitable for most things.
author lindquist
date Fri, 11 Jan 2008 17:57:40 +0100
parents
children
comparison
equal deleted inserted replaced
131:5825d48b27d1 132:1700239cab2e
1 /**
2 * Cylindrical Bessel functions of integral order.
3 *
4 * Copyright: Based on the CEPHES math library, which is
5 * Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com).
6 * License: BSD style: $(LICENSE)
7 * Authors: Stephen L. Moshier (original C code). Conversion to D by Don Clugston
8 */
9
10 module tango.math.Bessel;
11
12 import tango.math.Math;
13 private import tango.math.IEEE;
14
15
16 private { // Rational polynomial approximations to j0, y0, j1, y1.
17
18 // sqrt(j0^2(1/x^2) + y0^2(1/x^2)) = z P(z**2)/Q(z**2), z(x) = 1/sqrt(x)
19 // Peak error = 1.80e-20
20 const real j0modulusn[] = [ 0x1.154700ea96e79656p-7, 0x1.72244b6e998cd6fp-4,
21 0x1.6ebccf42e9c19fd2p-1, 0x1.6bd844e89cbd639ap+1, 0x1.e812b377c75ebc96p+2,
22 0x1.46d69ca24ce76686p+3, 0x1.b756f7234cc67146p+2, 0x1.943a7471eaa50ab2p-2
23 ];
24
25 const real j0modulusd[] = [ 0x1.5b84007c37011506p-7, 0x1.cfe76758639bdab4p-4,
26 0x1.cbfa09bf71bafc7ep-1, 0x1.c8eafb3836f2eeb4p+1, 0x1.339db78060eb706ep+3,
27 0x1.a06530916be8bc7ap+3, 0x1.23bfe7f67a54893p+3, 1.0
28 ];
29
30
31 // atan(y0(x)/j0(x)) = x - pi/4 + x P(x**2)/Q(x**2)
32 // Peak error = 2.80e-21. Relative error spread = 5.5e-1
33 const real j0phasen[] = [ -0x1.ccbaf3865bb0985ep-22, -0x1.3a6b175e64bdb82ep-14,
34 -0x1.06124b5310cdca28p-8, -0x1.3cebb7ab09cf1b14p-4, -0x1.00156ed209b43c6p-1,
35 -0x1.78aa9ba4254ca20cp-1
36 ];
37
38 const real j0phased[] = [ 0x1.ccbaf3865bb09918p-19, 0x1.3b5b0e12900d58b8p-11,
39 0x1.0897373ff9906f7ep-5, 0x1.450a5b8c552ade4ap-1, 0x1.123e263e7f0e96d2p+2,
40 0x1.d82ecca5654be7d2p+2, 1.0
41 ];
42
43
44 // j1(x) = (x^2-r0^2)(x^2-r1^2)(x^2-r2^2) x P(x**2)/Q(x**2), 0 <= x <= 9
45 // Peak error = 2e-21
46 const real j1n[] = [ -0x1.2f494fa4e623b1bp+58, 0x1.8289f0a5f1e1a784p+52,
47 -0x1.9d773ee29a52c6d8p+45, 0x1.e9394ff57a46071cp+37, -0x1.616c7939904a359p+29,
48 0x1.424414b9ee5671eap+20, -0x1.6db34a9892d653e6p+10, 0x1.dcd7412d90a0db86p-1,
49 -0x1.1444a1643199ee5ep-12
50 ];
51
52 const real j1d[] = [ 0x1.5a1e0a45eb67bacep+75, 0x1.35ee485d62f0ccbap+68,
53 0x1.11ee7aad4e4bcd8p+60, 0x1.3adde5dead800244p+51, 0x1.041c413dfbab693p+42,
54 0x1.4066d12193fcc082p+32, 0x1.24309d0dc2c4d42ep+22, 0x1.7115bea028dd75f2p+11,
55 1.0
56 ];
57
58 // sqrt(j1^2(1/x^2) + y1^2(1/x^2)) = z P(z**2)/Q(z**2), z(x) = 1/sqrt(x)
59 // Peak error = 1.35e=20, Relative error spread = 9.9e0
60 const real [] j1modulusn = [ 0x1.059262020bf7520ap-6, 0x1.012ffc4d1f5cdbc8p-3,
61 0x1.03c17ce18cae596p+0, 0x1.6e0414a7114ae3ccp+1, 0x1.cb047410d229cbc4p+2,
62 0x1.4385d04bb718faaap+1, 0x1.914074c30c746222p-2, -0x1.42abe77f6b307aa6p+2
63 ];
64
65 const real [] j1modulusd = [ 0x1.47d4e6ad98d8246ep-6, 0x1.42562f48058ff904p-3,
66 0x1.44985e2af35c6f9cp+0, 0x1.c6f4a03469c4ef6cp+1, 0x1.1829a060e8d604cp+3,
67 0x1.44111c892f9cc84p+1, -0x1.d7c36d7f1e5aef6ap-1, -0x1.8eeafb1ac81c4c06p+2,
68 1.0
69 ];
70
71 // atan(y1(x)/j1(x)) = x - 3pi/4 + z P(z**2)/Q(z**2), z(x) = 1/x
72 // Peak error = 4.83e-21. Relative error spread = 1.9e0
73 const real [] j1phasen = [ 0x1.ca9f612d83aaa818p-20, 0x1.2e82fcfb7d0fee9ep-12,
74 0x1.e28858c1e947506p-7, 0x1.12b8f96e5173d20ep-2, 0x1.965e6a013154c0ap+0,
75 0x1.0156a25eaa0dd78p+1
76 ];
77
78 const real [] j1phased = [ 0x1.31bf961e57c71ae4p-18, 0x1.9464d8f2abf750a6p-11,
79 0x1.446a786bac2131fp-5, 0x1.76caa8513919873cp-1, 0x1.2130b56bc1a563e4p+2,
80 0x1.b3cc1a865259dfc6p+2, 0x1p+0
81 ];
82
83 }
84
85 /***
86 * Bessel function of order zero
87 *
88 * Returns Bessel function of first kind, order zero of the argument.
89 */
90
91 /* The domain is divided into the intervals [0, 9] and
92 * (9, infinity). In the first interval the rational approximation
93 * is (x^2 - r^2) (x^2 - s^2) (x^2 - t^2) P7(x^2) / Q8(x^2),
94 * where r, s, t are the first three zeros of the function.
95 * In the second interval the expansion is in terms of the
96 * modulus M0(x) = sqrt(J0(x)^2 + Y0(x)^2) and phase P0(x)
97 * = atan(Y0(x)/J0(x)). M0 is approximated by sqrt(1/x)P7(1/x)/Q7(1/x).
98 * The approximation to J0 is M0 * cos(x - pi/4 + 1/x P5(1/x^2)/Q6(1/x^2)).
99 */
100 real cylBessel_j0(real x)
101 {
102
103 // j0(x) = (x^2-JZ1)(x^2-JZ2)(x^2-JZ3)P(x**2)/Q(x**2), 0 <= x <= 9
104 // Peak error = 8.49e-22. Relative error spread = 2.2e-3
105 const real j0n[] = [ -0x1.3e8ff72b890d72d8p+59, 0x1.cc86e3755a4c803p+53,
106 -0x1.0ea6f5bac6623616p+47, 0x1.532c6d94d36f2874p+39, -0x1.ef25a232f6c00118p+30,
107 0x1.aa0690536c11fc2p+21, -0x1.94e67651cc57535p+11, 0x1.4bfe47ac8411eeb2p+0
108 ];
109
110 const real j0d[] = [ 0x1.0096dec5f6560158p+73, 0x1.11705db14995fb9cp+66,
111 0x1.220a41c3daaa7a58p+58, 0x1.93c6b48d196c1082p+49, 0x1.9814684a10dbfda2p+40,
112 0x1.36f20ec527fccda4p+31, 0x1.634596b9247fc34p+21, 0x1.1d3eb73f90657bfcp+11,
113 1.0
114 ];
115 real xx, y, z, modulus, phase;
116
117 xx = x * x;
118 if ( xx < 81.0L ) {
119 const real [] JZ = [5.783185962946784521176L,
120 30.47126234366208639908L, 7.488700679069518344489e1L];
121 y = (xx - JZ[0]) * (xx - JZ[1]) * (xx - JZ[2]);
122 y *= rationalPoly( xx, j0n, j0d);
123 return y;
124 }
125
126 y = fabs(x);
127 xx = 1.0/xx;
128 phase = rationalPoly( xx, j0phasen, j0phased);
129
130 z = 1.0/y;
131 modulus = rationalPoly( z, j0modulusn, j0modulusd);
132
133 y = modulus * cos( y - PI_4 + z*phase) / sqrt(y);
134 return y;
135 }
136
137 /**
138 * Bessel function of the second kind, order zero
139 * Also known as the cylindrical Neumann function, order zero.
140 *
141 * Returns Bessel function of the second kind, of order
142 * zero, of the argument.
143 */
144 real cylBessel_y0(real x)
145 {
146 /* The domain is divided into the intervals [0, 5>, [5,9> and
147 * [9, infinity). In the first interval a rational approximation
148 * R(x) is employed to compute y0(x) = R(x) + 2/pi * log(x) * j0(x).
149 *
150 * In the second interval, the approximation is
151 * (x - p)(x - q)(x - r)(x - s)P7(x)/Q7(x)
152 * where p, q, r, s are zeros of y0(x).
153 *
154 * The third interval uses the same approximations to modulus
155 * and phase as j0(x), whence y0(x) = modulus * sin(phase).
156 */
157
158 // y0(x) = 2/pi * log(x) * j0(x) + P(z**2)/Q(z**2), 0 <= x <= 5
159 // Peak error = 8.55e-22. Relative error spread = 2.7e-1
160 const real y0n[] = [ -0x1.068026b402e2bf7ap+54, 0x1.3a2f7be8c4c8a03ep+55,
161 -0x1.89928488d6524792p+51, 0x1.3e3ea2846f756432p+46, -0x1.c8be8d9366867c78p+39,
162 0x1.43879530964e5fbap+32, -0x1.bee052fef72a5d8p+23, 0x1.e688c8fe417c24d8p+13
163 ];
164
165 const real y0d[] = [ 0x1.bc96c5351e564834p+57, 0x1.6821ac3b4c5209a6p+51,
166 0x1.27098b571836ce64p+44, 0x1.41870d2a9b90aa76p+36, 0x1.00394fd321f52f48p+28,
167 0x1.317ce3b16d65b27p+19, 0x1.0432b36efe4b20aep+10, 1.0
168 ];
169
170 // y0(x) = (x-Y0Z1)(x-Y0Z2)(x-Y0Z3)(x-Y0Z4)P(x)/Q(x), 4.5 <= x <= 9
171 // Peak error = 2.35e-20. Relative error spread = 7.8e-13
172 const real y059n[] = [ -0x1.0fce17d26a21f218p+19, -0x1.c6fc144765fdfaa8p+16,
173 0x1.3e20237c53c7180ep+19, 0x1.7d14055ff6a493c4p+17, 0x1.b8b694729689d1f4p+12,
174 -0x1.1e24596784b6c5cp+12, 0x1.35189cb3ece7ab46p+6, 0x1.9428b3f406b4aa08p+4,
175 -0x1.791187b68dd4240ep+0, 0x1.8417216d568b325ep-6
176 ];
177
178 const real y059d[] = [ 0x1.17af71a3d4167676p+30, 0x1.a36abbb668c79d6cp+31,
179 -0x1.4ff64a14ed73c4d6p+29, 0x1.9d427af195244ffep+26, -0x1.4e85bbbc8d2fd914p+23,
180 0x1.ac59b523ae0bd16cp+19, -0x1.8ebda33eaac74518p+15, 0x1.16194a051cd55a12p+11,
181 -0x1.f2d714ab48d1bd7ep+5, 1.0
182 ];
183
184
185 real xx, y, z, modulus, phase;
186
187 if ( x < 0.0 ) return -real.max;
188 xx = x * x;
189 if ( xx < 81.0L ) {
190 if ( xx < 20.25L ) {
191 y = M_2_PI * log(x) * cylBessel_j0(x);
192 y += rationalPoly( xx, y0n, y0d);
193 } else {
194 const real [] Y0Z = [3.957678419314857868376e0L, 7.086051060301772697624e0L,
195 1.022234504349641701900e1L, 1.336109747387276347827e1L];
196 y = (x - Y0Z[0])*(x - Y0Z[1])*(x - Y0Z[2])*(x - Y0Z[3]);
197 y *= rationalPoly( x, y059n, y059d);
198 }
199 return y;
200 }
201
202 y = fabs(x);
203 xx = 1.0/xx;
204 phase = rationalPoly( xx, j0phasen, j0phased);
205
206 z = 1.0/y;
207 modulus = rationalPoly( z, j0modulusn, j0modulusd);
208
209 y = modulus * sin( y - PI_4 + z*phase) / sqrt(y);
210 return y;
211 }
212
213 /**
214 * Bessel function of order one
215 *
216 * Returns Bessel function of order one of the argument.
217 */
218 real cylBessel_j1(real x)
219 {
220 /* The domain is divided into the intervals [0, 9] and
221 * (9, infinity). In the first interval the rational approximation
222 * is (x^2 - r^2) (x^2 - s^2) (x^2 - t^2) x P8(x^2) / Q8(x^2),
223 * where r, s, t are the first three zeros of the function.
224 * In the second interval the expansion is in terms of the
225 * modulus M1(x) = sqrt(J1(x)^2 + Y1(x)^2) and phase P1(x)
226 * = atan(Y1(x)/J1(x)). M1 is approximated by sqrt(1/x)P7(1/x)/Q8(1/x).
227 * The approximation to j1 is M1 * cos(x - 3 pi/4 + 1/x P5(1/x^2)/Q6(1/x^2)).
228 */
229
230 real xx, y, z, modulus, phase;
231
232 xx = x * x;
233 if ( xx < 81.0L ) {
234 const real [] JZ = [1.46819706421238932572e1L,
235 4.92184563216946036703e1L, 1.03499453895136580332e2L];
236 y = (xx - JZ[0]) * (xx - JZ[1]) * (xx - JZ[2]);
237 y *= x * poly( xx, j1n) / poly( xx, j1d);
238 return y;
239 }
240 y = fabs(x);
241 xx = 1.0/xx;
242 phase = rationalPoly( xx, j1phasen, j1phased);
243 z = 1.0/y;
244 modulus = rationalPoly( z, j1modulusn, j1modulusd);
245
246 const real M_3PI_4 = 3 * PI_4;
247
248 y = modulus * cos( y - M_3PI_4 + z*phase) / sqrt(y);
249 if( x < 0 )
250 y = -y;
251 return y;
252 }
253
254 /**
255 * Bessel function of the second kind, order zero
256 *
257 * Returns Bessel function of the second kind, of order
258 * zero, of the argument.
259 */
260 real cylBessel_y1(real x)
261 in {
262 assert(x>=0.0);
263 // TODO: should it return -infinity for x<0 ?
264 }
265 body {
266 /* The domain is divided into the intervals [0, 4.5>, [4.5,9> and
267 * [9, infinity). In the first interval a rational approximation
268 * R(x) is employed to compute y0(x) = R(x) + 2/pi * log(x) * j0(x).
269 *
270 * In the second interval, the approximation is
271 * (x - p)(x - q)(x - r)(x - s)P9(x)/Q10(x)
272 * where p, q, r, s are zeros of y1(x).
273 *
274 * The third interval uses the same approximations to modulus
275 * and phase as j1(x), whence y1(x) = modulus * sin(phase).
276 *
277 * ACCURACY:
278 *
279 * Absolute error, when y0(x) < 1; else relative error:
280 *
281 * arithmetic domain # trials peak rms
282 * IEEE 0, 30 36000 2.7e-19 5.3e-20
283 *
284 */
285
286 // y1(x) = 2/pi * (log(x) * j1(x) - 1/x) + R(x^2) z P(z**2)/Q(z**2)
287 // 0 <= x <= 4.5, z(x) = x
288 // Peak error = 7.25e-22. Relative error spread = 4.5e-2
289 const real [] y1n = [ -0x1.32cab2601090742p+54, 0x1.432ceb7a8eaeff16p+52,
290 -0x1.bcebec5a2484d3fap+47, 0x1.cc58f3cb54d6ac66p+41, -0x1.b1255e154d0eec0ep+34,
291 0x1.7a337df43298a7c8p+26, -0x1.f77a1afdeff0b62cp+16
292 ];
293
294 const real [] y1d = [ 0x1.8733bcfd7236e604p+56, 0x1.5af412c672fd18d4p+50,
295 0x1.394ba130685755ep+43, 0x1.7b3321523b24afcp+35, 0x1.52946dac22f61d0cp+27,
296 0x1.c9040c6053de5318p+18, 0x1.be5156e6771dba34p+9, 1.0
297 ];
298
299
300 // y1(x) = (x-YZ1)(x-YZ2)(x-YZ3)(x-YZ4)R(x) P(z)/Q(z)
301 // z(x) = x, 4.5 <= x <= 9
302 // Peak error = 3.27e-22. Relative error spread = 4.5e-2
303 const real y159n[] = [ 0x1.2fed87b1e60aa736p+18, -0x1.1a2b18cdb2d1ec5ep+20,
304 -0x1.b848827f47b47022p+20, -0x1.b2e422305ea19a86p+20,
305 -0x1.e3f82ac304534676p+16, 0x1.47a2cb5e852d657ep+14, 0x1.81b2fc6e44d7be8p+12,
306 -0x1.cd861d7b090dd22ep+9, 0x1.588897d683cbfbe2p+5, -0x1.5c7feccf76856bcap-1
307 ];
308
309 const real y159d[] = [ 0x1.9b64f2a4d5614462p+26, -0x1.17501e0e38db675ap+30,
310 0x1.fe88b567c2911c1cp+31, -0x1.86b1781e04e748d4p+29, 0x1.ccd7d4396f2edbcap+26,
311 -0x1.694110c682e5cbcap+23, 0x1.c20f7005b88c789ep+19, -0x1.983a5b4275ab7da8p+15,
312 0x1.17c60380490fa1fcp+11, -0x1.ee84c254392634d8p+5, 1.0
313 ];
314
315 real xx, y, z, modulus, phase;
316
317 z = 1.0/x;
318 xx = x * x;
319 if ( xx < 81.0L ) {
320 if ( xx < 20.25L ) {
321 y = M_2_PI * (log(x) * cylBessel_j1(x) - z);
322 y += x * poly( xx, y1n) / poly( xx, y1d);
323 } else {
324 const real [] Y1Z =
325 [ 2.19714132603101703515e0L, 5.42968104079413513277e0L,
326 8.59600586833116892643e0L, 1.17491548308398812434e1L];
327 y = (x - Y1Z[0])*(x - Y1Z[1])*(x - Y1Z[2])*(x - Y1Z[3]);
328 y *= rationalPoly( x, y159n, y159d);
329 }
330 return y;
331 }
332 xx = 1.0/xx;
333 phase = rationalPoly( xx, j1phasen, j1phased);
334 modulus = rationalPoly( z, j1modulusn, j1modulusd);
335
336 const real M_3PI_4 = 3 * PI_4;
337
338 z = modulus * sin( x - M_3PI_4 + z*phase) / sqrt(x);
339 return z;
340 }
341
342 /**
343 * Bessel function of integer order
344 *
345 * Returns Bessel function of order n, where n is a
346 * (possibly negative) integer.
347 *
348 * The ratio of jn(x) to j0(x) is computed by backward
349 * recurrence. First the ratio jn/jn-1 is found by a
350 * continued fraction expansion. Then the recurrence
351 * relating successive orders is applied until j0 or j1 is
352 * reached.
353 *
354 * If n = 0 or 1 the routine for j0 or j1 is called
355 * directly.
356 *
357 * BUGS: Not suitable for large n or x.
358 *
359 */
360 real cylBessel_jn(int n, real x )
361 {
362 real pkm2, pkm1, pk, xk, r, ans;
363 int k, sign;
364
365 if ( n < 0 ) {
366 n = -n;
367 if ( (n & 1) == 0 ) /* -1**n */
368 sign = 1;
369 else
370 sign = -1;
371 } else
372 sign = 1;
373
374 if ( x < 0.0L ) {
375 if ( n & 1 )
376 sign = -sign;
377 x = -x;
378 }
379
380 if ( n == 0 )
381 return sign * cylBessel_j0(x);
382 if ( n == 1 )
383 return sign * cylBessel_j1(x);
384 // BUG: This code from Cephes is fast, but it makes the Wronksian test fail.
385 // (accuracy is 8 bits lower).
386 // But, the problem might lie in the n = 2 case in cylBessel_yn().
387 // if ( n == 2 )
388 // return sign * (2.0L * cylBessel_j1(x) / x - cylBessel_j0(x));
389
390 if ( x < real.epsilon )
391 return 0;
392
393 /* continued fraction */
394 k = 53;
395 pk = 2 * (n + k);
396 ans = pk;
397 xk = x * x;
398
399 do {
400 pk -= 2.0L;
401 ans = pk - (xk/ans);
402 } while( --k > 0 );
403 ans = x/ans;
404
405 /* backward recurrence */
406
407 pk = 1.0L;
408 pkm1 = 1.0L/ans;
409 k = n-1;
410 r = 2 * k;
411
412 do {
413 pkm2 = (pkm1 * r - pk * x) / x;
414 pk = pkm1;
415 pkm1 = pkm2;
416 r -= 2.0L;
417 } while( --k > 0 );
418
419 if ( fabs(pk) > fabs(pkm1) )
420 ans = cylBessel_j1(x)/pk;
421 else
422 ans = cylBessel_j0(x)/pkm1;
423 return sign * ans;
424 }
425
426 /**
427 * Bessel function of second kind of integer order
428 *
429 * Returns Bessel function of order n, where n is a
430 * (possibly negative) integer.
431 *
432 * The function is evaluated by forward recurrence on
433 * n, starting with values computed by the routines
434 * cylBessel_y0() and cylBessel_y1().
435 *
436 * If n = 0 or 1 the routine for cylBessel_y0 or cylBessel_y1 is called
437 * directly.
438 */
439 real cylBessel_yn(int n, real x)
440 in {
441 assert(x>0); // TODO: should it return -infinity for x<=0 ?
442 }
443 body {
444 real an, r;
445 int k, sign;
446
447 if ( n < 0 ) {
448 n = -n;
449 if ( (n & 1) == 0 ) /* -1**n */
450 sign = 1;
451 else
452 sign = -1;
453 } else
454 sign = 1;
455
456 if ( n == 0 )
457 return sign * cylBessel_y0(x);
458 if ( n == 1 )
459 return sign * cylBessel_y1(x);
460
461 /* forward recurrence on n */
462 real anm2 = cylBessel_y0(x);
463 real anm1 = cylBessel_y1(x);
464 k = 1;
465 r = 2 * k;
466 do {
467 an = r * anm1 / x - anm2;
468 anm2 = anm1;
469 anm1 = an;
470 r += 2.0L;
471 ++k;
472 } while( k < n );
473 return sign * an;
474 }
475
476 private {
477 // Evaluate Chebyshev series
478 double evalCheby(double x, double [] poly)
479 {
480 double b0, b1, b2;
481
482 b0 = poly[$-1];
483 b1 = 0.0;
484 for (int i=poly.length-1; i>=0; --i) {
485 b2 = b1;
486 b1 = b0;
487 b0 = x * b1 - b2 + poly[i];
488 }
489 return 0.5*(b0-b2);
490 }
491 }
492
493 /**
494 * Modified Bessel function of order zero
495 *
496 * Returns modified Bessel function of order zero of the
497 * argument.
498 *
499 * The function is defined as i0(x) = j0( ix ).
500 *
501 * The range is partitioned into the two intervals [0,8] and
502 * (8, infinity). Chebyshev polynomial expansions are employed
503 * in each interval.
504 */
505 double cylBessel_i0(double x)
506 {
507 // Chebyshev coefficients for exp(-x) I0(x) in the interval [0,8].
508 // lim(x->0){ exp(-x) I0(x) } = 1.
509 const double [] A = [ 0x1.5a84e9035a22ap-1, -0x1.37febc057cd8dp-2,
510 0x1.5f7ac77ac88c0p-3, -0x1.84b70342d06eap-4, 0x1.93e8acea8a32dp-5,
511 -0x1.84e9ef121b6f0p-6, 0x1.59961f3dde3ddp-7, -0x1.1b65e201aa849p-8,
512 0x1.adc758a12100ep-10, -0x1.2e2fd1f15eb52p-11, 0x1.8b51b74107cabp-13,
513 -0x1.e2b2659c41d5ap-15, 0x1.13f58be9a2859p-16, -0x1.2866fcba56427p-18,
514 0x1.2bf24978cf4acp-20, -0x1.1ec638f227f8dp-22, 0x1.03b769d4d6435p-24,
515 -0x1.beaf68c0b30abp-27, 0x1.6d903a454cb34p-29, -0x1.1d4fe13ae9556p-31,
516 0x1.a98becc743c10p-34, -0x1.2fc957a946abcp-36, 0x1.9fe2fe19bd324p-39,
517 -0x1.1164c62ee1af0p-41, 0x1.59b464b262627p-44, -0x1.a5022c297fbebp-47,
518 0x1.ee6d893f65ebap-50, -0x1.184eb721ebbb4p-52, 0x1.33362977da589p-55,
519 -0x1.45cb72134d0efp-58 ];
520
521 // Chebyshev coefficients for exp(-x) sqrt(x) I0(x)
522 // in the inverted interval [8,infinity].
523 // lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi).
524 const double [] B = [ 0x1.9be62aca809cbp-1, 0x1.b998ca2e59049p-9,
525 0x1.20fa378999e52p-14, 0x1.8412bc101c586p-19, 0x1.b8007d9cd616ep-23,
526 0x1.8569280d6d56dp-26, 0x1.d2c64a9225b87p-29, 0x1.0f9ccc0f46f75p-31,
527 0x1.a24feabe8004fp-37, -0x1.1511d08397425p-35, -0x1.d0fd7357e7bf2p-37,
528 -0x1.f904303178d66p-40, 0x1.94347fa268cecp-41, 0x1.b1c8c6b83c073p-42,
529 0x1.156ff0d5fc545p-46, -0x1.75d99cf68bb32p-45, -0x1.583fe7e65629ap-47,
530 0x1.12a919094e6d7p-48, 0x1.fee7da3eafb1fp-50, -0x1.8aee7d908de38p-52,
531 -0x1.4600babd21fe4p-52, 0x1.3f3dd076041cdp-55, 0x1.9be1812d98421p-55,
532 -0x1.646da66119130p-58, -0x1.0adb754ca8b19p-57 ];
533
534 double y;
535
536 if (x < 0)
537 x = -x;
538 if (x <= 8.0) {
539 y = (x/2.0) - 2.0;
540 return exp(x) * evalCheby( y, A);
541 }
542 return exp(x) * evalCheby( 32.0/x - 2.0, B) / sqrt(x);
543 }
544
545 /**
546 * Modified Bessel function of order one
547 *
548 * Returns modified Bessel function of order one of the
549 * argument.
550 *
551 * The function is defined as i1(x) = -i j1( ix ).
552 *
553 * The range is partitioned into the two intervals [0,8] and
554 * (8, infinity). Chebyshev polynomial expansions are employed
555 * in each interval.
556 */
557 double cylBessel_i1(double x)
558 {
559 const double [] A = [ 0x1.02a63724a7ffap-2, -0x1.694d10469192ep-3,
560 0x1.a46dad536f53cp-4, -0x1.b1bbc537c9ebcp-5, 0x1.951e3e7bb2349p-6,
561 -0x1.5a29f7913a26ap-7, 0x1.1065349d3a1b4p-8, -0x1.8cc620b3cd4a4p-10,
562 0x1.0c95db6c6df7dp-11, -0x1.533cad3d694fep-13, 0x1.911b542c70d0bp-15,
563 -0x1.bd5f9b8debbcfp-17, 0x1.d1c4ed511afc5p-19, -0x1.cc0798363992ap-21,
564 0x1.ae344b347d108p-23, -0x1.7dd3e24b8c3e8p-25, 0x1.4258e02395010p-27,
565 -0x1.0361b28ea67e6p-29, 0x1.8ea34b43fdf6cp-32, -0x1.2510397eb07dep-34,
566 0x1.9cee2b21d3154p-37, -0x1.173835fb70366p-39, 0x1.6af784779d955p-42,
567 -0x1.c628e1c8f0b3bp-45, 0x1.11d7f0615290cp-47, -0x1.3eaaa7e0d1573p-50,
568 0x1.663e3e593bfacp-53, -0x1.857d0c38a0576p-56, 0x1.99f2a0c3c4014p-59
569 ];
570
571 // Chebyshev coefficients for exp(-x) sqrt(x) I1(x)
572 // in the inverted interval [8,infinity].
573 // lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi).
574 const double [] B = [ 0x1.8ea18b55b1514p-1, -0x1.3fda053fcdb4cp-7,
575 -0x1.cfd7f804aa9a6p-14, -0x1.048df49ca0373p-18, -0x1.0dbfd2e9e5443p-22,
576 -0x1.c415394bb46c1p-26, -0x1.0790b9ad53528p-28, -0x1.334ca5423dd80p-31,
577 -0x1.4dcf9d4504c0cp-36, 0x1.1e1a1f1587865p-35, 0x1.f101f653c457bp-37,
578 0x1.1e7d3f6439fa3p-39, -0x1.953e1076ab493p-41, -0x1.cbc458e73e255p-42,
579 -0x1.7a9482e6d22a0p-46, 0x1.80d3c26b3281ep-45, 0x1.776e1762d31e8p-47,
580 -0x1.12db5138afbc7p-48, -0x1.0efcd8bc4d22ap-49, 0x1.7d68e5f04a2d1p-52,
581 0x1.55915fceb588ap-52, -0x1.2806c9c773320p-55, -0x1.acea3b2532277p-55,
582 0x1.45b8aea87b950p-58, 0x1.1556db352e8e6p-57 ];
583
584 double y, z;
585
586 z = fabs(x);
587 if( z <= 8.0 ) {
588 y = (z/2.0) - 2.0;
589 z = evalCheby( y, A ) * z * exp(z);
590 } else {
591 z = exp(z) * evalCheby( 32.0/z - 2.0, B ) / sqrt(z);
592 }
593 if (x < 0.0 )
594 z = -z;
595 return z;
596 }
597
598 debug(UnitTest) {
599
600 unittest {
601 // argument, result1, result2, derivative. Correct result is result1+result2.
602 const real [4][] j0_test_points = [
603 [8.0L, 1.71646118164062500000E-1L, 4.68897349140609086941E-6L, -2.34636346853914624381E-1L],
604 [4.54541015625L, -3.09783935546875000000E-1L, 7.07472668157686463367E-6L, 2.42993657373627558460E-1L],
605 [2.85711669921875L, -2.07901000976562500000E-1L, 1.15237285263902751582E-5L, -3.90402225324501311651E-1L],
606 [2.0L, 2.23876953125000000000E-1L, 1.38260162356680518275E-5L, -5.76724807756873387202E-1L],
607 [1.16415321826934814453125e-10L, 9.99984741210937500000E-1L, 1.52587890624999966119E-5L,
608 9.99999999999999999997E-1L],
609 [-2.0L, 2.23876953125000000000E-1L,
610 1.38260162356680518275E-5L, 5.76724807756873387202E-1L]
611 ];
612
613 const real [4][] y0_test_points = [
614 [ 8.0L, 2.23510742187500000000E-1L, 1.07472000662205273234E-5L, 1.58060461731247494256E-1L],
615 [4.54541015625L, -2.08114624023437500000E-1L, 1.45018823856668874574E-5L, -2.88887645307401250876E-1L],
616 [2.85711669921875L, 4.20303344726562500000E-1L, 1.32781607563122276008E-5L, -2.82488638474982469213E-1],
617 [2.0L, 5.10360717773437500000E-1L, 1.49548763076195966066E-5L, 1.07032431540937546888E-1L],
618 [1.16415321826934814453125e-10L, -1.46357574462890625000E1L, 3.54110537011061127637E-6L,
619 5.46852220461145271913E9L]
620 ];
621
622 const real [4][] j1_test_points = [
623 [ 8.0L, 2.34634399414062500000E-1L, 1.94743985212438127665E-6L,1.42321263780814578043E-1],
624 [4.54541015625L, -2.42996215820312500000E-1L, 2.55844668494153980076E-6L, -2.56317734136211337012E-1],
625 [2.85711669921875L, 3.90396118164062500000E-1L, 6.10716043881165077013E-6L, -3.44531507106757980441E-1L],
626 [2.0L, 5.76721191406250000000E-1L, 3.61635062338720244824E-6L, -6.44716247372010255494E-2L],
627 [1.16415321826934814453125e-10L, 5.820677273504770710133016109466552734375e-11L,
628 8.881784197001251337312921818461805735896e-16L, 4.99999999999999999997E-1L],
629 [-2.0L, -5.76721191406250000000E-1L, -3.61635062338720244824E-6L, -6.44716247372010255494E-2L]
630 ];
631
632 const real [4][] y1_test_points = [
633 [8.0L, -1.58065795898437500000E-1L,
634 5.33416719000574444473E-6L, 2.43279047103972157309E-1L],
635 [4.54541015625L, 2.88879394531250000000E-1L,
636 8.25077615125087585195E-6L, -2.71656024771791736625E-1L],
637 [2.85711669921875L, 2.82485961914062500000E-1,
638 2.67656091996921314433E-6L, 3.21444694221532719737E-1],
639 [2.0L, -1.07040405273437500000E-1L,
640 7.97373249995311162923E-6L, 5.63891888420213893041E-1],
641 [1.16415321826934814453125e-10L, -5.46852220500000000000E9L,
642 3.88547280871200700671E-1L, 4.69742480525120196168E19L]
643 ];
644
645 foreach(real [4] t; j0_test_points) {
646 assert(feqrel(cylBessel_j0(t[0]), t[1]+t[2]) >=real.mant_dig-3);
647 }
648
649 foreach(real [4] t; y0_test_points) {
650 assert(feqrel(cylBessel_y0(t[0]), t[1]+t[2]) >=real.mant_dig-4);
651 }
652 foreach(real [4] t; j1_test_points) {
653 assert(feqrel(cylBessel_j1(t[0]), t[1]+t[2]) >=real.mant_dig-3);
654 }
655
656 foreach(real [4] t; y1_test_points) {
657 assert(feqrel(cylBessel_y1(t[0]), t[1]+t[2]) >=real.mant_dig-4);
658 }
659
660 // Values from MS Excel, of doubtful accuracy.
661 assert(fabs(-0.060_409_940_421_649 - cylBessel_j0(173.5)) < 0.000_000_000_1);
662 assert(fabs(-0.044_733_447_576_5866 - cylBessel_y0(313.25)) < 0.000_000_000_1);
663 assert(fabs(0.00391280088318945 - cylBessel_j1(123.25)) < 0.000_000_000_1);
664 assert(fabs(-0.0648628570878951 - cylBessel_j1(-91)) < 0.000_000_000_1);
665 assert(fabs(-0.0759578537652805 - cylBessel_y1(107.75)) < 0.000_000_000_1);
666
667 assert(fabs(13.442_456_516_6771-cylBessel_i0(4.2)) < 0.000_001);
668 assert(fabs(1.6500020842093e+28-cylBessel_i0(-68)) < 0.000_001e+28);
669 assert(fabs(4.02746515903173e+10-cylBessel_i1(27)) < 0.000_001e+10);
670 assert(fabs(-2.83613942886386e-02-cylBessel_i1(-0.0567)) < 0.000_000_001e-2);
671 }
672 }
673
674 debug(UnitTest) {
675
676 unittest {
677
678 // Wronksian test for Bessel functions
679 void testWronksian(int n, real x)
680 {
681 real Jnp1 = cylBessel_jn(n + 1, x);
682 real Jmn = cylBessel_jn(-n, x);
683 real Jn = cylBessel_jn(n, x);
684 real Jmnp1 = cylBessel_jn(-(n + 1), x);
685 /* This should be trivially zero. */
686 assert( fabs(Jnp1 * Jmn + Jn * Jmnp1) == 0);
687 if (x < 0.0) {
688 x = -x;
689 Jn = cylBessel_jn(n, x);
690 Jnp1 = cylBessel_jn(n + 1, x);
691 }
692 real Yn = cylBessel_yn(n, x);
693 real Ynp1 = cylBessel_yn(n + 1, x);
694 /* The Wronksian. */
695 real w1 = Jnp1 * Yn - Jn * Ynp1;
696 /* What the Wronksian should be. */
697 real w2 = 2.0 / (PI * x);
698
699 real reldif = feqrel(w1, w2);
700 assert(reldif >= real.mant_dig-6);
701 }
702
703 real delta;
704 int n, i, j;
705
706 delta = 0.6 / PI;
707 for (n = -30; n <= 30; n++) {
708 real x = -30.0;
709 while (x < 30.0) {
710 testWronksian (n, x);
711 x += delta;
712 }
713 delta += .00123456;
714 }
715 assert(cylBessel_jn(20, 1e-80)==0);
716
717 // NaN propagation
718 assert(isIdentical(cylBessel_i1(NaN(0xDEF)), NaN(0xDEF)));
719 assert(isIdentical(cylBessel_i0(NaN(0x846)), NaN(0x846)));
720
721 }
722
723 }