compiler_builtins/math/libm_math/jnf.rs
1/* origin: FreeBSD /usr/src/lib/msun/src/e_jnf.c */
2/*
3 * Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
4 */
5/*
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 *
9 * Developed at SunPro, a Sun Microsystems, Inc. business.
10 * Permission to use, copy, modify, and distribute this
11 * software is freely granted, provided that this notice
12 * is preserved.
13 * ====================================================
14 */
15
16use super::{fabsf, j0f, j1f, logf, y0f, y1f};
17
18/// Integer order of the [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the first kind (f32).
19#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
20pub fn jnf(n: i32, mut x: f32) -> f32 {
21 let mut ix: u32;
22 let mut nm1: i32;
23 let mut sign: bool;
24 let mut i: i32;
25 let mut a: f32;
26 let mut b: f32;
27 let mut temp: f32;
28
29 ix = x.to_bits();
30 sign = (ix >> 31) != 0;
31 ix &= 0x7fffffff;
32 if ix > 0x7f800000 {
33 /* nan */
34 return x;
35 }
36
37 /* J(-n,x) = J(n,-x), use |n|-1 to avoid overflow in -n */
38 if n == 0 {
39 return j0f(x);
40 }
41 if n < 0 {
42 nm1 = -(n + 1);
43 x = -x;
44 sign = !sign;
45 } else {
46 nm1 = n - 1;
47 }
48 if nm1 == 0 {
49 return j1f(x);
50 }
51
52 sign &= (n & 1) != 0; /* even n: 0, odd n: signbit(x) */
53 x = fabsf(x);
54 if ix == 0 || ix == 0x7f800000 {
55 /* if x is 0 or inf */
56 b = 0.0;
57 } else if (nm1 as f32) < x {
58 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
59 a = j0f(x);
60 b = j1f(x);
61 i = 0;
62 while i < nm1 {
63 i += 1;
64 temp = b;
65 b = b * (2.0 * (i as f32) / x) - a;
66 a = temp;
67 }
68 } else if ix < 0x35800000 {
69 /* x < 2**-20 */
70 /* x is tiny, return the first Taylor expansion of J(n,x)
71 * J(n,x) = 1/n!*(x/2)^n - ...
72 */
73 if nm1 > 8 {
74 /* underflow */
75 nm1 = 8;
76 }
77 temp = 0.5 * x;
78 b = temp;
79 a = 1.0;
80 i = 2;
81 while i <= nm1 + 1 {
82 a *= i as f32; /* a = n! */
83 b *= temp; /* b = (x/2)^n */
84 i += 1;
85 }
86 b = b / a;
87 } else {
88 /* use backward recurrence */
89 /* x x^2 x^2
90 * J(n,x)/J(n-1,x) = ---- ------ ------ .....
91 * 2n - 2(n+1) - 2(n+2)
92 *
93 * 1 1 1
94 * (for large x) = ---- ------ ------ .....
95 * 2n 2(n+1) 2(n+2)
96 * -- - ------ - ------ -
97 * x x x
98 *
99 * Let w = 2n/x and h=2/x, then the above quotient
100 * is equal to the continued fraction:
101 * 1
102 * = -----------------------
103 * 1
104 * w - -----------------
105 * 1
106 * w+h - ---------
107 * w+2h - ...
108 *
109 * To determine how many terms needed, let
110 * Q(0) = w, Q(1) = w(w+h) - 1,
111 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
112 * When Q(k) > 1e4 good for single
113 * When Q(k) > 1e9 good for double
114 * When Q(k) > 1e17 good for quadruple
115 */
116 /* determine k */
117 let mut t: f32;
118 let mut q0: f32;
119 let mut q1: f32;
120 let mut w: f32;
121 let h: f32;
122 let mut z: f32;
123 let mut tmp: f32;
124 let nf: f32;
125 let mut k: i32;
126
127 nf = (nm1 as f32) + 1.0;
128 w = 2.0 * nf / x;
129 h = 2.0 / x;
130 z = w + h;
131 q0 = w;
132 q1 = w * z - 1.0;
133 k = 1;
134 while q1 < 1.0e4 {
135 k += 1;
136 z += h;
137 tmp = z * q1 - q0;
138 q0 = q1;
139 q1 = tmp;
140 }
141 t = 0.0;
142 i = k;
143 while i >= 0 {
144 t = 1.0 / (2.0 * ((i as f32) + nf) / x - t);
145 i -= 1;
146 }
147 a = t;
148 b = 1.0;
149 /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
150 * Hence, if n*(log(2n/x)) > ...
151 * single 8.8722839355e+01
152 * double 7.09782712893383973096e+02
153 * long double 1.1356523406294143949491931077970765006170e+04
154 * then recurrent value may overflow and the result is
155 * likely underflow to zero
156 */
157 tmp = nf * logf(fabsf(w));
158 if tmp < 88.721679688 {
159 i = nm1;
160 while i > 0 {
161 temp = b;
162 b = 2.0 * (i as f32) * b / x - a;
163 a = temp;
164 i -= 1;
165 }
166 } else {
167 i = nm1;
168 while i > 0 {
169 temp = b;
170 b = 2.0 * (i as f32) * b / x - a;
171 a = temp;
172 /* scale b to avoid spurious overflow */
173 let x1p60 = f32::from_bits(0x5d800000); // 0x1p60 == 2^60
174 if b > x1p60 {
175 a /= b;
176 t /= b;
177 b = 1.0;
178 }
179 i -= 1;
180 }
181 }
182 z = j0f(x);
183 w = j1f(x);
184 if fabsf(z) >= fabsf(w) {
185 b = t * z / b;
186 } else {
187 b = t * w / a;
188 }
189 }
190
191 if sign { -b } else { b }
192}
193
194/// Integer order of the [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the second kind (f32).
195#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
196pub fn ynf(n: i32, x: f32) -> f32 {
197 let mut ix: u32;
198 let mut ib: u32;
199 let nm1: i32;
200 let mut sign: bool;
201 let mut i: i32;
202 let mut a: f32;
203 let mut b: f32;
204 let mut temp: f32;
205
206 ix = x.to_bits();
207 sign = (ix >> 31) != 0;
208 ix &= 0x7fffffff;
209 if ix > 0x7f800000 {
210 /* nan */
211 return x;
212 }
213 if sign && ix != 0 {
214 /* x < 0 */
215 return 0.0 / 0.0;
216 }
217 if ix == 0x7f800000 {
218 return 0.0;
219 }
220
221 if n == 0 {
222 return y0f(x);
223 }
224 if n < 0 {
225 nm1 = -(n + 1);
226 sign = (n & 1) != 0;
227 } else {
228 nm1 = n - 1;
229 sign = false;
230 }
231 if nm1 == 0 {
232 if sign {
233 return -y1f(x);
234 } else {
235 return y1f(x);
236 }
237 }
238
239 a = y0f(x);
240 b = y1f(x);
241 /* quit if b is -inf */
242 ib = b.to_bits();
243 i = 0;
244 while i < nm1 && ib != 0xff800000 {
245 i += 1;
246 temp = b;
247 b = (2.0 * (i as f32) / x) * b - a;
248 ib = b.to_bits();
249 a = temp;
250 }
251
252 if sign { -b } else { b }
253}