vendor/libm/src/math/jn.rs

   1 /* origin: FreeBSD /usr/src/lib/msun/src/e_jn.c */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunSoft, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12 /*
  13  * jn(n, x), yn(n, x)
  14  * floating point Bessel's function of the 1st and 2nd kind
  15  * of order n
  16  *
  17  * Special cases:
  18  *      y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
  19  *      y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
  20  * Note 2. About jn(n,x), yn(n,x)
  21  *      For n=0, j0(x) is called,
  22  *      for n=1, j1(x) is called,
  23  *      for n<=x, forward recursion is used starting
  24  *      from values of j0(x) and j1(x).
  25  *      for n>x, a continued fraction approximation to
  26  *      j(n,x)/j(n-1,x) is evaluated and then backward
  27  *      recursion is used starting from a supposed value
  28  *      for j(n,x). The resulting value of j(0,x) is
  29  *      compared with the actual value to correct the
  30  *      supposed value of j(n,x).
  31  *
  32  *      yn(n,x) is similar in all respects, except
  33  *      that forward recursion is used for all
  34  *      values of n>1.
  35  */
  36
  37 use super::{cos, fabs, get_high_word, get_low_word, j0, j1, log, sin, sqrt, y0, y1};
  38
  39 const INVSQRTPI: f64 = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
  40
  41 pub fn jn(n: i32, mut x: f64) -> f64 {
  42     let mut ix: u32;
  43     let lx: u32;
  44     let nm1: i32;
  45     let mut i: i32;
  46     let mut sign: bool;
  47     let mut a: f64;
  48     let mut b: f64;
  49     let mut temp: f64;
  50
  51     ix = get_high_word(x);
  52     lx = get_low_word(x);
  53     sign = (ix >> 31) != 0;
  54     ix &= 0x7fffffff;
  55
  56     // -lx == !lx + 1
  57     if (ix | (lx | ((!lx).wrapping_add(1))) >> 31) > 0x7ff00000 {
  58         /* nan */
  59         return x;
  60     }
  61
  62     /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
  63      * Thus, J(-n,x) = J(n,-x)
  64      */
  65     /* nm1 = |n|-1 is used instead of |n| to handle n==INT_MIN */
  66     if n == 0 {
  67         return j0(x);
  68     }
  69     if n < 0 {
  70         nm1 = -(n + 1);
  71         x = -x;
  72         sign = !sign;
  73     } else {
  74         nm1 = n - 1;
  75     }
  76     if nm1 == 0 {
  77         return j1(x);
  78     }
  79
  80     sign &= (n & 1) != 0; /* even n: 0, odd n: signbit(x) */
  81     x = fabs(x);
  82     if (ix | lx) == 0 || ix == 0x7ff00000 {
  83         /* if x is 0 or inf */
  84         b = 0.0;
  85     } else if (nm1 as f64) < x {
  86         /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
  87         if ix >= 0x52d00000 {
  88             /* x > 2**302 */
  89             /* (x >> n**2)
  90              *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
  91              *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
  92              *      Let s=sin(x), c=cos(x),
  93              *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
  94              *
  95              *             n    sin(xn)*sqt2    cos(xn)*sqt2
  96              *          ----------------------------------
  97              *             0     s-c             c+s
  98              *             1    -s-c            -c+s
  99              *             2    -s+c            -c-s
 100              *             3     s+c             c-s
 101              */
 102             temp = match nm1 & 3 {
 103                 0 => -cos(x) + sin(x),
 104                 1 => -cos(x) - sin(x),
 105                 2 => cos(x) - sin(x),
 106                 3 | _ => cos(x) + sin(x),
 107             };
 108             b = INVSQRTPI * temp / sqrt(x);
 109         } else {
 110             a = j0(x);
 111             b = j1(x);
 112             i = 0;
 113             while i < nm1 {
 114                 i += 1;
 115                 temp = b;
 116                 b = b * (2.0 * (i as f64) / x) - a; /* avoid underflow */
 117                 a = temp;
 118             }
 119         }
 120     } else {
 121         if ix < 0x3e100000 {
 122             /* x < 2**-29 */
 123             /* x is tiny, return the first Taylor expansion of J(n,x)
 124              * J(n,x) = 1/n!*(x/2)^n  - ...
 125              */
 126             if nm1 > 32 {
 127                 /* underflow */
 128                 b = 0.0;
 129             } else {
 130                 temp = x * 0.5;
 131                 b = temp;
 132                 a = 1.0;
 133                 i = 2;
 134                 while i <= nm1 + 1 {
 135                     a *= i as f64; /* a = n! */
 136                     b *= temp; /* b = (x/2)^n */
 137                     i += 1;
 138                 }
 139                 b = b / a;
 140             }
 141         } else {
 142             /* use backward recurrence */
 143             /*                      x      x^2      x^2
 144              *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
 145              *                      2n  - 2(n+1) - 2(n+2)
 146              *
 147              *                      1      1        1
 148              *  (for large x)   =  ----  ------   ------   .....
 149              *                      2n   2(n+1)   2(n+2)
 150              *                      -- - ------ - ------ -
 151              *                       x     x         x
 152              *
 153              * Let w = 2n/x and h=2/x, then the above quotient
 154              * is equal to the continued fraction:
 155              *                  1
 156              *      = -----------------------
 157              *                     1
 158              *         w - -----------------
 159              *                        1
 160              *              w+h - ---------
 161              *                     w+2h - ...
 162              *
 163              * To determine how many terms needed, let
 164              * Q(0) = w, Q(1) = w(w+h) - 1,
 165              * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
 166              * When Q(k) > 1e4      good for single
 167              * When Q(k) > 1e9      good for double
 168              * When Q(k) > 1e17     good for quadruple
 169              */
 170             /* determine k */
 171             let mut t: f64;
 172             let mut q0: f64;
 173             let mut q1: f64;
 174             let mut w: f64;
 175             let h: f64;
 176             let mut z: f64;
 177             let mut tmp: f64;
 178             let nf: f64;
 179
 180             let mut k: i32;
 181
 182             nf = (nm1 as f64) + 1.0;
 183             w = 2.0 * nf / x;
 184             h = 2.0 / x;
 185             z = w + h;
 186             q0 = w;
 187             q1 = w * z - 1.0;
 188             k = 1;
 189             while q1 < 1.0e9 {
 190                 k += 1;
 191                 z += h;
 192                 tmp = z * q1 - q0;
 193                 q0 = q1;
 194                 q1 = tmp;
 195             }
 196             t = 0.0;
 197             i = k;
 198             while i >= 0 {
 199                 t = 1.0 / (2.0 * ((i as f64) + nf) / x - t);
 200                 i -= 1;
 201             }
 202             a = t;
 203             b = 1.0;
 204             /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
 205              *  Hence, if n*(log(2n/x)) > ...
 206              *  single 8.8722839355e+01
 207              *  double 7.09782712893383973096e+02
 208              *  long double 1.1356523406294143949491931077970765006170e+04
 209              *  then recurrent value may overflow and the result is
 210              *  likely underflow to zero
 211              */
 212             tmp = nf * log(fabs(w));
 213             if tmp < 7.09782712893383973096e+02 {
 214                 i = nm1;
 215                 while i > 0 {
 216                     temp = b;
 217                     b = b * (2.0 * (i as f64)) / x - a;
 218                     a = temp;
 219                     i -= 1;
 220                 }
 221             } else {
 222                 i = nm1;
 223                 while i > 0 {
 224                     temp = b;
 225                     b = b * (2.0 * (i as f64)) / x - a;
 226                     a = temp;
 227                     /* scale b to avoid spurious overflow */
 228                     let x1p500 = f64::from_bits(0x5f30000000000000); // 0x1p500 == 2^500
 229                     if b > x1p500 {
 230                         a /= b;
 231                         t /= b;
 232                         b = 1.0;
 233                     }
 234                     i -= 1;
 235                 }
 236             }
 237             z = j0(x);
 238             w = j1(x);
 239             if fabs(z) >= fabs(w) {
 240                 b = t * z / b;
 241             } else {
 242                 b = t * w / a;
 243             }
 244         }
 245     }
 246
 247     if sign {
 248         -b
 249     } else {
 250         b
 251     }
 252 }
 253
 254 pub fn yn(n: i32, x: f64) -> f64 {
 255     let mut ix: u32;
 256     let lx: u32;
 257     let mut ib: u32;
 258     let nm1: i32;
 259     let mut sign: bool;
 260     let mut i: i32;
 261     let mut a: f64;
 262     let mut b: f64;
 263     let mut temp: f64;
 264
 265     ix = get_high_word(x);
 266     lx = get_low_word(x);
 267     sign = (ix >> 31) != 0;
 268     ix &= 0x7fffffff;
 269
 270     // -lx == !lx + 1
 271     if (ix | (lx | ((!lx).wrapping_add(1))) >> 31) > 0x7ff00000 {
 272         /* nan */
 273         return x;
 274     }
 275     if sign && (ix | lx) != 0 {
 276         /* x < 0 */
 277         return 0.0 / 0.0;
 278     }
 279     if ix == 0x7ff00000 {
 280         return 0.0;
 281     }
 282
 283     if n == 0 {
 284         return y0(x);
 285     }
 286     if n < 0 {
 287         nm1 = -(n + 1);
 288         sign = (n & 1) != 0;
 289     } else {
 290         nm1 = n - 1;
 291         sign = false;
 292     }
 293     if nm1 == 0 {
 294         if sign {
 295             return -y1(x);
 296         } else {
 297             return y1(x);
 298         }
 299     }
 300
 301     if ix >= 0x52d00000 {
 302         /* x > 2**302 */
 303         /* (x >> n**2)
 304          *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 305          *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 306          *      Let s=sin(x), c=cos(x),
 307          *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 308          *
 309          *             n    sin(xn)*sqt2    cos(xn)*sqt2
 310          *          ----------------------------------
 311          *             0     s-c             c+s
 312          *             1    -s-c            -c+s
 313          *             2    -s+c            -c-s
 314          *             3     s+c             c-s
 315          */
 316         temp = match nm1 & 3 {
 317             0 => -sin(x) - cos(x),
 318             1 => -sin(x) + cos(x),
 319             2 => sin(x) + cos(x),
 320             3 | _ => sin(x) - cos(x),
 321         };
 322         b = INVSQRTPI * temp / sqrt(x);
 323     } else {
 324         a = y0(x);
 325         b = y1(x);
 326         /* quit if b is -inf */
 327         ib = get_high_word(b);
 328         i = 0;
 329         while i < nm1 && ib != 0xfff00000 {
 330             i += 1;
 331             temp = b;
 332             b = (2.0 * (i as f64) / x) * b - a;
 333             ib = get_high_word(b);
 334             a = temp;
 335         }
 336     }
 337
 338     if sign {
 339         -b
 340     } else {
 341         b
 342     }
 343 }