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