Back to index

glibc  2.9
e_jnl.c
Go to the documentation of this file.
00001 /*
00002  * ====================================================
00003  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00004  *
00005  * Developed at SunPro, a Sun Microsystems, Inc. business.
00006  * Permission to use, copy, modify, and distribute this
00007  * software is freely granted, provided that this notice
00008  * is preserved.
00009  * ====================================================
00010  */
00011 
00012 /* Modifications for 128-bit long double are
00013    Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
00014    and are incorporated herein by permission of the author.  The author 
00015    reserves the right to distribute this material elsewhere under different
00016    copying permissions.  These modifications are distributed here under 
00017    the following terms:
00018 
00019     This library is free software; you can redistribute it and/or
00020     modify it under the terms of the GNU Lesser General Public
00021     License as published by the Free Software Foundation; either
00022     version 2.1 of the License, or (at your option) any later version.
00023 
00024     This library is distributed in the hope that it will be useful,
00025     but WITHOUT ANY WARRANTY; without even the implied warranty of
00026     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00027     Lesser General Public License for more details.
00028 
00029     You should have received a copy of the GNU Lesser General Public
00030     License along with this library; if not, write to the Free Software
00031     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00032 
00033 /*
00034  * __ieee754_jn(n, x), __ieee754_yn(n, x)
00035  * floating point Bessel's function of the 1st and 2nd kind
00036  * of order n
00037  *
00038  * Special cases:
00039  *     y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
00040  *     y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
00041  * Note 2. About jn(n,x), yn(n,x)
00042  *     For n=0, j0(x) is called,
00043  *     for n=1, j1(x) is called,
00044  *     for n<x, forward recursion us used starting
00045  *     from values of j0(x) and j1(x).
00046  *     for n>x, a continued fraction approximation to
00047  *     j(n,x)/j(n-1,x) is evaluated and then backward
00048  *     recursion is used starting from a supposed value
00049  *     for j(n,x). The resulting value of j(0,x) is
00050  *     compared with the actual value to correct the
00051  *     supposed value of j(n,x).
00052  *
00053  *     yn(n,x) is similar in all respects, except
00054  *     that forward recursion is used for all
00055  *     values of n>1.
00056  *
00057  */
00058 
00059 #include "math.h"
00060 #include "math_private.h"
00061 
00062 #ifdef __STDC__
00063 static const long double
00064 #else
00065 static long double
00066 #endif
00067   invsqrtpi = 5.6418958354775628694807945156077258584405E-1L,
00068   two = 2.0e0L,
00069   one = 1.0e0L,
00070   zero = 0.0L;
00071 
00072 
00073 #ifdef __STDC__
00074 long double
00075 __ieee754_jnl (int n, long double x)
00076 #else
00077 long double
00078 __ieee754_jnl (n, x)
00079      int n;
00080      long double x;
00081 #endif
00082 {
00083   u_int32_t se;
00084   int32_t i, ix, sgn;
00085   long double a, b, temp, di;
00086   long double z, w;
00087   ieee854_long_double_shape_type u;
00088 
00089 
00090   /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
00091    * Thus, J(-n,x) = J(n,-x)
00092    */
00093 
00094   u.value = x;
00095   se = u.parts32.w0;
00096   ix = se & 0x7fffffff;
00097 
00098   /* if J(n,NaN) is NaN */
00099   if (ix >= 0x7fff0000)
00100     {
00101       if ((u.parts32.w0 & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3)
00102        return x + x;
00103     }
00104 
00105   if (n < 0)
00106     {
00107       n = -n;
00108       x = -x;
00109       se ^= 0x80000000;
00110     }
00111   if (n == 0)
00112     return (__ieee754_j0l (x));
00113   if (n == 1)
00114     return (__ieee754_j1l (x));
00115   sgn = (n & 1) & (se >> 31);      /* even n -- 0, odd n -- sign(x) */
00116   x = fabsl (x);
00117 
00118   if (x == 0.0L || ix >= 0x7fff0000)      /* if x is 0 or inf */
00119     b = zero;
00120   else if ((long double) n <= x)
00121     {
00122       /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
00123       if (ix >= 0x412D0000)
00124        {                    /* x > 2**302 */
00125 
00126          /* ??? Could use an expansion for large x here.  */
00127 
00128          /* (x >> n**2)
00129           *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
00130           *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
00131           *      Let s=sin(x), c=cos(x),
00132           *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
00133           *
00134           *             n    sin(xn)*sqt2    cos(xn)*sqt2
00135           *          ----------------------------------
00136           *             0     s-c             c+s
00137           *             1    -s-c            -c+s
00138           *             2    -s+c            -c-s
00139           *             3     s+c             c-s
00140           */
00141          long double s;
00142          long double c;
00143          __sincosl (x, &s, &c);
00144          switch (n & 3)
00145            {
00146            case 0:
00147              temp = c + s;
00148              break;
00149            case 1:
00150              temp = -c + s;
00151              break;
00152            case 2:
00153              temp = -c - s;
00154              break;
00155            case 3:
00156              temp = c - s;
00157              break;
00158            }
00159          b = invsqrtpi * temp / __ieee754_sqrtl (x);
00160        }
00161       else
00162        {
00163          a = __ieee754_j0l (x);
00164          b = __ieee754_j1l (x);
00165          for (i = 1; i < n; i++)
00166            {
00167              temp = b;
00168              b = b * ((long double) (i + i) / x) - a;   /* avoid underflow */
00169              a = temp;
00170            }
00171        }
00172     }
00173   else
00174     {
00175       if (ix < 0x3fc60000)
00176        {                    /* x < 2**-57 */
00177          /* x is tiny, return the first Taylor expansion of J(n,x)
00178           * J(n,x) = 1/n!*(x/2)^n  - ...
00179           */
00180          if (n >= 400)             /* underflow, result < 10^-4952 */
00181            b = zero;
00182          else
00183            {
00184              temp = x * 0.5;
00185              b = temp;
00186              for (a = one, i = 2; i <= n; i++)
00187               {
00188                 a *= (long double) i;     /* a = n! */
00189                 b *= temp;  /* b = (x/2)^n */
00190               }
00191              b = b / a;
00192            }
00193        }
00194       else
00195        {
00196          /* use backward recurrence */
00197          /*                      x      x^2      x^2
00198           *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
00199           *                      2n  - 2(n+1) - 2(n+2)
00200           *
00201           *                      1      1        1
00202           *  (for large x)   =  ----  ------   ------   .....
00203           *                      2n   2(n+1)   2(n+2)
00204           *                      -- - ------ - ------ -
00205           *                       x     x         x
00206           *
00207           * Let w = 2n/x and h=2/x, then the above quotient
00208           * is equal to the continued fraction:
00209           *                  1
00210           *      = -----------------------
00211           *                     1
00212           *         w - -----------------
00213           *                        1
00214           *              w+h - ---------
00215           *                     w+2h - ...
00216           *
00217           * To determine how many terms needed, let
00218           * Q(0) = w, Q(1) = w(w+h) - 1,
00219           * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
00220           * When Q(k) > 1e4      good for single
00221           * When Q(k) > 1e9      good for double
00222           * When Q(k) > 1e17     good for quadruple
00223           */
00224          /* determine k */
00225          long double t, v;
00226          long double q0, q1, h, tmp;
00227          int32_t k, m;
00228          w = (n + n) / (long double) x;
00229          h = 2.0L / (long double) x;
00230          q0 = w;
00231          z = w + h;
00232          q1 = w * z - 1.0L;
00233          k = 1;
00234          while (q1 < 1.0e17L)
00235            {
00236              k += 1;
00237              z += h;
00238              tmp = z * q1 - q0;
00239              q0 = q1;
00240              q1 = tmp;
00241            }
00242          m = n + n;
00243          for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
00244            t = one / (i / x - t);
00245          a = t;
00246          b = one;
00247          /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
00248           *  Hence, if n*(log(2n/x)) > ...
00249           *  single 8.8722839355e+01
00250           *  double 7.09782712893383973096e+02
00251           *  long double 1.1356523406294143949491931077970765006170e+04
00252           *  then recurrent value may overflow and the result is
00253           *  likely underflow to zero
00254           */
00255          tmp = n;
00256          v = two / x;
00257          tmp = tmp * __ieee754_logl (fabsl (v * tmp));
00258 
00259          if (tmp < 1.1356523406294143949491931077970765006170e+04L)
00260            {
00261              for (i = n - 1, di = (long double) (i + i); i > 0; i--)
00262               {
00263                 temp = b;
00264                 b *= di;
00265                 b = b / x - a;
00266                 a = temp;
00267                 di -= two;
00268               }
00269            }
00270          else
00271            {
00272              for (i = n - 1, di = (long double) (i + i); i > 0; i--)
00273               {
00274                 temp = b;
00275                 b *= di;
00276                 b = b / x - a;
00277                 a = temp;
00278                 di -= two;
00279                 /* scale b to avoid spurious overflow */
00280                 if (b > 1e100L)
00281                   {
00282                     a /= b;
00283                     t /= b;
00284                     b = one;
00285                   }
00286               }
00287            }
00288          b = (t * __ieee754_j0l (x) / b);
00289        }
00290     }
00291   if (sgn == 1)
00292     return -b;
00293   else
00294     return b;
00295 }
00296 
00297 #ifdef __STDC__
00298 long double
00299 __ieee754_ynl (int n, long double x)
00300 #else
00301 long double
00302 __ieee754_ynl (n, x)
00303      int n;
00304      long double x;
00305 #endif
00306 {
00307   u_int32_t se;
00308   int32_t i, ix;
00309   int32_t sign;
00310   long double a, b, temp;
00311   ieee854_long_double_shape_type u;
00312 
00313   u.value = x;
00314   se = u.parts32.w0;
00315   ix = se & 0x7fffffff;
00316 
00317   /* if Y(n,NaN) is NaN */
00318   if (ix >= 0x7fff0000)
00319     {
00320       if ((u.parts32.w0 & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3)
00321        return x + x;
00322     }
00323   if (x <= 0.0L)
00324     {
00325       if (x == 0.0L)
00326        return -HUGE_VALL + x;
00327       if (se & 0x80000000)
00328        return zero / (zero * x);
00329     }
00330   sign = 1;
00331   if (n < 0)
00332     {
00333       n = -n;
00334       sign = 1 - ((n & 1) << 1);
00335     }
00336   if (n == 0)
00337     return (__ieee754_y0l (x));
00338   if (n == 1)
00339     return (sign * __ieee754_y1l (x));
00340   if (ix >= 0x7fff0000)
00341     return zero;
00342   if (ix >= 0x412D0000)
00343     {                       /* x > 2**302 */
00344 
00345       /* ??? See comment above on the possible futility of this.  */
00346 
00347       /* (x >> n**2)
00348        *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
00349        *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
00350        *      Let s=sin(x), c=cos(x),
00351        *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
00352        *
00353        *             n    sin(xn)*sqt2    cos(xn)*sqt2
00354        *          ----------------------------------
00355        *             0     s-c             c+s
00356        *             1    -s-c            -c+s
00357        *             2    -s+c            -c-s
00358        *             3     s+c             c-s
00359        */
00360       long double s;
00361       long double c;
00362       __sincosl (x, &s, &c);
00363       switch (n & 3)
00364        {
00365        case 0:
00366          temp = s - c;
00367          break;
00368        case 1:
00369          temp = -s - c;
00370          break;
00371        case 2:
00372          temp = -s + c;
00373          break;
00374        case 3:
00375          temp = s + c;
00376          break;
00377        }
00378       b = invsqrtpi * temp / __ieee754_sqrtl (x);
00379     }
00380   else
00381     {
00382       a = __ieee754_y0l (x);
00383       b = __ieee754_y1l (x);
00384       /* quit if b is -inf */
00385       u.value = b;
00386       se = u.parts32.w0 & 0xffff0000;
00387       for (i = 1; i < n && se != 0xffff0000; i++)
00388        {
00389          temp = b;
00390          b = ((long double) (i + i) / x) * b - a;
00391          u.value = b;
00392          se = u.parts32.w0 & 0xffff0000;
00393          a = temp;
00394        }
00395     }
00396   if (sign > 0)
00397     return b;
00398   else
00399     return -b;
00400 }