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