123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307 |
- package math
- func Jn(n int, x float64) float64 {
- const (
- TwoM29 = 1.0 / (1 << 29)
- Two302 = 1 << 302
- )
-
- switch {
- case IsNaN(x):
- return x
- case IsInf(x, 0):
- return 0
- }
-
-
- if n == 0 {
- return J0(x)
- }
- if x == 0 {
- return 0
- }
- if n < 0 {
- n, x = -n, -x
- }
- if n == 1 {
- return J1(x)
- }
- sign := false
- if x < 0 {
- x = -x
- if n&1 == 1 {
- sign = true
- }
- }
- var b float64
- if float64(n) <= x {
-
- if x >= Two302 {
-
-
-
-
-
-
-
-
-
-
-
-
- var temp float64
- switch n & 3 {
- case 0:
- temp = Cos(x) + Sin(x)
- case 1:
- temp = -Cos(x) + Sin(x)
- case 2:
- temp = -Cos(x) - Sin(x)
- case 3:
- temp = Cos(x) - Sin(x)
- }
- b = (1 / SqrtPi) * temp / Sqrt(x)
- } else {
- b = J1(x)
- for i, a := 1, J0(x); i < n; i++ {
- a, b = b, b*(float64(i+i)/x)-a
- }
- }
- } else {
- if x < TwoM29 {
-
-
- if n > 33 {
- b = 0
- } else {
- temp := x * 0.5
- b = temp
- a := 1.0
- for i := 2; i <= n; i++ {
- a *= float64(i)
- b *= temp
- }
- b /= a
- }
- } else {
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- w := float64(n+n) / x
- h := 2 / x
- q0 := w
- z := w + h
- q1 := w*z - 1
- k := 1
- for q1 < 1e9 {
- k += 1
- z += h
- q0, q1 = q1, z*q1-q0
- }
- m := n + n
- t := 0.0
- for i := 2 * (n + k); i >= m; i -= 2 {
- t = 1 / (float64(i)/x - t)
- }
- a := t
- b = 1
-
-
-
-
-
-
-
- tmp := float64(n)
- v := 2 / x
- tmp = tmp * Log(Abs(v*tmp))
- if tmp < 7.09782712893383973096e+02 {
- for i := n - 1; i > 0; i-- {
- di := float64(i + i)
- a, b = b, b*di/x-a
- di -= 2
- }
- } else {
- for i := n - 1; i > 0; i-- {
- di := float64(i + i)
- a, b = b, b*di/x-a
- di -= 2
-
- if b > 1e100 {
- a /= b
- t /= b
- b = 1
- }
- }
- }
- b = t * J0(x) / b
- }
- }
- if sign {
- return -b
- }
- return b
- }
- func Yn(n int, x float64) float64 {
- const Two302 = 1 << 302
-
- switch {
- case x < 0 || IsNaN(x):
- return NaN()
- case IsInf(x, 1):
- return 0
- }
- if n == 0 {
- return Y0(x)
- }
- if x == 0 {
- if n < 0 && n&1 == 1 {
- return Inf(1)
- }
- return Inf(-1)
- }
- sign := false
- if n < 0 {
- n = -n
- if n&1 == 1 {
- sign = true
- }
- }
- if n == 1 {
- if sign {
- return -Y1(x)
- }
- return Y1(x)
- }
- var b float64
- if x >= Two302 {
-
-
-
-
-
-
-
-
-
-
-
-
- var temp float64
- switch n & 3 {
- case 0:
- temp = Sin(x) - Cos(x)
- case 1:
- temp = -Sin(x) - Cos(x)
- case 2:
- temp = -Sin(x) + Cos(x)
- case 3:
- temp = Sin(x) + Cos(x)
- }
- b = (1 / SqrtPi) * temp / Sqrt(x)
- } else {
- a := Y0(x)
- b = Y1(x)
-
- for i := 1; i < n && !IsInf(b, -1); i++ {
- a, b = b, (float64(i+i)/x)*b-a
- }
- }
- if sign {
- return -b
- }
- return b
- }
|