Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 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 |
Brad Fitzpatrick | 5fea2cc | 2016-03-01 23:21:55 +0000 | [diff] [blame] | 13 | // came with this notice. The go code is a simplified |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 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 | ) |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 58 | // special cases |
| 59 | switch { |
Luuk van Dijk | 8dd3de4 | 2012-02-01 16:08:31 +0100 | [diff] [blame] | 60 | case IsNaN(x): |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 61 | return x |
Luuk van Dijk | 8dd3de4 | 2012-02-01 16:08:31 +0100 | [diff] [blame] | 62 | case IsInf(x, 0): |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 63 | return 0 |
| 64 | } |
Charles L. Dorian | 3c3e68b | 2010-04-09 14:37:33 -0700 | [diff] [blame] | 65 | // J(-n, x) = (-1)**n * J(n, x), J(n, -x) = (-1)**n * J(n, x) |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 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) |
Charles L. Dorian | 3c3e68b | 2010-04-09 14:37:33 -0700 | [diff] [blame] | 126 | // J(n,x) = 1/n!*(x/2)**n - ... |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 127 | |
| 128 | if n > 33 { // underflow |
| 129 | b = 0 |
| 130 | } else { |
| 131 | temp := x * 0.5 |
| 132 | b = temp |
Russ Cox | f2b5a07 | 2011-01-19 23:09:00 -0500 | [diff] [blame] | 133 | a := 1.0 |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 134 | for i := 2; i <= n; i++ { |
| 135 | a *= float64(i) // a = n! |
Charles L. Dorian | 3c3e68b | 2010-04-09 14:37:33 -0700 | [diff] [blame] | 136 | b *= temp // b = (x/2)**n |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 137 | } |
| 138 | b /= a |
| 139 | } |
| 140 | } else { |
| 141 | // use backward recurrence |
Charles L. Dorian | 3c3e68b | 2010-04-09 14:37:33 -0700 | [diff] [blame] | 142 | // x x**2 x**2 |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 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 |
Russ Cox | f2b5a07 | 2011-01-19 23:09:00 -0500 | [diff] [blame] | 182 | t := 0.0 |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 183 | for i := 2 * (n + k); i >= m; i -= 2 { |
| 184 | t = 1 / (float64(i)/x - t) |
| 185 | } |
| 186 | a := t |
| 187 | b = 1 |
Charles L. Dorian | 3c3e68b | 2010-04-09 14:37:33 -0700 | [diff] [blame] | 188 | // estimate log((2/x)**n*n!) = n*log(2/x)+n*ln(n) |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 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 |
Rob Pike | 1a13f9b | 2011-09-29 09:54:20 -0700 | [diff] [blame] | 198 | tmp = tmp * Log(Abs(v*tmp)) |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 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 |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 203 | } |
| 204 | } else { |
| 205 | for i := n - 1; i > 0; i-- { |
| 206 | di := float64(i + i) |
| 207 | a, b = b, b*di/x-a |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 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 |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 235 | // special cases |
| 236 | switch { |
Luuk van Dijk | 8dd3de4 | 2012-02-01 16:08:31 +0100 | [diff] [blame] | 237 | case x < 0 || IsNaN(x): |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 238 | return NaN() |
Luuk van Dijk | 8dd3de4 | 2012-02-01 16:08:31 +0100 | [diff] [blame] | 239 | case IsInf(x, 1): |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 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 |
Luuk van Dijk | 8dd3de4 | 2012-02-01 16:08:31 +0100 | [diff] [blame] | 296 | for i := 1; i < n && !IsInf(b, -1); i++ { |
Charles L. Dorian | 1ec91c8 | 2010-03-26 14:09:39 -0700 | [diff] [blame] | 297 | a, b = b, (float64(i+i)/x)*b-a |
| 298 | } |
| 299 | } |
| 300 | if sign { |
| 301 | return -b |
| 302 | } |
| 303 | return b |
| 304 | } |