1 /*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 */
11
12 /* Modifications for 128-bit long double are
13 Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
14 and are incorporated herein by permission of the author. The author
15 reserves the right to distribute this material elsewhere under different
16 copying permissions. These modifications are distributed here under
17 the following terms:
18
19 This library is free software; you can redistribute it and/or
20 modify it under the terms of the GNU Lesser General Public
21 License as published by the Free Software Foundation; either
22 version 2.1 of the License, or (at your option) any later version.
23
24 This library is distributed in the hope that it will be useful,
25 but WITHOUT ANY WARRANTY; without even the implied warranty of
26 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 Lesser General Public License for more details.
28
29 You should have received a copy of the GNU Lesser General Public
30 License along with this library; if not, see
31 <http://www.gnu.org/licenses/>. */
32
33 /*
34 * __ieee754_jn(n, x), __ieee754_yn(n, x)
35 * floating point Bessel's function of the 1st and 2nd kind
36 * of order n
37 *
38 * Special cases:
39 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
40 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
41 * Note 2. About jn(n,x), yn(n,x)
42 * For n=0, j0(x) is called,
43 * for n=1, j1(x) is called,
44 * for n<x, forward recursion us used starting
45 * from values of j0(x) and j1(x).
46 * for n>x, a continued fraction approximation to
47 * j(n,x)/j(n-1,x) is evaluated and then backward
48 * recursion is used starting from a supposed value
49 * for j(n,x). The resulting value of j(0,x) is
50 * compared with the actual value to correct the
51 * supposed value of j(n,x).
52 *
53 * yn(n,x) is similar in all respects, except
54 * that forward recursion is used for all
55 * values of n>1.
56 *
57 */
58
59 #include "quadmath-imp.h"
60
61 static const __float128
62 invsqrtpi = 5.6418958354775628694807945156077258584405E-1Q,
63 two = 2,
64 one = 1,
65 zero = 0;
66
67
68 __float128
jnq(int n,__float128 x)69 jnq (int n, __float128 x)
70 {
71 uint32_t se;
72 int32_t i, ix, sgn;
73 __float128 a, b, temp, di, ret;
74 __float128 z, w;
75 ieee854_float128 u;
76
77
78 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
79 * Thus, J(-n,x) = J(n,-x)
80 */
81
82 u.value = x;
83 se = u.words32.w0;
84 ix = se & 0x7fffffff;
85
86 /* if J(n,NaN) is NaN */
87 if (ix >= 0x7fff0000)
88 {
89 if ((u.words32.w0 & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3)
90 return x + x;
91 }
92
93 if (n < 0)
94 {
95 n = -n;
96 x = -x;
97 se ^= 0x80000000;
98 }
99 if (n == 0)
100 return (j0q (x));
101 if (n == 1)
102 return (j1q (x));
103 sgn = (n & 1) & (se >> 31); /* even n -- 0, odd n -- sign(x) */
104 x = fabsq (x);
105
106 {
107 SET_RESTORE_ROUNDF128 (FE_TONEAREST);
108 if (x == 0 || ix >= 0x7fff0000) /* if x is 0 or inf */
109 return sgn == 1 ? -zero : zero;
110 else if ((__float128) n <= x)
111 {
112 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
113 if (ix >= 0x412D0000)
114 { /* x > 2**302 */
115
116 /* ??? Could use an expansion for large x here. */
117
118 /* (x >> n**2)
119 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
120 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
121 * Let s=sin(x), c=cos(x),
122 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
123 *
124 * n sin(xn)*sqt2 cos(xn)*sqt2
125 * ----------------------------------
126 * 0 s-c c+s
127 * 1 -s-c -c+s
128 * 2 -s+c -c-s
129 * 3 s+c c-s
130 */
131 __float128 s;
132 __float128 c;
133 sincosq (x, &s, &c);
134 switch (n & 3)
135 {
136 case 0:
137 temp = c + s;
138 break;
139 case 1:
140 temp = -c + s;
141 break;
142 case 2:
143 temp = -c - s;
144 break;
145 case 3:
146 temp = c - s;
147 break;
148 }
149 b = invsqrtpi * temp / sqrtq (x);
150 }
151 else
152 {
153 a = j0q (x);
154 b = j1q (x);
155 for (i = 1; i < n; i++)
156 {
157 temp = b;
158 b = b * ((__float128) (i + i) / x) - a; /* avoid underflow */
159 a = temp;
160 }
161 }
162 }
163 else
164 {
165 if (ix < 0x3fc60000)
166 { /* x < 2**-57 */
167 /* x is tiny, return the first Taylor expansion of J(n,x)
168 * J(n,x) = 1/n!*(x/2)^n - ...
169 */
170 if (n >= 400) /* underflow, result < 10^-4952 */
171 b = zero;
172 else
173 {
174 temp = x * 0.5;
175 b = temp;
176 for (a = one, i = 2; i <= n; i++)
177 {
178 a *= (__float128) i; /* a = n! */
179 b *= temp; /* b = (x/2)^n */
180 }
181 b = b / a;
182 }
183 }
184 else
185 {
186 /* use backward recurrence */
187 /* x x^2 x^2
188 * J(n,x)/J(n-1,x) = ---- ------ ------ .....
189 * 2n - 2(n+1) - 2(n+2)
190 *
191 * 1 1 1
192 * (for large x) = ---- ------ ------ .....
193 * 2n 2(n+1) 2(n+2)
194 * -- - ------ - ------ -
195 * x x x
196 *
197 * Let w = 2n/x and h=2/x, then the above quotient
198 * is equal to the continued fraction:
199 * 1
200 * = -----------------------
201 * 1
202 * w - -----------------
203 * 1
204 * w+h - ---------
205 * w+2h - ...
206 *
207 * To determine how many terms needed, let
208 * Q(0) = w, Q(1) = w(w+h) - 1,
209 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
210 * When Q(k) > 1e4 good for single
211 * When Q(k) > 1e9 good for double
212 * When Q(k) > 1e17 good for quadruple
213 */
214 /* determine k */
215 __float128 t, v;
216 __float128 q0, q1, h, tmp;
217 int32_t k, m;
218 w = (n + n) / (__float128) x;
219 h = 2 / (__float128) x;
220 q0 = w;
221 z = w + h;
222 q1 = w * z - 1;
223 k = 1;
224 while (q1 < 1.0e17Q)
225 {
226 k += 1;
227 z += h;
228 tmp = z * q1 - q0;
229 q0 = q1;
230 q1 = tmp;
231 }
232 m = n + n;
233 for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
234 t = one / (i / x - t);
235 a = t;
236 b = one;
237 /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
238 * Hence, if n*(log(2n/x)) > ...
239 * single 8.8722839355e+01
240 * double 7.09782712893383973096e+02
241 * long double 1.1356523406294143949491931077970765006170e+04
242 * then recurrent value may overflow and the result is
243 * likely underflow to zero
244 */
245 tmp = n;
246 v = two / x;
247 tmp = tmp * logq (fabsq (v * tmp));
248
249 if (tmp < 1.1356523406294143949491931077970765006170e+04Q)
250 {
251 for (i = n - 1, di = (__float128) (i + i); i > 0; i--)
252 {
253 temp = b;
254 b *= di;
255 b = b / x - a;
256 a = temp;
257 di -= two;
258 }
259 }
260 else
261 {
262 for (i = n - 1, di = (__float128) (i + i); i > 0; i--)
263 {
264 temp = b;
265 b *= di;
266 b = b / x - a;
267 a = temp;
268 di -= two;
269 /* scale b to avoid spurious overflow */
270 if (b > 1e100Q)
271 {
272 a /= b;
273 t /= b;
274 b = one;
275 }
276 }
277 }
278 /* j0() and j1() suffer enormous loss of precision at and
279 * near zero; however, we know that their zero points never
280 * coincide, so just choose the one further away from zero.
281 */
282 z = j0q (x);
283 w = j1q (x);
284 if (fabsq (z) >= fabsq (w))
285 b = (t * z / b);
286 else
287 b = (t * w / a);
288 }
289 }
290 if (sgn == 1)
291 ret = -b;
292 else
293 ret = b;
294 }
295 if (ret == 0)
296 {
297 ret = copysignq (FLT128_MIN, ret) * FLT128_MIN;
298 errno = ERANGE;
299 }
300 else
301 math_check_force_underflow (ret);
302 return ret;
303 }
304
305
306 __float128
ynq(int n,__float128 x)307 ynq (int n, __float128 x)
308 {
309 uint32_t se;
310 int32_t i, ix;
311 int32_t sign;
312 __float128 a, b, temp, ret;
313 ieee854_float128 u;
314
315 u.value = x;
316 se = u.words32.w0;
317 ix = se & 0x7fffffff;
318
319 /* if Y(n,NaN) is NaN */
320 if (ix >= 0x7fff0000)
321 {
322 if ((u.words32.w0 & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3)
323 return x + x;
324 }
325 if (x <= 0)
326 {
327 if (x == 0)
328 return ((n < 0 && (n & 1) != 0) ? 1 : -1) / 0.0Q;
329 if (se & 0x80000000)
330 return zero / (zero * x);
331 }
332 sign = 1;
333 if (n < 0)
334 {
335 n = -n;
336 sign = 1 - ((n & 1) << 1);
337 }
338 if (n == 0)
339 return (y0q (x));
340 {
341 SET_RESTORE_ROUNDF128 (FE_TONEAREST);
342 if (n == 1)
343 {
344 ret = sign * y1q (x);
345 goto out;
346 }
347 if (ix >= 0x7fff0000)
348 return zero;
349 if (ix >= 0x412D0000)
350 { /* x > 2**302 */
351
352 /* ??? See comment above on the possible futility of this. */
353
354 /* (x >> n**2)
355 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
356 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
357 * Let s=sin(x), c=cos(x),
358 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
359 *
360 * n sin(xn)*sqt2 cos(xn)*sqt2
361 * ----------------------------------
362 * 0 s-c c+s
363 * 1 -s-c -c+s
364 * 2 -s+c -c-s
365 * 3 s+c c-s
366 */
367 __float128 s;
368 __float128 c;
369 sincosq (x, &s, &c);
370 switch (n & 3)
371 {
372 case 0:
373 temp = s - c;
374 break;
375 case 1:
376 temp = -s - c;
377 break;
378 case 2:
379 temp = -s + c;
380 break;
381 case 3:
382 temp = s + c;
383 break;
384 }
385 b = invsqrtpi * temp / sqrtq (x);
386 }
387 else
388 {
389 a = y0q (x);
390 b = y1q (x);
391 /* quit if b is -inf */
392 u.value = b;
393 se = u.words32.w0 & 0xffff0000;
394 for (i = 1; i < n && se != 0xffff0000; i++)
395 {
396 temp = b;
397 b = ((__float128) (i + i) / x) * b - a;
398 u.value = b;
399 se = u.words32.w0 & 0xffff0000;
400 a = temp;
401 }
402 }
403 /* If B is +-Inf, set up errno accordingly. */
404 if (! finiteq (b))
405 errno = ERANGE;
406 if (sign > 0)
407 ret = b;
408 else
409 ret = -b;
410 }
411 out:
412 if (isinfq (ret))
413 ret = copysignq (FLT128_MAX, ret) * FLT128_MAX;
414 return ret;
415 }
416