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