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