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