1/* @(#)e_jn.c 5.1 93/09/24 */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/*
14 * __ieee754_jn(n, x), __ieee754_yn(n, x)
15 * floating point Bessel's function of the 1st and 2nd kind
16 * of order n
17 *
18 * Special cases:
19 * y0(0)=y1(0)=yn(n,0) = -inf with overflow signal;
20 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
21 * Note 2. About jn(n,x), yn(n,x)
22 * For n=0, j0(x) is called,
23 * for n=1, j1(x) is called,
24 * for n<x, forward recursion us used starting
25 * from values of j0(x) and j1(x).
26 * for n>x, a continued fraction approximation to
27 * j(n,x)/j(n-1,x) is evaluated and then backward
28 * recursion is used starting from a supposed value
29 * for j(n,x). The resulting value of j(0,x) is
30 * compared with the actual value to correct the
31 * supposed value of j(n,x).
32 *
33 * yn(n,x) is similar in all respects, except
34 * that forward recursion is used for all
35 * values of n>1.
36 *
37 */
38
39#include <errno.h>
40#include <float.h>
41#include <math.h>
42#include <math-narrow-eval.h>
43#include <math_private.h>
44#include <fenv_private.h>
45#include <math-underflow.h>
46
47static const double
48 invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
49 two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
50 one = 1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
51
52static const double zero = 0.00000000000000000000e+00;
53
54double
55__ieee754_jn (int n, double x)
56{
57 int32_t i, hx, ix, lx, sgn;
58 double a, b, temp, di, ret;
59 double z, w;
60
61 /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
62 * Thus, J(-n,x) = J(n,-x)
63 */
64 EXTRACT_WORDS (hx, lx, x);
65 ix = 0x7fffffff & hx;
66 /* if J(n,NaN) is NaN */
67 if (__glibc_unlikely ((ix | ((uint32_t) (lx | -lx)) >> 31) > 0x7ff00000))
68 return x + x;
69 if (n < 0)
70 {
71 n = -n;
72 x = -x;
73 hx ^= 0x80000000;
74 }
75 if (n == 0)
76 return (__ieee754_j0 (x));
77 if (n == 1)
78 return (__ieee754_j1 (x));
79 sgn = (n & 1) & (hx >> 31); /* even n -- 0, odd n -- sign(x) */
80 x = fabs (x);
81 {
82 SET_RESTORE_ROUND (FE_TONEAREST);
83 if (__glibc_unlikely ((ix | lx) == 0 || ix >= 0x7ff00000))
84 /* if x is 0 or inf */
85 return sgn == 1 ? -zero : zero;
86 else if ((double) n <= x)
87 {
88 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
89 if (ix >= 0x52D00000) /* x > 2**302 */
90 { /* (x >> n**2)
91 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
92 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
93 * Let s=sin(x), c=cos(x),
94 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
95 *
96 * n sin(xn)*sqt2 cos(xn)*sqt2
97 * ----------------------------------
98 * 0 s-c c+s
99 * 1 -s-c -c+s
100 * 2 -s+c -c-s
101 * 3 s+c c-s
102 */
103 double s;
104 double c;
105 __sincos (x, &s, &c);
106 switch (n & 3)
107 {
108 case 0: temp = c + s; break;
109 case 1: temp = -c + s; break;
110 case 2: temp = -c - s; break;
111 case 3: temp = c - s; break;
112 default: __builtin_unreachable ();
113 }
114 b = invsqrtpi * temp / sqrt (x);
115 }
116 else
117 {
118 a = __ieee754_j0 (x);
119 b = __ieee754_j1 (x);
120 for (i = 1; i < n; i++)
121 {
122 temp = b;
123 b = b * ((double) (i + i) / x) - a; /* avoid underflow */
124 a = temp;
125 }
126 }
127 }
128 else
129 {
130 if (ix < 0x3e100000) /* x < 2**-29 */
131 { /* x is tiny, return the first Taylor expansion of J(n,x)
132 * J(n,x) = 1/n!*(x/2)^n - ...
133 */
134 if (n > 33) /* underflow */
135 b = zero;
136 else
137 {
138 temp = x * 0.5; b = temp;
139 for (a = one, i = 2; i <= n; i++)
140 {
141 a *= (double) i; /* a = n! */
142 b *= temp; /* b = (x/2)^n */
143 }
144 b = b / a;
145 }
146 }
147 else
148 {
149 /* use backward recurrence */
150 /* x x^2 x^2
151 * J(n,x)/J(n-1,x) = ---- ------ ------ .....
152 * 2n - 2(n+1) - 2(n+2)
153 *
154 * 1 1 1
155 * (for large x) = ---- ------ ------ .....
156 * 2n 2(n+1) 2(n+2)
157 * -- - ------ - ------ -
158 * x x x
159 *
160 * Let w = 2n/x and h=2/x, then the above quotient
161 * is equal to the continued fraction:
162 * 1
163 * = -----------------------
164 * 1
165 * w - -----------------
166 * 1
167 * w+h - ---------
168 * w+2h - ...
169 *
170 * To determine how many terms needed, let
171 * Q(0) = w, Q(1) = w(w+h) - 1,
172 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
173 * When Q(k) > 1e4 good for single
174 * When Q(k) > 1e9 good for double
175 * When Q(k) > 1e17 good for quadruple
176 */
177 /* determine k */
178 double t, v;
179 double q0, q1, h, tmp; int32_t k, m;
180 w = (n + n) / (double) x; h = 2.0 / (double) x;
181 q0 = w; z = w + h; q1 = w * z - 1.0; k = 1;
182 while (q1 < 1.0e9)
183 {
184 k += 1; z += h;
185 tmp = z * q1 - q0;
186 q0 = q1;
187 q1 = tmp;
188 }
189 m = n + n;
190 for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
191 t = one / (i / x - t);
192 a = t;
193 b = one;
194 /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
195 * Hence, if n*(log(2n/x)) > ...
196 * single 8.8722839355e+01
197 * double 7.09782712893383973096e+02
198 * long double 1.1356523406294143949491931077970765006170e+04
199 * then recurrent value may overflow and the result is
200 * likely underflow to zero
201 */
202 tmp = n;
203 v = two / x;
204 tmp = tmp * __ieee754_log (fabs (v * tmp));
205 if (tmp < 7.09782712893383973096e+02)
206 {
207 for (i = n - 1, di = (double) (i + i); i > 0; i--)
208 {
209 temp = b;
210 b *= di;
211 b = b / x - a;
212 a = temp;
213 di -= two;
214 }
215 }
216 else
217 {
218 for (i = n - 1, di = (double) (i + i); i > 0; i--)
219 {
220 temp = b;
221 b *= di;
222 b = b / x - a;
223 a = temp;
224 di -= two;
225 /* scale b to avoid spurious overflow */
226 if (b > 1e100)
227 {
228 a /= b;
229 t /= b;
230 b = one;
231 }
232 }
233 }
234 /* j0() and j1() suffer enormous loss of precision at and
235 * near zero; however, we know that their zero points never
236 * coincide, so just choose the one further away from zero.
237 */
238 z = __ieee754_j0 (x);
239 w = __ieee754_j1 (x);
240 if (fabs (z) >= fabs (w))
241 b = (t * z / b);
242 else
243 b = (t * w / a);
244 }
245 }
246 if (sgn == 1)
247 ret = -b;
248 else
249 ret = b;
250 ret = math_narrow_eval (ret);
251 }
252 if (ret == 0)
253 {
254 ret = math_narrow_eval (copysign (DBL_MIN, ret) * DBL_MIN);
255 __set_errno (ERANGE);
256 }
257 else
258 math_check_force_underflow (ret);
259 return ret;
260}
261strong_alias (__ieee754_jn, __jn_finite)
262
263double
264__ieee754_yn (int n, double x)
265{
266 int32_t i, hx, ix, lx;
267 int32_t sign;
268 double a, b, temp, ret;
269
270 EXTRACT_WORDS (hx, lx, x);
271 ix = 0x7fffffff & hx;
272 /* if Y(n,NaN) is NaN */
273 if (__glibc_unlikely ((ix | ((uint32_t) (lx | -lx)) >> 31) > 0x7ff00000))
274 return x + x;
275 sign = 1;
276 if (n < 0)
277 {
278 n = -n;
279 sign = 1 - ((n & 1) << 1);
280 }
281 if (n == 0)
282 return (__ieee754_y0 (x));
283 if (__glibc_unlikely ((ix | lx) == 0))
284 return -sign / zero;
285 /* -inf and overflow exception. */;
286 if (__glibc_unlikely (hx < 0))
287 return zero / (zero * x);
288 {
289 SET_RESTORE_ROUND (FE_TONEAREST);
290 if (n == 1)
291 {
292 ret = sign * __ieee754_y1 (x);
293 goto out;
294 }
295 if (__glibc_unlikely (ix == 0x7ff00000))
296 return zero;
297 if (ix >= 0x52D00000) /* x > 2**302 */
298 { /* (x >> n**2)
299 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
300 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
301 * Let s=sin(x), c=cos(x),
302 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
303 *
304 * n sin(xn)*sqt2 cos(xn)*sqt2
305 * ----------------------------------
306 * 0 s-c c+s
307 * 1 -s-c -c+s
308 * 2 -s+c -c-s
309 * 3 s+c c-s
310 */
311 double c;
312 double s;
313 __sincos (x, &s, &c);
314 switch (n & 3)
315 {
316 case 0: temp = s - c; break;
317 case 1: temp = -s - c; break;
318 case 2: temp = -s + c; break;
319 case 3: temp = s + c; break;
320 default: __builtin_unreachable ();
321 }
322 b = invsqrtpi * temp / sqrt (x);
323 }
324 else
325 {
326 uint32_t high;
327 a = __ieee754_y0 (x);
328 b = __ieee754_y1 (x);
329 /* quit if b is -inf */
330 GET_HIGH_WORD (high, b);
331 for (i = 1; i < n && high != 0xfff00000; i++)
332 {
333 temp = b;
334 b = ((double) (i + i) / x) * b - a;
335 GET_HIGH_WORD (high, b);
336 a = temp;
337 }
338 /* If B is +-Inf, set up errno accordingly. */
339 if (!isfinite (b))
340 __set_errno (ERANGE);
341 }
342 if (sign > 0)
343 ret = b;
344 else
345 ret = -b;
346 }
347 out:
348 if (isinf (ret))
349 ret = copysign (DBL_MAX, ret) * DBL_MAX;
350 return ret;
351}
352strong_alias (__ieee754_yn, __yn_finite)
353