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         // TODO(rsc): Remove manual inlining of IsNaN, IsInf
  59         // when compiler does it for us
  60         // special cases
  61         switch {
  62         case x != x: // IsNaN(x)
  63                 return x
  64         case x < -MaxFloat64 || x > MaxFloat64: // IsInf(x, 0):
  65                 return 0
  66         }
  67         // J(-n, x) = (-1)**n * J(n, x), J(n, -x) = (-1)**n * J(n, x)
  68         // Thus, J(-n, x) = J(n, -x)
  69
  70         if n == 0 {
  71                 return J0(x)
  72         }
  73         if x == 0 {
  74                 return 0
  75         }
  76         if n < 0 {
  77                 n, x = -n, -x
  78         }
  79         if n == 1 {
  80                 return J1(x)
  81         }
  82         sign := false
  83         if x < 0 {
  84                 x = -x
  85                 if n&1 == 1 {
  86                         sign = true // odd n and negative x
  87                 }
  88         }
  89         var b float64
  90         if float64(n) <= x {
  91                 // Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x)
  92                 if x >= Two302 { // x > 2**302
  93
  94                         // (x >> n**2)
  95                         //          Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
  96                         //          Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
  97                         //          Let s=sin(x), c=cos(x),
  98                         //              xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
  99                         //
 100                         //                 n    sin(xn)*sqt2    cos(xn)*sqt2
 101                         //              ----------------------------------
 102                         //                 0     s-c             c+s
 103                         //                 1    -s-c            -c+s
 104                         //                 2    -s+c            -c-s
 105                         //                 3     s+c             c-s
 106
 107                         var temp float64
 108                         switch n & 3 {
 109                         case 0:
 110                                 temp = Cos(x) + Sin(x)
 111                         case 1:
 112                                 temp = -Cos(x) + Sin(x)
 113                         case 2:
 114                                 temp = -Cos(x) - Sin(x)
 115                         case 3:
 116                                 temp = Cos(x) - Sin(x)
 117                         }
 118                         b = (1 / SqrtPi) * temp / Sqrt(x)
 119                 } else {
 120                         b = J1(x)
 121                         for i, a := 1, J0(x); i < n; i++ {
 122                                 a, b = b, b*(float64(i+i)/x)-a // avoid underflow
 123                         }
 124                 }
 125         } else {
 126                 if x < TwoM29 { // x < 2**-29
 127                         // x is tiny, return the first Taylor expansion of J(n,x)
 128                         // J(n,x) = 1/n!*(x/2)**n  - ...
 129
 130                         if n > 33 { // underflow
 131                                 b = 0
 132                         } else {
 133                                 temp := x * 0.5
 134                                 b = temp
 135                                 a := 1.0
 136                                 for i := 2; i <= n; i++ {
 137                                         a *= float64(i) // a = n!
 138                                         b *= temp       // b = (x/2)**n
 139                                 }
 140                                 b /= a
 141                         }
 142                 } else {
 143                         // use backward recurrence
 144                         //                      x      x**2      x**2
 145                         //  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
 146                         //                      2n  - 2(n+1) - 2(n+2)
 147                         //
 148                         //                      1      1        1
 149                         //  (for large x)   =  ----  ------   ------   .....
 150                         //                      2n   2(n+1)   2(n+2)
 151                         //                      -- - ------ - ------ -
 152                         //                       x     x         x
 153                         //
 154                         // Let w = 2n/x and h=2/x, then the above quotient
 155                         // is equal to the continued fraction:
 156                         //                  1
 157                         //      = -----------------------
 158                         //                     1
 159                         //         w - -----------------
 160                         //                        1
 161                         //              w+h - ---------
 162                         //                     w+2h - ...
 163                         //
 164                         // To determine how many terms needed, let
 165                         // Q(0) = w, Q(1) = w(w+h) - 1,
 166                         // Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
 167                         // When Q(k) > 1e4      good for single
 168                         // When Q(k) > 1e9      good for double
 169                         // When Q(k) > 1e17     good for quadruple
 170
 171                         // determine k
 172                         w := float64(n+n) / x
 173                         h := 2 / x
 174                         q0 := w
 175                         z := w + h
 176                         q1 := w*z - 1
 177                         k := 1
 178                         for q1 < 1e9 {
 179                                 k += 1
 180                                 z += h
 181                                 q0, q1 = q1, z*q1-q0
 182                         }
 183                         m := n + n
 184                         t := 0.0
 185                         for i := 2 * (n + k); i >= m; i -= 2 {
 186                                 t = 1 / (float64(i)/x - t)
 187                         }
 188                         a := t
 189                         b = 1
 190                         //  estimate log((2/x)**n*n!) = n*log(2/x)+n*ln(n)
 191                         //  Hence, if n*(log(2n/x)) > ...
 192                         //  single 8.8722839355e+01
 193                         //  double 7.09782712893383973096e+02
 194                         //  long double 1.1356523406294143949491931077970765006170e+04
 195                         //  then recurrent value may overflow and the result is
 196                         //  likely underflow to zero
 197
 198                         tmp := float64(n)
 199                         v := 2 / x
 200                         tmp = tmp * Log(Abs(v*tmp))
 201                         if tmp < 7.09782712893383973096e+02 {
 202                                 for i := n - 1; i > 0; i-- {
 203                                         di := float64(i + i)
 204                                         a, b = b, b*di/x-a
 205                                         di -= 2
 206                                 }
 207                         } else {
 208                                 for i := n - 1; i > 0; i-- {
 209                                         di := float64(i + i)
 210                                         a, b = b, b*di/x-a
 211                                         di -= 2
 212                                         // scale b to avoid spurious overflow
 213                                         if b > 1e100 {
 214                                                 a /= b
 215                                                 t /= b
 216                                                 b = 1
 217                                         }
 218                                 }
 219                         }
 220                         b = t * J0(x) / b
 221                 }
 222         }
 223         if sign {
 224                 return -b
 225         }
 226         return b
 227 }
 228
 229 // Yn returns the order-n Bessel function of the second kind.
 230 //
 231 // Special cases are:
 232 //      Yn(n, +Inf) = 0
 233 //      Yn(n > 0, 0) = -Inf
 234 //      Yn(n < 0, 0) = +Inf if n is odd, -Inf if n is even
 235 //      Y1(n, x < 0) = NaN
 236 //      Y1(n, NaN) = NaN
 237 func Yn(n int, x float64) float64 {
 238         const Two302 = 1 << 302 // 2**302 0x52D0000000000000
 239         // TODO(rsc): Remove manual inlining of IsNaN, IsInf
 240         // when compiler does it for us
 241         // special cases
 242         switch {
 243         case x < 0 || x != x: // x < 0 || IsNaN(x):
 244                 return NaN()
 245         case x > MaxFloat64: // IsInf(x, 1)
 246                 return 0
 247         }
 248
 249         if n == 0 {
 250                 return Y0(x)
 251         }
 252         if x == 0 {
 253                 if n < 0 && n&1 == 1 {
 254                         return Inf(1)
 255                 }
 256                 return Inf(-1)
 257         }
 258         sign := false
 259         if n < 0 {
 260                 n = -n
 261                 if n&1 == 1 {
 262                         sign = true // sign true if n < 0 && |n| odd
 263                 }
 264         }
 265         if n == 1 {
 266                 if sign {
 267                         return -Y1(x)
 268                 }
 269                 return Y1(x)
 270         }
 271         var b float64
 272         if x >= Two302 { // x > 2**302
 273                 // (x >> n**2)
 274                 //          Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 275                 //          Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 276                 //          Let s=sin(x), c=cos(x),
 277                 //              xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
 278                 //
 279                 //                 n    sin(xn)*sqt2    cos(xn)*sqt2
 280                 //              ----------------------------------
 281                 //                 0     s-c             c+s
 282                 //                 1    -s-c            -c+s
 283                 //                 2    -s+c            -c-s
 284                 //                 3     s+c             c-s
 285
 286                 var temp float64
 287                 switch n & 3 {
 288                 case 0:
 289                         temp = Sin(x) - Cos(x)
 290                 case 1:
 291                         temp = -Sin(x) - Cos(x)
 292                 case 2:
 293                         temp = -Sin(x) + Cos(x)
 294                 case 3:
 295                         temp = Sin(x) + Cos(x)
 296                 }
 297                 b = (1 / SqrtPi) * temp / Sqrt(x)
 298         } else {
 299                 a := Y0(x)
 300                 b = Y1(x)
 301                 // quit if b is -inf
 302                 for i := 1; i < n && b >= -MaxFloat64; i++ { // for i := 1; i < n && !IsInf(b, -1); i++ {
 303                         a, b = b, (float64(i+i)/x)*b-a
 304                 }
 305         }
 306         if sign {
 307                 return -b
 308         }
 309         return b
 310 }