Back to index

lightning-sunbird  0.9+nobinonly
e_jn.c
Go to the documentation of this file.
00001 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
00002  *
00003  * ***** BEGIN LICENSE BLOCK *****
00004  * Version: MPL 1.1/GPL 2.0/LGPL 2.1
00005  *
00006  * The contents of this file are subject to the Mozilla Public License Version
00007  * 1.1 (the "License"); you may not use this file except in compliance with
00008  * the License. You may obtain a copy of the License at
00009  * http://www.mozilla.org/MPL/
00010  *
00011  * Software distributed under the License is distributed on an "AS IS" basis,
00012  * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
00013  * for the specific language governing rights and limitations under the
00014  * License.
00015  *
00016  * The Original Code is Mozilla Communicator client code, released
00017  * March 31, 1998.
00018  *
00019  * The Initial Developer of the Original Code is
00020  * Sun Microsystems, Inc.
00021  * Portions created by the Initial Developer are Copyright (C) 1998
00022  * the Initial Developer. All Rights Reserved.
00023  *
00024  * Contributor(s):
00025  *
00026  * Alternatively, the contents of this file may be used under the terms of
00027  * either of the GNU General Public License Version 2 or later (the "GPL"),
00028  * or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
00029  * in which case the provisions of the GPL or the LGPL are applicable instead
00030  * of those above. If you wish to allow use of your version of this file only
00031  * under the terms of either the GPL or the LGPL, and not to allow others to
00032  * use your version of this file under the terms of the MPL, indicate your
00033  * decision by deleting the provisions above and replace them with the notice
00034  * and other provisions required by the GPL or the LGPL. If you do not delete
00035  * the provisions above, a recipient may use your version of this file under
00036  * the terms of any one of the MPL, the GPL or the LGPL.
00037  *
00038  * ***** END LICENSE BLOCK ***** */
00039 
00040 /* @(#)e_jn.c 1.4 95/01/18 */
00041 /*
00042  * ====================================================
00043  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00044  *
00045  * Developed at SunSoft, a Sun Microsystems, Inc. business.
00046  * Permission to use, copy, modify, and distribute this
00047  * software is freely granted, provided that this notice 
00048  * is preserved.
00049  * ====================================================
00050  */
00051 
00052 /*
00053  * __ieee754_jn(n, x), __ieee754_yn(n, x)
00054  * floating point Bessel's function of the 1st and 2nd kind
00055  * of order n
00056  *          
00057  * Special cases:
00058  *     y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
00059  *     y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
00060  * Note 2. About jn(n,x), yn(n,x)
00061  *     For n=0, j0(x) is called,
00062  *     for n=1, j1(x) is called,
00063  *     for n<x, forward recursion us used starting
00064  *     from values of j0(x) and j1(x).
00065  *     for n>x, a continued fraction approximation to
00066  *     j(n,x)/j(n-1,x) is evaluated and then backward
00067  *     recursion is used starting from a supposed value
00068  *     for j(n,x). The resulting value of j(0,x) is
00069  *     compared with the actual value to correct the
00070  *     supposed value of j(n,x).
00071  *
00072  *     yn(n,x) is similar in all respects, except
00073  *     that forward recursion is used for all
00074  *     values of n>1.
00075  *     
00076  */
00077 
00078 #include "fdlibm.h"
00079 
00080 #ifdef __STDC__
00081 static const double
00082 #else
00083 static double
00084 #endif
00085 invsqrtpi=  5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
00086 two   =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
00087 one   =  1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
00088 
00089 static double zero  =  0.00000000000000000000e+00;
00090 
00091 #ifdef __STDC__
00092        double __ieee754_jn(int n, double x)
00093 #else
00094        double __ieee754_jn(n,x)
00095        int n; double x;
00096 #endif
00097 {
00098         fd_twoints u;
00099        int i,hx,ix,lx, sgn;
00100        double a, b, temp, di;
00101        double z, w;
00102 
00103     /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
00104      * Thus, J(-n,x) = J(n,-x)
00105      */
00106         u.d = x;
00107        hx = __HI(u);
00108        ix = 0x7fffffff&hx;
00109        lx = __LO(u);
00110     /* if J(n,NaN) is NaN */
00111        if((ix|((unsigned)(lx|-lx))>>31)>0x7ff00000) return x+x;
00112        if(n<0){             
00113               n = -n;
00114               x = -x;
00115               hx ^= 0x80000000;
00116        }
00117        if(n==0) return(__ieee754_j0(x));
00118        if(n==1) return(__ieee754_j1(x));
00119        sgn = (n&1)&(hx>>31);       /* even n -- 0, odd n -- sign(x) */
00120        x = fd_fabs(x);
00121        if((ix|lx)==0||ix>=0x7ff00000)     /* if x is 0 or inf */
00122            b = zero;
00123        else if((double)n<=x) {   
00124               /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
00125            if(ix>=0x52D00000) { /* x > 2**302 */
00126     /* (x >> n**2) 
00127      *     Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
00128      *     Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
00129      *     Let s=sin(x), c=cos(x), 
00130      *        xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
00131      *
00132      *           n   sin(xn)*sqt2  cos(xn)*sqt2
00133      *        ----------------------------------
00134      *           0    s-c           c+s
00135      *           1   -s-c          -c+s
00136      *           2   -s+c          -c-s
00137      *           3    s+c           c-s
00138      */
00139               switch(n&3) {
00140                   case 0: temp =  fd_cos(x)+fd_sin(x); break;
00141                   case 1: temp = -fd_cos(x)+fd_sin(x); break;
00142                   case 2: temp = -fd_cos(x)-fd_sin(x); break;
00143                   case 3: temp =  fd_cos(x)-fd_sin(x); break;
00144               }
00145               b = invsqrtpi*temp/fd_sqrt(x);
00146            } else {  
00147                a = __ieee754_j0(x);
00148                b = __ieee754_j1(x);
00149                for(i=1;i<n;i++){
00150                   temp = b;
00151                   b = b*((double)(i+i)/x) - a; /* avoid underflow */
00152                   a = temp;
00153                }
00154            }
00155        } else {
00156            if(ix<0x3e100000) {     /* x < 2**-29 */
00157     /* x is tiny, return the first Taylor expansion of J(n,x) 
00158      * J(n,x) = 1/n!*(x/2)^n  - ...
00159      */
00160               if(n>33)      /* underflow */
00161                   b = zero;
00162               else {
00163                   temp = x*0.5; b = temp;
00164                   for (a=one,i=2;i<=n;i++) {
00165                      a *= (double)i;             /* a = n! */
00166                      b *= temp;           /* b = (x/2)^n */
00167                   }
00168                   b = b/a;
00169               }
00170            } else {
00171               /* use backward recurrence */
00172               /*                   x      x^2      x^2       
00173                *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
00174                *                   2n  - 2(n+1) - 2(n+2)
00175                *
00176                *                   1      1        1       
00177                *  (for large x)   =  ----  ------   ------   .....
00178                *                   2n   2(n+1)   2(n+2)
00179                *                   -- - ------ - ------ - 
00180                *                    x     x         x
00181                *
00182                * Let w = 2n/x and h=2/x, then the above quotient
00183                * is equal to the continued fraction:
00184                *                1
00185                *     = -----------------------
00186                *                   1
00187                *        w - -----------------
00188                *                     1
00189                *             w+h - ---------
00190                *                   w+2h - ...
00191                *
00192                * To determine how many terms needed, let
00193                * Q(0) = w, Q(1) = w(w+h) - 1,
00194                * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
00195                * When Q(k) > 1e4   good for single 
00196                * When Q(k) > 1e9   good for double 
00197                * When Q(k) > 1e17  good for quadruple 
00198                */
00199            /* determine k */
00200               double t,v;
00201               double q0,q1,h,tmp; int k,m;
00202               w  = (n+n)/(double)x; h = 2.0/(double)x;
00203               q0 = w;  z = w+h; q1 = w*z - 1.0; k=1;
00204               while(q1<1.0e9) {
00205                      k += 1; z += h;
00206                      tmp = z*q1 - q0;
00207                      q0 = q1;
00208                      q1 = tmp;
00209               }
00210               m = n+n;
00211               for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
00212               a = t;
00213               b = one;
00214               /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
00215                *  Hence, if n*(log(2n/x)) > ...
00216                *  single 8.8722839355e+01
00217                *  double 7.09782712893383973096e+02
00218                *  long double 1.1356523406294143949491931077970765006170e+04
00219                *  then recurrent value may overflow and the result is 
00220                *  likely underflow to zero
00221                */
00222               tmp = n;
00223               v = two/x;
00224               tmp = tmp*__ieee754_log(fd_fabs(v*tmp));
00225               if(tmp<7.09782712893383973096e+02) {
00226                   for(i=n-1,di=(double)(i+i);i>0;i--){
00227                       temp = b;
00228                      b *= di;
00229                      b  = b/x - a;
00230                       a = temp;
00231                      di -= two;
00232                   }
00233               } else {
00234                   for(i=n-1,di=(double)(i+i);i>0;i--){
00235                       temp = b;
00236                      b *= di;
00237                      b  = b/x - a;
00238                       a = temp;
00239                      di -= two;
00240                   /* scale b to avoid spurious overflow */
00241                      if(b>1e100) {
00242                          a /= b;
00243                          t /= b;
00244                          b  = one;
00245                      }
00246                   }
00247               }
00248               b = (t*__ieee754_j0(x)/b);
00249            }
00250        }
00251        if(sgn==1) return -b; else return b;
00252 }
00253 
00254 #ifdef __STDC__
00255        double __ieee754_yn(int n, double x) 
00256 #else
00257        double __ieee754_yn(n,x) 
00258        int n; double x;
00259 #endif
00260 {
00261         fd_twoints u;
00262        int i,hx,ix,lx;
00263        int sign;
00264        double a, b, temp;
00265 
00266         u.d = x;
00267        hx = __HI(u);
00268        ix = 0x7fffffff&hx;
00269        lx = __LO(u);
00270     /* if Y(n,NaN) is NaN */
00271        if((ix|((unsigned)(lx|-lx))>>31)>0x7ff00000) return x+x;
00272        if((ix|lx)==0) return -one/zero;
00273        if(hx<0) return zero/zero;
00274        sign = 1;
00275        if(n<0){
00276               n = -n;
00277               sign = 1 - ((n&1)<<1);
00278        }
00279        if(n==0) return(__ieee754_y0(x));
00280        if(n==1) return(sign*__ieee754_y1(x));
00281        if(ix==0x7ff00000) return zero;
00282        if(ix>=0x52D00000) { /* x > 2**302 */
00283     /* (x >> n**2) 
00284      *     Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
00285      *     Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
00286      *     Let s=sin(x), c=cos(x), 
00287      *        xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
00288      *
00289      *           n   sin(xn)*sqt2  cos(xn)*sqt2
00290      *        ----------------------------------
00291      *           0    s-c           c+s
00292      *           1   -s-c          -c+s
00293      *           2   -s+c          -c-s
00294      *           3    s+c           c-s
00295      */
00296               switch(n&3) {
00297                   case 0: temp =  fd_sin(x)-fd_cos(x); break;
00298                   case 1: temp = -fd_sin(x)-fd_cos(x); break;
00299                   case 2: temp = -fd_sin(x)+fd_cos(x); break;
00300                   case 3: temp =  fd_sin(x)+fd_cos(x); break;
00301               }
00302               b = invsqrtpi*temp/fd_sqrt(x);
00303        } else {
00304            a = __ieee754_y0(x);
00305            b = __ieee754_y1(x);
00306        /* quit if b is -inf */
00307             u.d = b;
00308            for(i=1;i<n&&(__HI(u) != 0xfff00000);i++){ 
00309               temp = b;
00310               b = ((double)(i+i)/x)*b - a;
00311               a = temp;
00312            }
00313        }
00314        if(sign>0) return b; else return -b;
00315 }