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 >= 0x7ff00000)
00100     {
00101       if ((u.parts32.w0 & 0xfffff) | u.parts32.w1
00102          | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3)
00103        return x + x;
00104     }
00105 
00106   if (n < 0)
00107     {
00108       n = -n;
00109       x = -x;
00110       se ^= 0x80000000;
00111     }
00112   if (n == 0)
00113     return (__ieee754_j0l (x));
00114   if (n == 1)
00115     return (__ieee754_j1l (x));
00116   sgn = (n & 1) & (se >> 31);      /* even n -- 0, odd n -- sign(x) */
00117   x = fabsl (x);
00118 
00119   if (x == 0.0L || ix >= 0x7ff00000)      /* if x is 0 or inf */
00120     b = zero;
00121   else if ((long double) n <= x)
00122     {
00123       /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
00124       if (ix >= 0x52d00000)
00125        {                    /* x > 2**302 */
00126 
00127          /* ??? Could use an expansion for large x here.  */
00128 
00129          /* (x >> n**2)
00130           *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
00131           *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
00132           *      Let s=sin(x), c=cos(x),
00133           *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
00134           *
00135           *             n    sin(xn)*sqt2    cos(xn)*sqt2
00136           *          ----------------------------------
00137           *             0     s-c             c+s
00138           *             1    -s-c            -c+s
00139           *             2    -s+c            -c-s
00140           *             3     s+c             c-s
00141           */
00142          long double s;
00143          long double c;
00144          __sincosl (x, &s, &c);
00145          switch (n & 3)
00146            {
00147            case 0:
00148              temp = c + s;
00149              break;
00150            case 1:
00151              temp = -c + s;
00152              break;
00153            case 2:
00154              temp = -c - s;
00155              break;
00156            case 3:
00157              temp = c - s;
00158              break;
00159            }
00160          b = invsqrtpi * temp / __ieee754_sqrtl (x);
00161        }
00162       else
00163        {
00164          a = __ieee754_j0l (x);
00165          b = __ieee754_j1l (x);
00166          for (i = 1; i < n; i++)
00167            {
00168              temp = b;
00169              b = b * ((long double) (i + i) / x) - a;   /* avoid underflow */
00170              a = temp;
00171            }
00172        }
00173     }
00174   else
00175     {
00176       if (ix < 0x3e100000)
00177        {                    /* x < 2**-29 */
00178          /* x is tiny, return the first Taylor expansion of J(n,x)
00179           * J(n,x) = 1/n!*(x/2)^n  - ...
00180           */
00181          if (n >= 33)              /* underflow, result < 10^-300 */
00182            b = zero;
00183          else
00184            {
00185              temp = x * 0.5;
00186              b = temp;
00187              for (a = one, i = 2; i <= n; i++)
00188               {
00189                 a *= (long double) i;     /* a = n! */
00190                 b *= temp;  /* b = (x/2)^n */
00191               }
00192              b = b / a;
00193            }
00194        }
00195       else
00196        {
00197          /* use backward recurrence */
00198          /*                      x      x^2      x^2
00199           *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
00200           *                      2n  - 2(n+1) - 2(n+2)
00201           *
00202           *                      1      1        1
00203           *  (for large x)   =  ----  ------   ------   .....
00204           *                      2n   2(n+1)   2(n+2)
00205           *                      -- - ------ - ------ -
00206           *                       x     x         x
00207           *
00208           * Let w = 2n/x and h=2/x, then the above quotient
00209           * is equal to the continued fraction:
00210           *                  1
00211           *      = -----------------------
00212           *                     1
00213           *         w - -----------------
00214           *                        1
00215           *              w+h - ---------
00216           *                     w+2h - ...
00217           *
00218           * To determine how many terms needed, let
00219           * Q(0) = w, Q(1) = w(w+h) - 1,
00220           * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
00221           * When Q(k) > 1e4      good for single
00222           * When Q(k) > 1e9      good for double
00223           * When Q(k) > 1e17     good for quadruple
00224           */
00225          /* determine k */
00226          long double t, v;
00227          long double q0, q1, h, tmp;
00228          int32_t k, m;
00229          w = (n + n) / (long double) x;
00230          h = 2.0L / (long double) x;
00231          q0 = w;
00232          z = w + h;
00233          q1 = w * z - 1.0L;
00234          k = 1;
00235          while (q1 < 1.0e17L)
00236            {
00237              k += 1;
00238              z += h;
00239              tmp = z * q1 - q0;
00240              q0 = q1;
00241              q1 = tmp;
00242            }
00243          m = n + n;
00244          for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
00245            t = one / (i / x - t);
00246          a = t;
00247          b = one;
00248          /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
00249           *  Hence, if n*(log(2n/x)) > ...
00250           *  single 8.8722839355e+01
00251           *  double 7.09782712893383973096e+02
00252           *  long double 1.1356523406294143949491931077970765006170e+04
00253           *  then recurrent value may overflow and the result is
00254           *  likely underflow to zero
00255           */
00256          tmp = n;
00257          v = two / x;
00258          tmp = tmp * __ieee754_logl (fabsl (v * tmp));
00259 
00260          if (tmp < 1.1356523406294143949491931077970765006170e+04L)
00261            {
00262              for (i = n - 1, di = (long double) (i + i); i > 0; i--)
00263               {
00264                 temp = b;
00265                 b *= di;
00266                 b = b / x - a;
00267                 a = temp;
00268                 di -= two;
00269               }
00270            }
00271          else
00272            {
00273              for (i = n - 1, di = (long double) (i + i); i > 0; i--)
00274               {
00275                 temp = b;
00276                 b *= di;
00277                 b = b / x - a;
00278                 a = temp;
00279                 di -= two;
00280                 /* scale b to avoid spurious overflow */
00281                 if (b > 1e100L)
00282                   {
00283                     a /= b;
00284                     t /= b;
00285                     b = one;
00286                   }
00287               }
00288            }
00289          b = (t * __ieee754_j0l (x) / b);
00290        }
00291     }
00292   if (sgn == 1)
00293     return -b;
00294   else
00295     return b;
00296 }
00297 
00298 #ifdef __STDC__
00299 long double
00300 __ieee754_ynl (int n, long double x)
00301 #else
00302 long double
00303 __ieee754_ynl (n, x)
00304      int n;
00305      long double x;
00306 #endif
00307 {
00308   u_int32_t se;
00309   int32_t i, ix;
00310   int32_t sign;
00311   long double a, b, temp;
00312   ieee854_long_double_shape_type u;
00313 
00314   u.value = x;
00315   se = u.parts32.w0;
00316   ix = se & 0x7fffffff;
00317 
00318   /* if Y(n,NaN) is NaN */
00319   if (ix >= 0x7ff00000)
00320     {
00321       if ((u.parts32.w0 & 0xfffff) | u.parts32.w1
00322          | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3)
00323        return x + x;
00324     }
00325   if (x <= 0.0L)
00326     {
00327       if (x == 0.0L)
00328        return -HUGE_VALL + x;
00329       if (se & 0x80000000)
00330        return zero / (zero * x);
00331     }
00332   sign = 1;
00333   if (n < 0)
00334     {
00335       n = -n;
00336       sign = 1 - ((n & 1) << 1);
00337     }
00338   if (n == 0)
00339     return (__ieee754_y0l (x));
00340   if (n == 1)
00341     return (sign * __ieee754_y1l (x));
00342   if (ix >= 0x7ff00000)
00343     return zero;
00344   if (ix >= 0x52D00000)
00345     {                       /* x > 2**302 */
00346 
00347       /* ??? See comment above on the possible futility of this.  */
00348 
00349       /* (x >> n**2)
00350        *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
00351        *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
00352        *      Let s=sin(x), c=cos(x),
00353        *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
00354        *
00355        *             n    sin(xn)*sqt2    cos(xn)*sqt2
00356        *          ----------------------------------
00357        *             0     s-c             c+s
00358        *             1    -s-c            -c+s
00359        *             2    -s+c            -c-s
00360        *             3     s+c             c-s
00361        */
00362       long double s;
00363       long double c;
00364       __sincosl (x, &s, &c);
00365       switch (n & 3)
00366        {
00367        case 0:
00368          temp = s - c;
00369          break;
00370        case 1:
00371          temp = -s - c;
00372          break;
00373        case 2:
00374          temp = -s + c;
00375          break;
00376        case 3:
00377          temp = s + c;
00378          break;
00379        }
00380       b = invsqrtpi * temp / __ieee754_sqrtl (x);
00381     }
00382   else
00383     {
00384       a = __ieee754_y0l (x);
00385       b = __ieee754_y1l (x);
00386       /* quit if b is -inf */
00387       u.value = b;
00388       se = u.parts32.w0 & 0xfff00000;
00389       for (i = 1; i < n && se != 0xfff00000; i++)
00390        {
00391          temp = b;
00392          b = ((long double) (i + i) / x) * b - a;
00393          u.value = b;
00394          se = u.parts32.w0 & 0xfff00000;
00395          a = temp;
00396        }
00397     }
00398   if (sign > 0)
00399     return b;
00400   else
00401     return -b;
00402 }