libgo/go/math/jn.go

   1 // Copyright 2010 The Go Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package math
   6
   7 /*
   8         Bessel function of the first and second kinds of order n.
   9 */
  10
  11 // The original C code and the long comment below are
  12 // from FreeBSD's /usr/src/lib/msun/src/e_jn.c and
  13 // came with this notice.  The go code is a simplified
  14 // version of the original C.
  15 //
  16 // ====================================================
  17 // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  18 //
  19 // Developed at SunPro, a Sun Microsystems, Inc. business.
  20 // Permission to use, copy, modify, and distribute this
  21 // software is freely granted, provided that this notice
  22 // is preserved.
  23 // ====================================================
  24 //
  25 // __ieee754_jn(n, x), __ieee754_yn(n, x)
  26 // floating point Bessel's function of the 1st and 2nd kind
  27 // of order n
  28 //
  29 // Special cases:
  30 //      y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
  31 //      y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
  32 // Note 2. About jn(n,x), yn(n,x)
  33 //      For n=0, j0(x) is called,
  34 //      for n=1, j1(x) is called,
  35 //      for n<x, forward recursion is used starting
  36 //      from values of j0(x) and j1(x).
  37 //      for n>x, a continued fraction approximation to
  38 //      j(n,x)/j(n-1,x) is evaluated and then backward
  39 //      recursion is used starting from a supposed value
  40 //      for j(n,x). The resulting value of j(0,x) is
  41 //      compared with the actual value to correct the
  42 //      supposed value of j(n,x).
  43 //
  44 //      yn(n,x) is similar in all respects, except
  45 //      that forward recursion is used for all
  46 //      values of n>1.
  47
  48 // Jn returns the order-n Bessel function of the first kind.
  49 //
  50 // Special cases are:
  51 //      Jn(n, ±Inf) = 0
  52 //      Jn(n, NaN) = NaN
  53 func Jn(n int, x float64) float64 {
  54         const (
  55                 TwoM29 = 1.0 / (1 << 29) // 2**-29 0x3e10000000000000
  56                 Two302 = 1 << 302        // 2**302 0x52D0000000000000
  57         )
  58         // special cases
  59         switch {
  60         case IsNaN(x):
  61                 return x
  62         case IsInf(x, 0):
  63                 return 0
  64         }
  65         // J(-n, x) = (-1)**n * J(n, x), J(n, -x) = (-1)**n * J(n, x)
  66         // Thus, J(-n, x) = J(n, -x)
  67
  68         if n == 0 {
  69                 return J0(x)
  70         }
  71         if x == 0 {
  72                 return 0
  73         }
  74         if n < 0 {
  75                 n, x = -n, -x
  76         }
  77         if n == 1 {
  78                 return J1(x)
  79         }
  80         sign := false
  81         if x < 0 {
  82                 x = -x
  83                 if n&1 == 1 {
  84                         sign = true // odd n and negative x
  85                 }
  86         }
  87         var b float64
  88         if float64(n) <= x {
  89                 // Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x)
  90                 if x >= Two302 { // x > 2**302
  91
  92                         // (x >> n**2)
  93                         //          Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
  94                         //          Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
  95                         //          Let s=sin(x), c=cos(x),
  96                         //              xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
  97                         //
  98                         //                 n    sin(xn)*sqt2    cos(xn)*sqt2
  99                         //              ----------------------------------
 100                         //                 0     s-c             c+s
 101                         //                 1    -s-c            -c+s
 102                         //                 2    -s+c            -c-s
 103                         //                 3     s+c             c-s
 104
 105                         var temp float64
 106                         switch n & 3 {
 107                         case 0:
 108                                 temp = Cos(x) + Sin(x)
 109                         case 1:
 110                                 temp = -Cos(x) + Sin(x)
 111                         case 2:
 112                                 temp = -Cos(x) - Sin(x)
 113                         case 3:
 114                                 temp = Cos(x) - Sin(x)
 115                         }
 116                         b = (1 / SqrtPi) * temp / Sqrt(x)
 117                 } else {
 118                         b = J1(x)
 119                         for i, a := 1, J0(x); i < n; i++ {
 120                                 a, b = b, b*(float64(i+i)/x)-a // avoid underflow
 121                         }
 122                 }
 123         } else {
 124                 if x < TwoM29 { // x < 2**-29
 125                         // x is tiny, return the first Taylor expansion of J(n,x)
 126                         // J(n,x) = 1/n!*(x/2)**n  - ...
 127
 128                         if n > 33 { // underflow
 129                                 b = 0
 130                         } else {
 131                                 temp := x * 0.5
 132                                 b = temp
 133                                 a := 1.0
 134                                 for i := 2; i <= n; i++ {
 135                                         a *= float64(i) // a = n!
 136                                         b *= temp       // b = (x/2)**n
 137                                 }
 138                                 b /= a
 139                         }
 140                 } else {
 141                         // use backward recurrence
 142                         //                      x      x**2      x**2
 143                         //  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
 144                         //                      2n  - 2(n+1) - 2(n+2)
 145                         //
 146                         //                      1      1        1
 147                         //  (for large x)   =  ----  ------   ------   .....
 148                         //                      2n   2(n+1)   2(n+2)
 149                         //                      -- - ------ - ------ -
 150                         //                       x     x         x
 151                         //
 152                         // Let w = 2n/x and h=2/x, then the above quotient
 153                         // is equal to the continued fraction:
 154                         //                  1
 155                         //      = -----------------------
 156                         //                     1
 157                         //         w - -----------------
 158                         //                        1
 159                         //              w+h - ---------
 160                         //                     w+2h - ...
 161                         //
 162                         // To determine how many terms needed, let
 163                         // Q(0) = w, Q(1) = w(w+h) - 1,
 164                         // Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
 165                         // When Q(k) > 1e4      good for single
 166                         // When Q(k) > 1e9      good for double
 167                         // When Q(k) > 1e17     good for quadruple
 168
 169                         // determine k
 170                         w := float64(n+n) / x
 171                         h := 2 / x
 172                         q0 := w
 173                         z := w + h
 174                         q1 := w*z - 1
 175                         k := 1
 176                         for q1 < 1e9 {
 177                                 k += 1
 178                                 z += h
 179                                 q0, q1 = q1, z*q1-q0
 180                         }
 181                         m := n + n
 182                         t := 0.0
 183                         for i := 2 * (n + k); i >= m; i -= 2 {
 184                                 t = 1 / (float64(i)/x - t)
 185                         }
 186                         a := t
 187                         b = 1
 188                         //  estimate log((2/x)**n*n!) = n*log(2/x)+n*ln(n)
 189                         //  Hence, if n*(log(2n/x)) > ...
 190                         //  single 8.8722839355e+01
 191                         //  double 7.09782712893383973096e+02
 192                         //  long double 1.1356523406294143949491931077970765006170e+04
 193                         //  then recurrent value may overflow and the result is
 194                         //  likely underflow to zero
 195
 196                         tmp := float64(n)
 197                         v := 2 / x
 198                         tmp = tmp * Log(Abs(v*tmp))
 199                         if tmp < 7.09782712893383973096e+02 {
 200                                 for i := n - 1; i > 0; i-- {
 201                                         di := float64(i + i)
 202                                         a, b = b, b*di/x-a
 203                                 }
 204                         } else {
 205                                 for i := n - 1; i > 0; i-- {
 206                                         di := float64(i + i)
 207                                         a, b = b, b*di/x-a
 208                                         // scale b to avoid spurious overflow
 209                                         if b > 1e100 {
 210                                                 a /= b
 211                                                 t /= b
 212                                                 b = 1
 213                                         }
 214                                 }
 215                         }
 216                         b = t * J0(x) / b
 217                 }
 218         }
 219         if sign {
 220                 return -b
 221         }
 222         return b
 223 }
 224
 225 // Yn returns the order-n Bessel function of the second kind.
 226 //
 227 // Special cases are:
 228 //      Yn(n, +Inf) = 0
 229 //      Yn(n > 0, 0) = -Inf
 230 //      Yn(n < 0, 0) = +Inf if n is odd, -Inf if n is even
 231 //      Y1(n, x < 0) = NaN
 232 //      Y1(n, NaN) = NaN
 233 func Yn(n int, x float64) float64 {
 234         const Two302 = 1 << 302 // 2**302 0x52D0000000000000
 235         // special cases
 236         switch {
 237         case x < 0 || IsNaN(x):
 238                 return NaN()
 239         case IsInf(x, 1):
 240                 return 0
 241         }
 242
 243         if n == 0 {
 244                 return Y0(x)
 245         }
 246         if x == 0 {
 247                 if n < 0 && n&1 == 1 {
 248                         return Inf(1)
 249                 }
 250                 return Inf(-1)
 251         }
 252         sign := false
 253         if n < 0 {
 254                 n = -n
 255                 if n&1 == 1 {
 256                         sign = true // sign true if n < 0 && |n| odd
 257                 }
 258         }
 259         if n == 1 {
 260                 if sign {
 261                         return -Y1(x)
 262                 }
 263                 return Y1(x)
 264         }
 265         var b float64
 266         if x >= Two302 { // x > 2**302
 267                 // (x >> n**2)
 268                 //          Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 269                 //          Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 270                 //          Let s=sin(x), c=cos(x),
 271                 //              xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 272                 //
 273                 //                 n    sin(xn)*sqt2    cos(xn)*sqt2
 274                 //              ----------------------------------
 275                 //                 0     s-c             c+s
 276                 //                 1    -s-c            -c+s
 277                 //                 2    -s+c            -c-s
 278                 //                 3     s+c             c-s
 279
 280                 var temp float64
 281                 switch n & 3 {
 282                 case 0:
 283                         temp = Sin(x) - Cos(x)
 284                 case 1:
 285                         temp = -Sin(x) - Cos(x)
 286                 case 2:
 287                         temp = -Sin(x) + Cos(x)
 288                 case 3:
 289                         temp = Sin(x) + Cos(x)
 290                 }
 291                 b = (1 / SqrtPi) * temp / Sqrt(x)
 292         } else {
 293                 a := Y0(x)
 294                 b = Y1(x)
 295                 // quit if b is -inf
 296                 for i := 1; i < n && !IsInf(b, -1); i++ {
 297                         a, b = b, (float64(i+i)/x)*b-a
 298                 }
 299         }
 300         if sign {
 301                 return -b
 302         }
 303         return b
 304 }