(__ieee754_sqrt): Remove unused variables b and n.
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / ldbl-96 / 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 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.64189583547756286948079e-1L, two = 2.0e0L, one = 1.0e0L;
50
51 #ifdef __STDC__
52 static const long double zero = 0.0L;
53 #else
54 static long double zero = 0.0L;
55 #endif
56
57 #ifdef __STDC__
58 long double
59 __ieee754_jnl (int n, long double x)
60 #else
61 long double
62 __ieee754_jnl (n, x)
63      int n;
64      long double x;
65 #endif
66 {
67   u_int32_t se, i0, i1;
68   int32_t i, ix, sgn;
69   long double a, b, temp, di;
70   long double z, w;
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   GET_LDOUBLE_WORDS (se, i0, i1, x);
77   ix = se & 0x7fff;
78
79   /* if J(n,NaN) is NaN */
80   if ((ix == 0x7fff) && ((i0 & 0x7fffffff) != 0))
81     return x + x;
82   if (n < 0)
83     {
84       n = -n;
85       x = -x;
86       se ^= 0x8000;
87     }
88   if (n == 0)
89     return (__ieee754_j0l (x));
90   if (n == 1)
91     return (__ieee754_j1l (x));
92   sgn = (n & 1) & (se >> 15);   /* even n -- 0, odd n -- sign(x) */
93   x = fabsl (x);
94   if ((ix | i0 | i1) == 0 || ix >= 0x7fff)      /* if x is 0 or inf */
95     b = zero;
96   else if ((long double) n <= x)
97     {
98       /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
99       if (ix >= 0x412D)
100         {                       /* x > 2**302 */
101
102           /* ??? This might be a futile gesture.
103              If x exceeds X_TLOSS anyway, the wrapper function
104              will set the result to zero. */
105
106           /* (x >> n**2)
107            *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
108            *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
109            *      Let s=sin(x), c=cos(x),
110            *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
111            *
112            *             n    sin(xn)*sqt2    cos(xn)*sqt2
113            *          ----------------------------------
114            *             0     s-c             c+s
115            *             1    -s-c            -c+s
116            *             2    -s+c            -c-s
117            *             3     s+c             c-s
118            */
119           long double s;
120           long double c;
121           __sincosl (x, &s, &c);
122           switch (n & 3)
123             {
124             case 0:
125               temp = c + s;
126               break;
127             case 1:
128               temp = -c + s;
129               break;
130             case 2:
131               temp = -c - s;
132               break;
133             case 3:
134               temp = c - s;
135               break;
136             }
137           b = invsqrtpi * temp / __ieee754_sqrtl (x);
138         }
139       else
140         {
141           a = __ieee754_j0l (x);
142           b = __ieee754_j1l (x);
143           for (i = 1; i < n; i++)
144             {
145               temp = b;
146               b = b * ((long double) (i + i) / x) - a;  /* avoid underflow */
147               a = temp;
148             }
149         }
150     }
151   else
152     {
153       if (ix < 0x3fde)
154         {                       /* x < 2**-33 */
155           /* x is tiny, return the first Taylor expansion of J(n,x)
156            * J(n,x) = 1/n!*(x/2)^n  - ...
157            */
158           if (n >= 400)         /* underflow, result < 10^-4952 */
159             b = zero;
160           else
161             {
162               temp = x * 0.5;
163               b = temp;
164               for (a = one, i = 2; i <= n; i++)
165                 {
166                   a *= (long double) i; /* a = n! */
167                   b *= temp;    /* b = (x/2)^n */
168                 }
169               b = b / a;
170             }
171         }
172       else
173         {
174           /* use backward recurrence */
175           /*                      x      x^2      x^2
176            *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
177            *                      2n  - 2(n+1) - 2(n+2)
178            *
179            *                      1      1        1
180            *  (for large x)   =  ----  ------   ------   .....
181            *                      2n   2(n+1)   2(n+2)
182            *                      -- - ------ - ------ -
183            *                       x     x         x
184            *
185            * Let w = 2n/x and h=2/x, then the above quotient
186            * is equal to the continued fraction:
187            *                  1
188            *      = -----------------------
189            *                     1
190            *         w - -----------------
191            *                        1
192            *              w+h - ---------
193            *                     w+2h - ...
194            *
195            * To determine how many terms needed, let
196            * Q(0) = w, Q(1) = w(w+h) - 1,
197            * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
198            * When Q(k) > 1e4      good for single
199            * When Q(k) > 1e9      good for double
200            * When Q(k) > 1e17     good for quadruple
201            */
202           /* determine k */
203           long double t, v;
204           long double q0, q1, h, tmp;
205           int32_t k, m;
206           w = (n + n) / (long double) x;
207           h = 2.0L / (long double) x;
208           q0 = w;
209           z = w + h;
210           q1 = w * z - 1.0L;
211           k = 1;
212           while (q1 < 1.0e11L)
213             {
214               k += 1;
215               z += h;
216               tmp = z * q1 - q0;
217               q0 = q1;
218               q1 = tmp;
219             }
220           m = n + n;
221           for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
222             t = one / (i / x - t);
223           a = t;
224           b = one;
225           /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
226            *  Hence, if n*(log(2n/x)) > ...
227            *  single 8.8722839355e+01
228            *  double 7.09782712893383973096e+02
229            *  long double 1.1356523406294143949491931077970765006170e+04
230            *  then recurrent value may overflow and the result is
231            *  likely underflow to zero
232            */
233           tmp = n;
234           v = two / x;
235           tmp = tmp * __ieee754_logl (fabsl (v * tmp));
236
237           if (tmp < 1.1356523406294143949491931077970765006170e+04L)
238             {
239               for (i = n - 1, di = (long double) (i + i); i > 0; i--)
240                 {
241                   temp = b;
242                   b *= di;
243                   b = b / x - a;
244                   a = temp;
245                   di -= two;
246                 }
247             }
248           else
249             {
250               for (i = n - 1, di = (long double) (i + i); i > 0; i--)
251                 {
252                   temp = b;
253                   b *= di;
254                   b = b / x - a;
255                   a = temp;
256                   di -= two;
257                   /* scale b to avoid spurious overflow */
258                   if (b > 1e100L)
259                     {
260                       a /= b;
261                       t /= b;
262                       b = one;
263                     }
264                 }
265             }
266           b = (t * __ieee754_j0l (x) / b);
267         }
268     }
269   if (sgn == 1)
270     return -b;
271   else
272     return b;
273 }
274
275 #ifdef __STDC__
276 long double
277 __ieee754_ynl (int n, long double x)
278 #else
279 long double
280 __ieee754_ynl (n, x)
281      int n;
282      long double x;
283 #endif
284 {
285   u_int32_t se, i0, i1;
286   int32_t i, ix;
287   int32_t sign;
288   long double a, b, temp;
289
290
291   GET_LDOUBLE_WORDS (se, i0, i1, x);
292   ix = se & 0x7fff;
293   /* if Y(n,NaN) is NaN */
294   if ((ix == 0x7fff) && ((i0 & 0x7fffffff) != 0))
295     return x + x;
296   if ((ix | i0 | i1) == 0)
297     return -one / zero;
298   if (se & 0x8000)
299     return zero / zero;
300   sign = 1;
301   if (n < 0)
302     {
303       n = -n;
304       sign = 1 - ((n & 1) << 1);
305     }
306   if (n == 0)
307     return (__ieee754_y0l (x));
308   if (n == 1)
309     return (sign * __ieee754_y1l (x));
310   if (ix == 0x7fff)
311     return zero;
312   if (ix >= 0x412D)
313     {                           /* x > 2**302 */
314
315       /* ??? See comment above on the possible futility of this.  */
316
317       /* (x >> n**2)
318        *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
319        *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
320        *      Let s=sin(x), c=cos(x),
321        *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
322        *
323        *             n    sin(xn)*sqt2    cos(xn)*sqt2
324        *          ----------------------------------
325        *             0     s-c             c+s
326        *             1    -s-c            -c+s
327        *             2    -s+c            -c-s
328        *             3     s+c             c-s
329        */
330       long double s;
331       long double c;
332       __sincosl (x, &s, &c);
333       switch (n & 3)
334         {
335         case 0:
336           temp = s - c;
337           break;
338         case 1:
339           temp = -s - c;
340           break;
341         case 2:
342           temp = -s + c;
343           break;
344         case 3:
345           temp = s + c;
346           break;
347         }
348       b = invsqrtpi * temp / __ieee754_sqrtl (x);
349     }
350   else
351     {
352       a = __ieee754_y0l (x);
353       b = __ieee754_y1l (x);
354       /* quit if b is -inf */
355       GET_LDOUBLE_WORDS (se, i0, i1, b);
356       for (i = 1; i < n && se != 0xffff; i++)
357         {
358           temp = b;
359           b = ((long double) (i + i) / x) * b - a;
360           GET_LDOUBLE_WORDS (se, i0, i1, b);
361           a = temp;
362         }
363     }
364   if (sign > 0)
365     return b;
366   else
367     return -b;
368 }