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 <errno.h>
60#include <float.h>
61#include <math.h>
62#include <math_private.h>
63#include <fenv_private.h>
64#include <math-underflow.h>
65
66static const _Float128
67 invsqrtpi = L(5.6418958354775628694807945156077258584405E-1),
68 two = 2,
69 one = 1,
70 zero = 0;
71
72
73_Float128
74__ieee754_jnl (int n, _Float128 x)
75{
76 uint32_t se;
77 int32_t i, ix, sgn;
78 _Float128 a, b, temp, di, ret;
79 _Float128 z, w;
80 ieee854_long_double_shape_type u;
81
82
83 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
84 * Thus, J(-n,x) = J(n,-x)
85 */
86
87 u.value = x;
88 se = u.parts32.w0;
89 ix = se & 0x7fffffff;
90
91 /* if J(n,NaN) is NaN */
92 if (ix >= 0x7fff0000)
93 {
94 if ((u.parts32.w0 & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3)
95 return x + x;
96 }
97
98 if (n < 0)
99 {
100 n = -n;
101 x = -x;
102 se ^= 0x80000000;
103 }
104 if (n == 0)
105 return (__ieee754_j0l (x));
106 if (n == 1)
107 return (__ieee754_j1l (x));
108 sgn = (n & 1) & (se >> 31); /* even n -- 0, odd n -- sign(x) */
109 x = fabsl (x);
110
111 {
112 SET_RESTORE_ROUNDL (FE_TONEAREST);
113 if (x == 0 || ix >= 0x7fff0000) /* if x is 0 or inf */
114 return sgn == 1 ? -zero : zero;
115 else if ((_Float128) n <= x)
116 {
117 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
118 if (ix >= 0x412D0000)
119 { /* x > 2**302 */
120
121 /* ??? Could use an expansion for large x here. */
122
123 /* (x >> n**2)
124 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
125 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
126 * Let s=sin(x), c=cos(x),
127 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
128 *
129 * n sin(xn)*sqt2 cos(xn)*sqt2
130 * ----------------------------------
131 * 0 s-c c+s
132 * 1 -s-c -c+s
133 * 2 -s+c -c-s
134 * 3 s+c c-s
135 */
136 _Float128 s;
137 _Float128 c;
138 __sincosl (x, &s, &c);
139 switch (n & 3)
140 {
141 case 0:
142 temp = c + s;
143 break;
144 case 1:
145 temp = -c + s;
146 break;
147 case 2:
148 temp = -c - s;
149 break;
150 case 3:
151 temp = c - s;
152 break;
153 default:
154 __builtin_unreachable ();
155 }
156 b = invsqrtpi * temp / sqrtl (x);
157 }
158 else
159 {
160 a = __ieee754_j0l (x);
161 b = __ieee754_j1l (x);
162 for (i = 1; i < n; i++)
163 {
164 temp = b;
165 b = b * ((_Float128) (i + i) / x) - a; /* avoid underflow */
166 a = temp;
167 }
168 }
169 }
170 else
171 {
172 if (ix < 0x3fc60000)
173 { /* x < 2**-57 */
174 /* x is tiny, return the first Taylor expansion of J(n,x)
175 * J(n,x) = 1/n!*(x/2)^n - ...
176 */
177 if (n >= 400) /* underflow, result < 10^-4952 */
178 b = zero;
179 else
180 {
181 temp = x * 0.5;
182 b = temp;
183 for (a = one, i = 2; i <= n; i++)
184 {
185 a *= (_Float128) i; /* a = n! */
186 b *= temp; /* b = (x/2)^n */
187 }
188 b = b / a;
189 }
190 }
191 else
192 {
193 /* use backward recurrence */
194 /* x x^2 x^2
195 * J(n,x)/J(n-1,x) = ---- ------ ------ .....
196 * 2n - 2(n+1) - 2(n+2)
197 *
198 * 1 1 1
199 * (for large x) = ---- ------ ------ .....
200 * 2n 2(n+1) 2(n+2)
201 * -- - ------ - ------ -
202 * x x x
203 *
204 * Let w = 2n/x and h=2/x, then the above quotient
205 * is equal to the continued fraction:
206 * 1
207 * = -----------------------
208 * 1
209 * w - -----------------
210 * 1
211 * w+h - ---------
212 * w+2h - ...
213 *
214 * To determine how many terms needed, let
215 * Q(0) = w, Q(1) = w(w+h) - 1,
216 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
217 * When Q(k) > 1e4 good for single
218 * When Q(k) > 1e9 good for double
219 * When Q(k) > 1e17 good for quadruple
220 */
221 /* determine k */
222 _Float128 t, v;
223 _Float128 q0, q1, h, tmp;
224 int32_t k, m;
225 w = (n + n) / (_Float128) x;
226 h = 2 / (_Float128) x;
227 q0 = w;
228 z = w + h;
229 q1 = w * z - 1;
230 k = 1;
231 while (q1 < L(1.0e17))
232 {
233 k += 1;
234 z += h;
235 tmp = z * q1 - q0;
236 q0 = q1;
237 q1 = tmp;
238 }
239 m = n + n;
240 for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
241 t = one / (i / x - t);
242 a = t;
243 b = one;
244 /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
245 * Hence, if n*(log(2n/x)) > ...
246 * single 8.8722839355e+01
247 * double 7.09782712893383973096e+02
248 * long double 1.1356523406294143949491931077970765006170e+04
249 * then recurrent value may overflow and the result is
250 * likely underflow to zero
251 */
252 tmp = n;
253 v = two / x;
254 tmp = tmp * __ieee754_logl (fabsl (v * tmp));
255
256 if (tmp < L(1.1356523406294143949491931077970765006170e+04))
257 {
258 for (i = n - 1, di = (_Float128) (i + i); i > 0; i--)
259 {
260 temp = b;
261 b *= di;
262 b = b / x - a;
263 a = temp;
264 di -= two;
265 }
266 }
267 else
268 {
269 for (i = n - 1, di = (_Float128) (i + i); i > 0; i--)
270 {
271 temp = b;
272 b *= di;
273 b = b / x - a;
274 a = temp;
275 di -= two;
276 /* scale b to avoid spurious overflow */
277 if (b > L(1e100))
278 {
279 a /= b;
280 t /= b;
281 b = one;
282 }
283 }
284 }
285 /* j0() and j1() suffer enormous loss of precision at and
286 * near zero; however, we know that their zero points never
287 * coincide, so just choose the one further away from zero.
288 */
289 z = __ieee754_j0l (x);
290 w = __ieee754_j1l (x);
291 if (fabsl (z) >= fabsl (w))
292 b = (t * z / b);
293 else
294 b = (t * w / a);
295 }
296 }
297 if (sgn == 1)
298 ret = -b;
299 else
300 ret = b;
301 }
302 if (ret == 0)
303 {
304 ret = copysignl (LDBL_MIN, ret) * LDBL_MIN;
305 __set_errno (ERANGE);
306 }
307 else
308 math_check_force_underflow (ret);
309 return ret;
310}
311strong_alias (__ieee754_jnl, __jnl_finite)
312
313_Float128
314__ieee754_ynl (int n, _Float128 x)
315{
316 uint32_t se;
317 int32_t i, ix;
318 int32_t sign;
319 _Float128 a, b, temp, ret;
320 ieee854_long_double_shape_type u;
321
322 u.value = x;
323 se = u.parts32.w0;
324 ix = se & 0x7fffffff;
325
326 /* if Y(n,NaN) is NaN */
327 if (ix >= 0x7fff0000)
328 {
329 if ((u.parts32.w0 & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3)
330 return x + x;
331 }
332 if (x <= 0)
333 {
334 if (x == 0)
335 return ((n < 0 && (n & 1) != 0) ? 1 : -1) / L(0.0);
336 if (se & 0x80000000)
337 return zero / (zero * x);
338 }
339 sign = 1;
340 if (n < 0)
341 {
342 n = -n;
343 sign = 1 - ((n & 1) << 1);
344 }
345 if (n == 0)
346 return (__ieee754_y0l (x));
347 {
348 SET_RESTORE_ROUNDL (FE_TONEAREST);
349 if (n == 1)
350 {
351 ret = sign * __ieee754_y1l (x);
352 goto out;
353 }
354 if (ix >= 0x7fff0000)
355 return zero;
356 if (ix >= 0x412D0000)
357 { /* x > 2**302 */
358
359 /* ??? See comment above on the possible futility of this. */
360
361 /* (x >> n**2)
362 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
363 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
364 * Let s=sin(x), c=cos(x),
365 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
366 *
367 * n sin(xn)*sqt2 cos(xn)*sqt2
368 * ----------------------------------
369 * 0 s-c c+s
370 * 1 -s-c -c+s
371 * 2 -s+c -c-s
372 * 3 s+c c-s
373 */
374 _Float128 s;
375 _Float128 c;
376 __sincosl (x, &s, &c);
377 switch (n & 3)
378 {
379 case 0:
380 temp = s - c;
381 break;
382 case 1:
383 temp = -s - c;
384 break;
385 case 2:
386 temp = -s + c;
387 break;
388 case 3:
389 temp = s + c;
390 break;
391 default:
392 __builtin_unreachable ();
393 }
394 b = invsqrtpi * temp / sqrtl (x);
395 }
396 else
397 {
398 a = __ieee754_y0l (x);
399 b = __ieee754_y1l (x);
400 /* quit if b is -inf */
401 u.value = b;
402 se = u.parts32.w0 & 0xffff0000;
403 for (i = 1; i < n && se != 0xffff0000; i++)
404 {
405 temp = b;
406 b = ((_Float128) (i + i) / x) * b - a;
407 u.value = b;
408 se = u.parts32.w0 & 0xffff0000;
409 a = temp;
410 }
411 }
412 /* If B is +-Inf, set up errno accordingly. */
413 if (! isfinite (b))
414 __set_errno (ERANGE);
415 if (sign > 0)
416 ret = b;
417 else
418 ret = -b;
419 }
420 out:
421 if (isinf (ret))
422 ret = copysignl (LDBL_MAX, ret) * LDBL_MAX;
423 return ret;
424}
425strong_alias (__ieee754_ynl, __ynl_finite)
426