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