1/* origin: FreeBSD /usr/src/lib/msun/src/e_jn.c */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunSoft, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12/*
13 * jn(n, x), yn(n, x)
14 * floating point Bessel's function of the 1st and 2nd kind
15 * of order n
16 *
17 * Special cases:
18 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
19 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
20 * Note 2. About jn(n,x), yn(n,x)
21 * For n=0, j0(x) is called,
22 * for n=1, j1(x) is called,
23 * for n<=x, forward recursion is used starting
24 * from values of j0(x) and j1(x).
25 * for n>x, a continued fraction approximation to
26 * j(n,x)/j(n-1,x) is evaluated and then backward
27 * recursion is used starting from a supposed value
28 * for j(n,x). The resulting value of j(0,x) is
29 * compared with the actual value to correct the
30 * supposed value of j(n,x).
31 *
32 * yn(n,x) is similar in all respects, except
33 * that forward recursion is used for all
34 * values of n>1.
35 */
36
37use super::{cos, fabs, get_high_word, get_low_word, j0, j1, log, sin, sqrt, y0, y1};
38
39const INVSQRTPI: f64 = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
40
41/// Integer order of the [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the first kind (f64).
42pub fn jn(n: i32, mut x: f64) -> f64 {
43 let mut ix: u32;
44 let lx: u32;
45 let nm1: i32;
46 let mut i: i32;
47 let mut sign: bool;
48 let mut a: f64;
49 let mut b: f64;
50 let mut temp: f64;
51
52 ix = get_high_word(x);
53 lx = get_low_word(x);
54 sign = (ix >> 31) != 0;
55 ix &= 0x7fffffff;
56
57 // -lx == !lx + 1
58 if (ix | (lx | ((!lx).wrapping_add(1))) >> 31) > 0x7ff00000 {
59 /* nan */
60 return x;
61 }
62
63 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
64 * Thus, J(-n,x) = J(n,-x)
65 */
66 /* nm1 = |n|-1 is used instead of |n| to handle n==INT_MIN */
67 if n == 0 {
68 return j0(x);
69 }
70 if n < 0 {
71 nm1 = -(n + 1);
72 x = -x;
73 sign = !sign;
74 } else {
75 nm1 = n - 1;
76 }
77 if nm1 == 0 {
78 return j1(x);
79 }
80
81 sign &= (n & 1) != 0; /* even n: 0, odd n: signbit(x) */
82 x = fabs(x);
83 if (ix | lx) == 0 || ix == 0x7ff00000 {
84 /* if x is 0 or inf */
85 b = 0.0;
86 } else if (nm1 as f64) < x {
87 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
88 if ix >= 0x52d00000 {
89 /* x > 2**302 */
90 /* (x >> n**2)
91 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
92 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
93 * Let s=sin(x), c=cos(x),
94 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
95 *
96 * n sin(xn)*sqt2 cos(xn)*sqt2
97 * ----------------------------------
98 * 0 s-c c+s
99 * 1 -s-c -c+s
100 * 2 -s+c -c-s
101 * 3 s+c c-s
102 */
103 temp = match nm1 & 3 {
104 0 => -cos(x) + sin(x),
105 1 => -cos(x) - sin(x),
106 2 => cos(x) - sin(x),
107 3 | _ => cos(x) + sin(x),
108 };
109 b = INVSQRTPI * temp / sqrt(x);
110 } else {
111 a = j0(x);
112 b = j1(x);
113 i = 0;
114 while i < nm1 {
115 i += 1;
116 temp = b;
117 b = b * (2.0 * (i as f64) / x) - a; /* avoid underflow */
118 a = temp;
119 }
120 }
121 } else {
122 if ix < 0x3e100000 {
123 /* x < 2**-29 */
124 /* x is tiny, return the first Taylor expansion of J(n,x)
125 * J(n,x) = 1/n!*(x/2)^n - ...
126 */
127 if nm1 > 32 {
128 /* underflow */
129 b = 0.0;
130 } else {
131 temp = x * 0.5;
132 b = temp;
133 a = 1.0;
134 i = 2;
135 while i <= nm1 + 1 {
136 a *= i as f64; /* a = n! */
137 b *= temp; /* b = (x/2)^n */
138 i += 1;
139 }
140 b = b / a;
141 }
142 } else {
143 /* use backward recurrence */
144 /* x x^2 x^2
145 * J(n,x)/J(n-1,x) = ---- ------ ------ .....
146 * 2n - 2(n+1) - 2(n+2)
147 *
148 * 1 1 1
149 * (for large x) = ---- ------ ------ .....
150 * 2n 2(n+1) 2(n+2)
151 * -- - ------ - ------ -
152 * x x x
153 *
154 * Let w = 2n/x and h=2/x, then the above quotient
155 * is equal to the continued fraction:
156 * 1
157 * = -----------------------
158 * 1
159 * w - -----------------
160 * 1
161 * w+h - ---------
162 * w+2h - ...
163 *
164 * To determine how many terms needed, let
165 * Q(0) = w, Q(1) = w(w+h) - 1,
166 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
167 * When Q(k) > 1e4 good for single
168 * When Q(k) > 1e9 good for double
169 * When Q(k) > 1e17 good for quadruple
170 */
171 /* determine k */
172 let mut t: f64;
173 let mut q0: f64;
174 let mut q1: f64;
175 let mut w: f64;
176 let h: f64;
177 let mut z: f64;
178 let mut tmp: f64;
179 let nf: f64;
180
181 let mut k: i32;
182
183 nf = (nm1 as f64) + 1.0;
184 w = 2.0 * nf / x;
185 h = 2.0 / x;
186 z = w + h;
187 q0 = w;
188 q1 = w * z - 1.0;
189 k = 1;
190 while q1 < 1.0e9 {
191 k += 1;
192 z += h;
193 tmp = z * q1 - q0;
194 q0 = q1;
195 q1 = tmp;
196 }
197 t = 0.0;
198 i = k;
199 while i >= 0 {
200 t = 1.0 / (2.0 * ((i as f64) + nf) / x - t);
201 i -= 1;
202 }
203 a = t;
204 b = 1.0;
205 /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
206 * Hence, if n*(log(2n/x)) > ...
207 * single 8.8722839355e+01
208 * double 7.09782712893383973096e+02
209 * long double 1.1356523406294143949491931077970765006170e+04
210 * then recurrent value may overflow and the result is
211 * likely underflow to zero
212 */
213 tmp = nf * log(fabs(w));
214 if tmp < 7.09782712893383973096e+02 {
215 i = nm1;
216 while i > 0 {
217 temp = b;
218 b = b * (2.0 * (i as f64)) / x - a;
219 a = temp;
220 i -= 1;
221 }
222 } else {
223 i = nm1;
224 while i > 0 {
225 temp = b;
226 b = b * (2.0 * (i as f64)) / x - a;
227 a = temp;
228 /* scale b to avoid spurious overflow */
229 let x1p500 = f64::from_bits(0x5f30000000000000); // 0x1p500 == 2^500
230 if b > x1p500 {
231 a /= b;
232 t /= b;
233 b = 1.0;
234 }
235 i -= 1;
236 }
237 }
238 z = j0(x);
239 w = j1(x);
240 if fabs(z) >= fabs(w) {
241 b = t * z / b;
242 } else {
243 b = t * w / a;
244 }
245 }
246 }
247
248 if sign { -b } else { b }
249}
250
251/// Integer order of the [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the second kind (f64).
252pub fn yn(n: i32, x: f64) -> f64 {
253 let mut ix: u32;
254 let lx: u32;
255 let mut ib: u32;
256 let nm1: i32;
257 let mut sign: bool;
258 let mut i: i32;
259 let mut a: f64;
260 let mut b: f64;
261 let mut temp: f64;
262
263 ix = get_high_word(x);
264 lx = get_low_word(x);
265 sign = (ix >> 31) != 0;
266 ix &= 0x7fffffff;
267
268 // -lx == !lx + 1
269 if (ix | (lx | ((!lx).wrapping_add(1))) >> 31) > 0x7ff00000 {
270 /* nan */
271 return x;
272 }
273 if sign && (ix | lx) != 0 {
274 /* x < 0 */
275 return 0.0 / 0.0;
276 }
277 if ix == 0x7ff00000 {
278 return 0.0;
279 }
280
281 if n == 0 {
282 return y0(x);
283 }
284 if n < 0 {
285 nm1 = -(n + 1);
286 sign = (n & 1) != 0;
287 } else {
288 nm1 = n - 1;
289 sign = false;
290 }
291 if nm1 == 0 {
292 if sign {
293 return -y1(x);
294 } else {
295 return y1(x);
296 }
297 }
298
299 if ix >= 0x52d00000 {
300 /* x > 2**302 */
301 /* (x >> n**2)
302 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
303 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
304 * Let s=sin(x), c=cos(x),
305 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
306 *
307 * n sin(xn)*sqt2 cos(xn)*sqt2
308 * ----------------------------------
309 * 0 s-c c+s
310 * 1 -s-c -c+s
311 * 2 -s+c -c-s
312 * 3 s+c c-s
313 */
314 temp = match nm1 & 3 {
315 0 => -sin(x) - cos(x),
316 1 => -sin(x) + cos(x),
317 2 => sin(x) + cos(x),
318 3 | _ => sin(x) - cos(x),
319 };
320 b = INVSQRTPI * temp / sqrt(x);
321 } else {
322 a = y0(x);
323 b = y1(x);
324 /* quit if b is -inf */
325 ib = get_high_word(b);
326 i = 0;
327 while i < nm1 && ib != 0xfff00000 {
328 i += 1;
329 temp = b;
330 b = (2.0 * (i as f64) / x) * b - a;
331 ib = get_high_word(b);
332 a = temp;
333 }
334 }
335
336 if sign { -b } else { b }
337}
338