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, write to the Free Software
31 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
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 "quadmath-imp.h"
61
62 static const __float128
63 invsqrtpi = 5.6418958354775628694807945156077258584405E-1Q,
64 two = 2.0e0Q,
65 one = 1.0e0Q,
66 zero = 0.0Q;
67
68
69 __float128
jnq(int n,__float128 x)70 jnq (int n, __float128 x)
71 {
72 uint32_t se;
73 int32_t i, ix, sgn;
74 __float128 a, b, temp, di;
75 __float128 z, w;
76 ieee854_float128 u;
77
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 u.value = x;
84 se = u.words32.w0;
85 ix = se & 0x7fffffff;
86
87 /* if J(n,NaN) is NaN */
88 if (ix >= 0x7fff0000)
89 {
90 if ((u.words32.w0 & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3)
91 return x + x;
92 }
93
94 if (n < 0)
95 {
96 n = -n;
97 x = -x;
98 se ^= 0x80000000;
99 }
100 if (n == 0)
101 return (j0q (x));
102 if (n == 1)
103 return (j1q (x));
104 sgn = (n & 1) & (se >> 31); /* even n -- 0, odd n -- sign(x) */
105 x = fabsq (x);
106
107 if (x == 0.0Q || ix >= 0x7fff0000) /* if x is 0 or inf */
108 b = zero;
109 else if ((__float128) n <= x)
110 {
111 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
112 if (ix >= 0x412D0000)
113 { /* x > 2**302 */
114
115 /* ??? Could use an expansion for large x here. */
116
117 /* (x >> n**2)
118 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
119 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
120 * Let s=sin(x), c=cos(x),
121 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
122 *
123 * n sin(xn)*sqt2 cos(xn)*sqt2
124 * ----------------------------------
125 * 0 s-c c+s
126 * 1 -s-c -c+s
127 * 2 -s+c -c-s
128 * 3 s+c c-s
129 */
130 __float128 s;
131 __float128 c;
132 sincosq (x, &s, &c);
133 switch (n & 3)
134 {
135 case 0:
136 temp = c + s;
137 break;
138 case 1:
139 temp = -c + s;
140 break;
141 case 2:
142 temp = -c - s;
143 break;
144 case 3:
145 temp = c - s;
146 break;
147 }
148 b = invsqrtpi * temp / sqrtq (x);
149 }
150 else
151 {
152 a = j0q (x);
153 b = j1q (x);
154 for (i = 1; i < n; i++)
155 {
156 temp = b;
157 b = b * ((__float128) (i + i) / x) - a; /* avoid underflow */
158 a = temp;
159 }
160 }
161 }
162 else
163 {
164 if (ix < 0x3fc60000)
165 { /* x < 2**-57 */
166 /* x is tiny, return the first Taylor expansion of J(n,x)
167 * J(n,x) = 1/n!*(x/2)^n - ...
168 */
169 if (n >= 400) /* underflow, result < 10^-4952 */
170 b = zero;
171 else
172 {
173 temp = x * 0.5;
174 b = temp;
175 for (a = one, i = 2; i <= n; i++)
176 {
177 a *= (__float128) i; /* a = n! */
178 b *= temp; /* b = (x/2)^n */
179 }
180 b = b / a;
181 }
182 }
183 else
184 {
185 /* use backward recurrence */
186 /* x x^2 x^2
187 * J(n,x)/J(n-1,x) = ---- ------ ------ .....
188 * 2n - 2(n+1) - 2(n+2)
189 *
190 * 1 1 1
191 * (for large x) = ---- ------ ------ .....
192 * 2n 2(n+1) 2(n+2)
193 * -- - ------ - ------ -
194 * x x x
195 *
196 * Let w = 2n/x and h=2/x, then the above quotient
197 * is equal to the continued fraction:
198 * 1
199 * = -----------------------
200 * 1
201 * w - -----------------
202 * 1
203 * w+h - ---------
204 * w+2h - ...
205 *
206 * To determine how many terms needed, let
207 * Q(0) = w, Q(1) = w(w+h) - 1,
208 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
209 * When Q(k) > 1e4 good for single
210 * When Q(k) > 1e9 good for double
211 * When Q(k) > 1e17 good for quadruple
212 */
213 /* determine k */
214 __float128 t, v;
215 __float128 q0, q1, h, tmp;
216 int32_t k, m;
217 w = (n + n) / (__float128) x;
218 h = 2.0Q / (__float128) x;
219 q0 = w;
220 z = w + h;
221 q1 = w * z - 1.0Q;
222 k = 1;
223 while (q1 < 1.0e17Q)
224 {
225 k += 1;
226 z += h;
227 tmp = z * q1 - q0;
228 q0 = q1;
229 q1 = tmp;
230 }
231 m = n + n;
232 for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
233 t = one / (i / x - t);
234 a = t;
235 b = one;
236 /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
237 * Hence, if n*(log(2n/x)) > ...
238 * single 8.8722839355e+01
239 * double 7.09782712893383973096e+02
240 * __float128 1.1356523406294143949491931077970765006170e+04
241 * then recurrent value may overflow and the result is
242 * likely underflow to zero
243 */
244 tmp = n;
245 v = two / x;
246 tmp = tmp * logq (fabsq (v * tmp));
247
248 if (tmp < 1.1356523406294143949491931077970765006170e+04Q)
249 {
250 for (i = n - 1, di = (__float128) (i + i); i > 0; i--)
251 {
252 temp = b;
253 b *= di;
254 b = b / x - a;
255 a = temp;
256 di -= two;
257 }
258 }
259 else
260 {
261 for (i = n - 1, di = (__float128) (i + i); i > 0; i--)
262 {
263 temp = b;
264 b *= di;
265 b = b / x - a;
266 a = temp;
267 di -= two;
268 /* scale b to avoid spurious overflow */
269 if (b > 1e100Q)
270 {
271 a /= b;
272 t /= b;
273 b = one;
274 }
275 }
276 }
277 /* j0() and j1() suffer enormous loss of precision at and
278 * near zero; however, we know that their zero points never
279 * coincide, so just choose the one further away from zero.
280 */
281 z = j0q (x);
282 w = j1q (x);
283 if (fabsq (z) >= fabsq (w))
284 b = (t * z / b);
285 else
286 b = (t * w / a);
287 }
288 }
289 if (sgn == 1)
290 return -b;
291 else
292 return b;
293 }
294
295 __float128
ynq(int n,__float128 x)296 ynq (int n, __float128 x)
297 {
298 uint32_t se;
299 int32_t i, ix;
300 int32_t sign;
301 __float128 a, b, temp;
302 ieee854_float128 u;
303
304 u.value = x;
305 se = u.words32.w0;
306 ix = se & 0x7fffffff;
307
308 /* if Y(n,NaN) is NaN */
309 if (ix >= 0x7fff0000)
310 {
311 if ((u.words32.w0 & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3)
312 return x + x;
313 }
314 if (x <= 0.0Q)
315 {
316 if (x == 0.0Q)
317 return -HUGE_VALQ + x;
318 if (se & 0x80000000)
319 return zero / (zero * x);
320 }
321 sign = 1;
322 if (n < 0)
323 {
324 n = -n;
325 sign = 1 - ((n & 1) << 1);
326 }
327 if (n == 0)
328 return (y0q (x));
329 if (n == 1)
330 return (sign * y1q (x));
331 if (ix >= 0x7fff0000)
332 return zero;
333 if (ix >= 0x412D0000)
334 { /* x > 2**302 */
335
336 /* ??? See comment above on the possible futility of this. */
337
338 /* (x >> n**2)
339 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
340 * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
341 * Let s=sin(x), c=cos(x),
342 * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
343 *
344 * n sin(xn)*sqt2 cos(xn)*sqt2
345 * ----------------------------------
346 * 0 s-c c+s
347 * 1 -s-c -c+s
348 * 2 -s+c -c-s
349 * 3 s+c c-s
350 */
351 __float128 s;
352 __float128 c;
353 sincosq (x, &s, &c);
354 switch (n & 3)
355 {
356 case 0:
357 temp = s - c;
358 break;
359 case 1:
360 temp = -s - c;
361 break;
362 case 2:
363 temp = -s + c;
364 break;
365 case 3:
366 temp = s + c;
367 break;
368 }
369 b = invsqrtpi * temp / sqrtq (x);
370 }
371 else
372 {
373 a = y0q (x);
374 b = y1q (x);
375 /* quit if b is -inf */
376 u.value = b;
377 se = u.words32.w0 & 0xffff0000;
378 for (i = 1; i < n && se != 0xffff0000; i++)
379 {
380 temp = b;
381 b = ((__float128) (i + i) / x) * b - a;
382 u.value = b;
383 se = u.words32.w0 & 0xffff0000;
384 a = temp;
385 }
386 }
387 /* If B is +-Inf, set up errno accordingly. */
388 if (! finiteq (b))
389 errno = ERANGE;
390 if (sign > 0)
391 return b;
392 else
393 return -b;
394 }
395