128-bit long double Bessel functions jn and yn.
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / ldbl-128 / e_jnl.c
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 contributed by
13    Stephen L. Moshier <moshier@na-net.ornl.gov>  */
14
15 /*
16  * __ieee754_jn(n, x), __ieee754_yn(n, x)
17  * floating point Bessel's function of the 1st and 2nd kind
18  * of order n
19  *
20  * Special cases:
21  *      y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
22  *      y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
23  * Note 2. About jn(n,x), yn(n,x)
24  *      For n=0, j0(x) is called,
25  *      for n=1, j1(x) is called,
26  *      for n<x, forward recursion us used starting
27  *      from values of j0(x) and j1(x).
28  *      for n>x, a continued fraction approximation to
29  *      j(n,x)/j(n-1,x) is evaluated and then backward
30  *      recursion is used starting from a supposed value
31  *      for j(n,x). The resulting value of j(0,x) is
32  *      compared with the actual value to correct the
33  *      supposed value of j(n,x).
34  *
35  *      yn(n,x) is similar in all respects, except
36  *      that forward recursion is used for all
37  *      values of n>1.
38  *
39  */
40
41 #include "math.h"
42 #include "math_private.h"
43
44 #ifdef __STDC__
45 static const long double
46 #else
47 static long double
48 #endif
49   invsqrtpi = 5.6418958354775628694807945156077258584405E-1L,
50   two = 2.0e0L,
51   one = 1.0e0L,
52   zero = 0.0L;
53
54
55 #ifdef __STDC__
56 long double
57 __ieee754_jnl (int n, long double x)
58 #else
59 long double
60 __ieee754_jnl (n, x)
61      int n;
62      long double x;
63 #endif
64 {
65   u_int32_t se;
66   int32_t i, ix, sgn;
67   long double a, b, temp, di;
68   long double z, w;
69   ieee854_long_double_shape_type u;
70
71
72   /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
73    * Thus, J(-n,x) = J(n,-x)
74    */
75
76   u.value = x;
77   se = u.parts32.w0;
78   ix = se & 0x7fffffff;
79
80   /* if J(n,NaN) is NaN */
81   if (ix >= 0x7fff0000)
82     {
83       if ((u.parts32.w0 & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3)
84         return x + x;
85     }
86
87   if (n < 0)
88     {
89       n = -n;
90       x = -x;
91       se ^= 0x80000000;
92     }
93   if (n == 0)
94     return (__ieee754_j0l (x));
95   if (n == 1)
96     return (__ieee754_j1l (x));
97   sgn = (n & 1) & (se >> 31);   /* even n -- 0, odd n -- sign(x) */
98   x = fabsl (x);
99
100   if (x == 0.0L || ix >= 0x7fff0000)    /* if x is 0 or inf */
101     b = zero;
102   else if ((long double) n <= x)
103     {
104       /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
105       if (ix >= 0x412D0000)
106         {                       /* x > 2**302 */
107
108           /* ??? Could use an expansion for large x here.  */
109
110           /* (x >> n**2)
111            *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
112            *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
113            *      Let s=sin(x), c=cos(x),
114            *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
115            *
116            *             n    sin(xn)*sqt2    cos(xn)*sqt2
117            *          ----------------------------------
118            *             0     s-c             c+s
119            *             1    -s-c            -c+s
120            *             2    -s+c            -c-s
121            *             3     s+c             c-s
122            */
123           long double s;
124           long double c;
125           __sincosl (x, &s, &c);
126           switch (n & 3)
127             {
128             case 0:
129               temp = c + s;
130               break;
131             case 1:
132               temp = -c + s;
133               break;
134             case 2:
135               temp = -c - s;
136               break;
137             case 3:
138               temp = c - s;
139               break;
140             }
141           b = invsqrtpi * temp / __ieee754_sqrtl (x);
142         }
143       else
144         {
145           a = __ieee754_j0l (x);
146           b = __ieee754_j1l (x);
147           for (i = 1; i < n; i++)
148             {
149               temp = b;
150               b = b * ((long double) (i + i) / x) - a;  /* avoid underflow */
151               a = temp;
152             }
153         }
154     }
155   else
156     {
157       if (ix < 0x3fc60000)
158         {                       /* x < 2**-57 */
159           /* x is tiny, return the first Taylor expansion of J(n,x)
160            * J(n,x) = 1/n!*(x/2)^n  - ...
161            */
162           if (n >= 400)         /* underflow, result < 10^-4952 */
163             b = zero;
164           else
165             {
166               temp = x * 0.5;
167               b = temp;
168               for (a = one, i = 2; i <= n; i++)
169                 {
170                   a *= (long double) i; /* a = n! */
171                   b *= temp;    /* b = (x/2)^n */
172                 }
173               b = b / a;
174             }
175         }
176       else
177         {
178           /* use backward recurrence */
179           /*                      x      x^2      x^2
180            *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
181            *                      2n  - 2(n+1) - 2(n+2)
182            *
183            *                      1      1        1
184            *  (for large x)   =  ----  ------   ------   .....
185            *                      2n   2(n+1)   2(n+2)
186            *                      -- - ------ - ------ -
187            *                       x     x         x
188            *
189            * Let w = 2n/x and h=2/x, then the above quotient
190            * is equal to the continued fraction:
191            *                  1
192            *      = -----------------------
193            *                     1
194            *         w - -----------------
195            *                        1
196            *              w+h - ---------
197            *                     w+2h - ...
198            *
199            * To determine how many terms needed, let
200            * Q(0) = w, Q(1) = w(w+h) - 1,
201            * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
202            * When Q(k) > 1e4      good for single
203            * When Q(k) > 1e9      good for double
204            * When Q(k) > 1e17     good for quadruple
205            */
206           /* determine k */
207           long double t, v;
208           long double q0, q1, h, tmp;
209           int32_t k, m;
210           w = (n + n) / (long double) x;
211           h = 2.0L / (long double) x;
212           q0 = w;
213           z = w + h;
214           q1 = w * z - 1.0L;
215           k = 1;
216           while (q1 < 1.0e17L)
217             {
218               k += 1;
219               z += h;
220               tmp = z * q1 - q0;
221               q0 = q1;
222               q1 = tmp;
223             }
224           m = n + n;
225           for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
226             t = one / (i / x - t);
227           a = t;
228           b = one;
229           /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
230            *  Hence, if n*(log(2n/x)) > ...
231            *  single 8.8722839355e+01
232            *  double 7.09782712893383973096e+02
233            *  long double 1.1356523406294143949491931077970765006170e+04
234            *  then recurrent value may overflow and the result is
235            *  likely underflow to zero
236            */
237           tmp = n;
238           v = two / x;
239           tmp = tmp * __ieee754_logl (fabsl (v * tmp));
240
241           if (tmp < 1.1356523406294143949491931077970765006170e+04L)
242             {
243               for (i = n - 1, di = (long double) (i + i); i > 0; i--)
244                 {
245                   temp = b;
246                   b *= di;
247                   b = b / x - a;
248                   a = temp;
249                   di -= two;
250                 }
251             }
252           else
253             {
254               for (i = n - 1, di = (long double) (i + i); i > 0; i--)
255                 {
256                   temp = b;
257                   b *= di;
258                   b = b / x - a;
259                   a = temp;
260                   di -= two;
261                   /* scale b to avoid spurious overflow */
262                   if (b > 1e100L)
263                     {
264                       a /= b;
265                       t /= b;
266                       b = one;
267                     }
268                 }
269             }
270           b = (t * __ieee754_j0l (x) / b);
271         }
272     }
273   if (sgn == 1)
274     return -b;
275   else
276     return b;
277 }
278
279 #ifdef __STDC__
280 long double
281 __ieee754_ynl (int n, long double x)
282 #else
283 long double
284 __ieee754_ynl (n, x)
285      int n;
286      long double x;
287 #endif
288 {
289   u_int32_t se;
290   int32_t i, ix;
291   int32_t sign;
292   long double a, b, temp;
293   ieee854_long_double_shape_type u;
294
295   u.value = x;
296   se = u.parts32.w0;
297   ix = se & 0x7fffffff;
298
299   /* if Y(n,NaN) is NaN */
300   if (ix >= 0x7fff0000)
301     {
302       if ((u.parts32.w0 & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3)
303         return x + x;
304     }
305   if (x <= 0.0L)
306     {
307       if (x == 0.0L)
308         return -one / zero;
309       if (se & 0x80000000)
310         return zero / zero;
311     }
312   sign = 1;
313   if (n < 0)
314     {
315       n = -n;
316       sign = 1 - ((n & 1) << 1);
317     }
318   if (n == 0)
319     return (__ieee754_y0l (x));
320   if (n == 1)
321     return (sign * __ieee754_y1l (x));
322   if (ix >= 0x7fff0000)
323     return zero;
324   if (ix >= 0x412D0000)
325     {                           /* x > 2**302 */
326
327       /* ??? See comment above on the possible futility of this.  */
328
329       /* (x >> n**2)
330        *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
331        *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
332        *      Let s=sin(x), c=cos(x),
333        *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
334        *
335        *             n    sin(xn)*sqt2    cos(xn)*sqt2
336        *          ----------------------------------
337        *             0     s-c             c+s
338        *             1    -s-c            -c+s
339        *             2    -s+c            -c-s
340        *             3     s+c             c-s
341        */
342       long double s;
343       long double c;
344       __sincosl (x, &s, &c);
345       switch (n & 3)
346         {
347         case 0:
348           temp = s - c;
349           break;
350         case 1:
351           temp = -s - c;
352           break;
353         case 2:
354           temp = -s + c;
355           break;
356         case 3:
357           temp = s + c;
358           break;
359         }
360       b = invsqrtpi * temp / __ieee754_sqrtl (x);
361     }
362   else
363     {
364       a = __ieee754_y0l (x);
365       b = __ieee754_y1l (x);
366       /* quit if b is -inf */
367       u.value = b;
368       se = u.parts32.w0 & 0xffff0000;
369       for (i = 1; i < n && se != 0xffff0000; i++)
370         {
371           temp = b;
372           b = ((long double) (i + i) / x) * b - a;
373           u.value = b;
374           se = u.parts32.w0 & 0xffff0000;
375           a = temp;
376         }
377     }
378   if (sign > 0)
379     return b;
380   else
381     return -b;
382 }