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 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 overflow 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 long double
67 invsqrtpi = 5.64189583547756286948079e-1L, two = 2.0e0L, one = 1.0e0L;
68
69static const long double zero = 0.0L;
70
71long double
72__ieee754_jnl (int n, long double x)
73{
74 uint32_t se, i0, i1;
75 int32_t i, ix, sgn;
76 long double a, b, temp, di, ret;
77 long double z, w;
78
79 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
80 * Thus, J(-n,x) = J(n,-x)
81 */
82
83 GET_LDOUBLE_WORDS (se, i0, i1, x);
84 ix = se & 0x7fff;
85
86 /* if J(n,NaN) is NaN */
87 if (__glibc_unlikely ((ix == 0x7fff) && ((i0 & 0x7fffffff) != 0)))
88 return x + x;
89 if (n < 0)
90 {
91 n = -n;
92 x = -x;
93 se ^= 0x8000;
94 }
95 if (n == 0)
96 return (__ieee754_j0l (x));
97 if (n == 1)
98 return (__ieee754_j1l (x));
99 sgn = (n & 1) & (se >> 15); /* even n -- 0, odd n -- sign(x) */
100 x = fabsl (x);
101 {
102 SET_RESTORE_ROUNDL (FE_TONEAREST);
103 if (__glibc_unlikely ((ix | i0 | i1) == 0 || ix >= 0x7fff))
104 /* if x is 0 or inf */
105 return sgn == 1 ? -zero : zero;
106 else if ((long double) n <= x)
107 {
108 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
109 if (ix >= 0x412D)
110 { /* x > 2**302 */
111
112 /* ??? This might be a futile gesture.
113 If x exceeds X_TLOSS anyway, the wrapper function
114 will set the result to zero. */
115
116 /* (x >> n**2)
117 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
118 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
119 * Let s=sin(x), c=cos(x),
120 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
121 *
122 * n sin(xn)*sqt2 cos(xn)*sqt2
123 * ----------------------------------
124 * 0 s-c c+s
125 * 1 -s-c -c+s
126 * 2 -s+c -c-s
127 * 3 s+c c-s
128 */
129 long double s;
130 long double c;
131 __sincosl (x, &s, &c);
132 switch (n & 3)
133 {
134 case 0:
135 temp = c + s;
136 break;
137 case 1:
138 temp = -c + s;
139 break;
140 case 2:
141 temp = -c - s;
142 break;
143 case 3:
144 temp = c - s;
145 break;
146 default:
147 __builtin_unreachable ();
148 }
149 b = invsqrtpi * temp / sqrtl (x);
150 }
151 else
152 {
153 a = __ieee754_j0l (x);
154 b = __ieee754_j1l (x);
155 for (i = 1; i < n; i++)
156 {
157 temp = b;
158 b = b * ((long double) (i + i) / x) - a; /* avoid underflow */
159 a = temp;
160 }
161 }
162 }
163 else
164 {
165 if (ix < 0x3fde)
166 { /* x < 2**-33 */
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 *= (long double) 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 long double t, v;
216 long double q0, q1, h, tmp;
217 int32_t k, m;
218 w = (n + n) / (long double) x;
219 h = 2.0L / (long double) x;
220 q0 = w;
221 z = w + h;
222 q1 = w * z - 1.0L;
223 k = 1;
224 while (q1 < 1.0e11L)
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 * __ieee754_logl (fabsl (v * tmp));
248
249 if (tmp < 1.1356523406294143949491931077970765006170e+04L)
250 {
251 for (i = n - 1, di = (long double) (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 = (long double) (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 > 1e100L)
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 = __ieee754_j0l (x);
283 w = __ieee754_j1l (x);
284 if (fabsl (z) >= fabsl (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 = copysignl (LDBL_MIN, ret) * LDBL_MIN;
298 __set_errno (ERANGE);
299 }
300 else
301 math_check_force_underflow (ret);
302 return ret;
303}
304strong_alias (__ieee754_jnl, __jnl_finite)
305
306long double
307__ieee754_ynl (int n, long double x)
308{
309 uint32_t se, i0, i1;
310 int32_t i, ix;
311 int32_t sign;
312 long double a, b, temp, ret;
313
314
315 GET_LDOUBLE_WORDS (se, i0, i1, x);
316 ix = se & 0x7fff;
317 /* if Y(n,NaN) is NaN */
318 if (__builtin_expect ((ix == 0x7fff) && ((i0 & 0x7fffffff) != 0), 0))
319 return x + x;
320 if (__builtin_expect ((ix | i0 | i1) == 0, 0))
321 /* -inf or inf and divide-by-zero exception. */
322 return ((n < 0 && (n & 1) != 0) ? 1.0L : -1.0L) / 0.0L;
323 if (__builtin_expect (se & 0x8000, 0))
324 return zero / (zero * x);
325 sign = 1;
326 if (n < 0)
327 {
328 n = -n;
329 sign = 1 - ((n & 1) << 1);
330 }
331 if (n == 0)
332 return (__ieee754_y0l (x));
333 {
334 SET_RESTORE_ROUNDL (FE_TONEAREST);
335 if (n == 1)
336 {
337 ret = sign * __ieee754_y1l (x);
338 goto out;
339 }
340 if (__glibc_unlikely (ix == 0x7fff))
341 return zero;
342 if (ix >= 0x412D)
343 { /* x > 2**302 */
344
345 /* ??? See comment above on the possible futility of this. */
346
347 /* (x >> n**2)
348 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
349 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
350 * Let s=sin(x), c=cos(x),
351 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
352 *
353 * n sin(xn)*sqt2 cos(xn)*sqt2
354 * ----------------------------------
355 * 0 s-c c+s
356 * 1 -s-c -c+s
357 * 2 -s+c -c-s
358 * 3 s+c c-s
359 */
360 long double s;
361 long double c;
362 __sincosl (x, &s, &c);
363 switch (n & 3)
364 {
365 case 0:
366 temp = s - c;
367 break;
368 case 1:
369 temp = -s - c;
370 break;
371 case 2:
372 temp = -s + c;
373 break;
374 case 3:
375 temp = s + c;
376 break;
377 default:
378 __builtin_unreachable ();
379 }
380 b = invsqrtpi * temp / sqrtl (x);
381 }
382 else
383 {
384 a = __ieee754_y0l (x);
385 b = __ieee754_y1l (x);
386 /* quit if b is -inf */
387 GET_LDOUBLE_WORDS (se, i0, i1, b);
388 /* Use 0xffffffff since GET_LDOUBLE_WORDS sign-extends SE. */
389 for (i = 1; i < n && se != 0xffffffff; i++)
390 {
391 temp = b;
392 b = ((long double) (i + i) / x) * b - a;
393 GET_LDOUBLE_WORDS (se, i0, i1, b);
394 a = temp;
395 }
396 }
397 /* If B is +-Inf, set up errno accordingly. */
398 if (! isfinite (b))
399 __set_errno (ERANGE);
400 if (sign > 0)
401 ret = b;
402 else
403 ret = -b;
404 }
405 out:
406 if (isinf (ret))
407 ret = copysignl (LDBL_MAX, ret) * LDBL_MAX;
408 return ret;
409}
410strong_alias (__ieee754_ynl, __ynl_finite)
411