src/dom/js/fdlibm/e_jn.c

   1 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
   2  *
   3  * ***** BEGIN LICENSE BLOCK *****
   4  * Version: MPL 1.1/GPL 2.0/LGPL 2.1
   5  *
   6  * The contents of this file are subject to the Mozilla Public License Version
   7  * 1.1 (the "License"); you may not use this file except in compliance with
   8  * the License. You may obtain a copy of the License at
   9  * http://www.mozilla.org/MPL/
  10  *
  11  * Software distributed under the License is distributed on an "AS IS" basis,
  12  * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
  13  * for the specific language governing rights and limitations under the
  14  * License.
  15  *
  16  * The Original Code is Mozilla Communicator client code, released
  17  * March 31, 1998.
  18  *
  19  * The Initial Developer of the Original Code is
  20  * Sun Microsystems, Inc.
  21  * Portions created by the Initial Developer are Copyright (C) 1998
  22  * the Initial Developer. All Rights Reserved.
  23  *
  24  * Contributor(s):
  25  *
  26  * Alternatively, the contents of this file may be used under the terms of
  27  * either of the GNU General Public License Version 2 or later (the "GPL"),
  28  * or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
  29  * in which case the provisions of the GPL or the LGPL are applicable instead
  30  * of those above. If you wish to allow use of your version of this file only
  31  * under the terms of either the GPL or the LGPL, and not to allow others to
  32  * use your version of this file under the terms of the MPL, indicate your
  33  * decision by deleting the provisions above and replace them with the notice
  34  * and other provisions required by the GPL or the LGPL. If you do not delete
  35  * the provisions above, a recipient may use your version of this file under
  36  * the terms of any one of the MPL, the GPL or the LGPL.
  37  *
  38  * ***** END LICENSE BLOCK ***** */
  39
  40 /* @(#)e_jn.c 1.4 95/01/18 */
  41 /*
  42  * ====================================================
  43  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  44  *
  45  * Developed at SunSoft, a Sun Microsystems, Inc. business.
  46  * Permission to use, copy, modify, and distribute this
  47  * software is freely granted, provided that this notice
  48  * is preserved.
  49  * ====================================================
  50  */
  51
  52 /*
  53  * __ieee754_jn(n, x), __ieee754_yn(n, x)
  54  * floating point Bessel's function of the 1st and 2nd kind
  55  * of order n
  56  *
  57  * Special cases:
  58  *      y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
  59  *      y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
  60  * Note 2. About jn(n,x), yn(n,x)
  61  *      For n=0, j0(x) is called,
  62  *      for n=1, j1(x) is called,
  63  *      for n<x, forward recursion us used starting
  64  *      from values of j0(x) and j1(x).
  65  *      for n>x, a continued fraction approximation to
  66  *      j(n,x)/j(n-1,x) is evaluated and then backward
  67  *      recursion is used starting from a supposed value
  68  *      for j(n,x). The resulting value of j(0,x) is
  69  *      compared with the actual value to correct the
  70  *      supposed value of j(n,x).
  71  *
  72  *      yn(n,x) is similar in all respects, except
  73  *      that forward recursion is used for all
  74  *      values of n>1.
  75  *
  76  */
  77
  78 #include "fdlibm.h"
  79
  80 #ifdef __STDC__
  81 static const double
  82 #else
  83 static double
  84 #endif
  85 invsqrtpi=  5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
  86 two   =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
  87 one   =  1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
  88
  89 static double zero  =  0.00000000000000000000e+00;
  90
  91 #ifdef __STDC__
  92         double __ieee754_jn(int n, double x)
  93 #else
  94         double __ieee754_jn(n,x)
  95         int n; double x;
  96 #endif
  97 {
  98         fd_twoints u;
  99         int i,hx,ix,lx, sgn;
 100         double a, b, temp, di;
 101         double z, w;
 102
 103     /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
 104      * Thus, J(-n,x) = J(n,-x)
 105      */
 106         u.d = x;
 107         hx = __HI(u);
 108         ix = 0x7fffffff&hx;
 109         lx = __LO(u);
 110     /* if J(n,NaN) is NaN */
 111         if((ix|((unsigned)(lx|-lx))>>31)>0x7ff00000) return x+x;
 112         if(n<0){
 113                 n = -n;
 114                 x = -x;
 115                 hx ^= 0x80000000;
 116         }
 117         if(n==0) return(__ieee754_j0(x));
 118         if(n==1) return(__ieee754_j1(x));
 119         sgn = (n&1)&(hx>>31);   /* even n -- 0, odd n -- sign(x) */
 120         x = fd_fabs(x);
 121         if((ix|lx)==0||ix>=0x7ff00000)  /* if x is 0 or inf */
 122             b = zero;
 123         else if((double)n<=x) {
 124                 /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
 125             if(ix>=0x52D00000) { /* x > 2**302 */
 126     /* (x >> n**2)
 127      *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 128      *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 129      *      Let s=sin(x), c=cos(x),
 130      *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 131      *
 132      *             n    sin(xn)*sqt2    cos(xn)*sqt2
 133      *          ----------------------------------
 134      *             0     s-c             c+s
 135      *             1    -s-c            -c+s
 136      *             2    -s+c            -c-s
 137      *             3     s+c             c-s
 138      */
 139                 switch(n&3) {
 140                     case 0: temp =  fd_cos(x)+fd_sin(x); break;
 141                     case 1: temp = -fd_cos(x)+fd_sin(x); break;
 142                     case 2: temp = -fd_cos(x)-fd_sin(x); break;
 143                     case 3: temp =  fd_cos(x)-fd_sin(x); break;
 144                 }
 145                 b = invsqrtpi*temp/fd_sqrt(x);
 146             } else {
 147                 a = __ieee754_j0(x);
 148                 b = __ieee754_j1(x);
 149                 for(i=1;i<n;i++){
 150                     temp = b;
 151                     b = b*((double)(i+i)/x) - a; /* avoid underflow */
 152                     a = temp;
 153                 }
 154             }
 155         } else {
 156             if(ix<0x3e100000) { /* x < 2**-29 */
 157     /* x is tiny, return the first Taylor expansion of J(n,x)
 158      * J(n,x) = 1/n!*(x/2)^n  - ...
 159      */
 160                 if(n>33)        /* underflow */
 161                     b = zero;
 162                 else {
 163                     temp = x*0.5; b = temp;
 164                     for (a=one,i=2;i<=n;i++) {
 165                         a *= (double)i;         /* a = n! */
 166                         b *= temp;              /* b = (x/2)^n */
 167                     }
 168                     b = b/a;
 169                 }
 170             } else {
 171                 /* use backward recurrence */
 172                 /*                      x      x^2      x^2
 173                  *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
 174                  *                      2n  - 2(n+1) - 2(n+2)
 175                  *
 176                  *                      1      1        1
 177                  *  (for large x)   =  ----  ------   ------   .....
 178                  *                      2n   2(n+1)   2(n+2)
 179                  *                      -- - ------ - ------ -
 180                  *                       x     x         x
 181                  *
 182                  * Let w = 2n/x and h=2/x, then the above quotient
 183                  * is equal to the continued fraction:
 184                  *                  1
 185                  *      = -----------------------
 186                  *                     1
 187                  *         w - -----------------
 188                  *                        1
 189                  *              w+h - ---------
 190                  *                     w+2h - ...
 191                  *
 192                  * To determine how many terms needed, let
 193                  * Q(0) = w, Q(1) = w(w+h) - 1,
 194                  * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
 195                  * When Q(k) > 1e4      good for single
 196                  * When Q(k) > 1e9      good for double
 197                  * When Q(k) > 1e17     good for quadruple
 198                  */
 199             /* determine k */
 200                 double t,v;
 201                 double q0,q1,h,tmp; int k,m;
 202                 w  = (n+n)/(double)x; h = 2.0/(double)x;
 203                 q0 = w;  z = w+h; q1 = w*z - 1.0; k=1;
 204                 while(q1<1.0e9) {
 205                         k += 1; z += h;
 206                         tmp = z*q1 - q0;
 207                         q0 = q1;
 208                         q1 = tmp;
 209                 }
 210                 m = n+n;
 211                 for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
 212                 a = t;
 213                 b = one;
 214                 /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
 215                  *  Hence, if n*(log(2n/x)) > ...
 216                  *  single 8.8722839355e+01
 217                  *  double 7.09782712893383973096e+02
 218                  *  long double 1.1356523406294143949491931077970765006170e+04
 219                  *  then recurrent value may overflow and the result is
 220                  *  likely underflow to zero
 221                  */
 222                 tmp = n;
 223                 v = two/x;
 224                 tmp = tmp*__ieee754_log(fd_fabs(v*tmp));
 225                 if(tmp<7.09782712893383973096e+02) {
 226                     for(i=n-1,di=(double)(i+i);i>0;i--){
 227                         temp = b;
 228                         b *= di;
 229                         b  = b/x - a;
 230                         a = temp;
 231                         di -= two;
 232                     }
 233                 } else {
 234                     for(i=n-1,di=(double)(i+i);i>0;i--){
 235                         temp = b;
 236                         b *= di;
 237                         b  = b/x - a;
 238                         a = temp;
 239                         di -= two;
 240                     /* scale b to avoid spurious overflow */
 241                         if(b>1e100) {
 242                             a /= b;
 243                             t /= b;
 244                             b  = one;
 245                         }
 246                     }
 247                 }
 248                 b = (t*__ieee754_j0(x)/b);
 249             }
 250         }
 251         if(sgn==1) return -b; else return b;
 252 }
 253
 254 #ifdef __STDC__
 255         double __ieee754_yn(int n, double x)
 256 #else
 257         double __ieee754_yn(n,x)
 258         int n; double x;
 259 #endif
 260 {
 261         fd_twoints u;
 262         int i,hx,ix,lx;
 263         int sign;
 264         double a, b, temp;
 265
 266         u.d = x;
 267         hx = __HI(u);
 268         ix = 0x7fffffff&hx;
 269         lx = __LO(u);
 270     /* if Y(n,NaN) is NaN */
 271         if((ix|((unsigned)(lx|-lx))>>31)>0x7ff00000) return x+x;
 272         if((ix|lx)==0) return -one/zero;
 273         if(hx<0) return zero/zero;
 274         sign = 1;
 275         if(n<0){
 276                 n = -n;
 277                 sign = 1 - ((n&1)<<1);
 278         }
 279         if(n==0) return(__ieee754_y0(x));
 280         if(n==1) return(sign*__ieee754_y1(x));
 281         if(ix==0x7ff00000) return zero;
 282         if(ix>=0x52D00000) { /* x > 2**302 */
 283     /* (x >> n**2)
 284      *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 285      *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 286      *      Let s=sin(x), c=cos(x),
 287      *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 288      *
 289      *             n    sin(xn)*sqt2    cos(xn)*sqt2
 290      *          ----------------------------------
 291      *             0     s-c             c+s
 292      *             1    -s-c            -c+s
 293      *             2    -s+c            -c-s
 294      *             3     s+c             c-s
 295      */
 296                 switch(n&3) {
 297                     case 0: temp =  fd_sin(x)-fd_cos(x); break;
 298                     case 1: temp = -fd_sin(x)-fd_cos(x); break;
 299                     case 2: temp = -fd_sin(x)+fd_cos(x); break;
 300                     case 3: temp =  fd_sin(x)+fd_cos(x); break;
 301                 }
 302                 b = invsqrtpi*temp/fd_sqrt(x);
 303         } else {
 304             a = __ieee754_y0(x);
 305             b = __ieee754_y1(x);
 306         /* quit if b is -inf */
 307             u.d = b;
 308             for(i=1;i<n&&(__HI(u) != 0xfff00000);i++){
 309                 temp = b;
 310                 b = ((double)(i+i)/x)*b - a;
 311                 a = temp;
 312             }
 313         }
 314         if(sign>0) return b; else return -b;
 315 }