|
132
|
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 }
|
