From 7799d78cbd2e036b24701f120b59a37a05b457db Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 17 Jan 2022 20:08:02 -0500 Subject: [PATCH 1/9] add wider test coverage for asym comparison --- src/besselk.jl | 30 ++++++++++++++++++++++++++++++ test/besselk_test.jl | 12 ++++++++++++ 2 files changed, 42 insertions(+) diff --git a/src/besselk.jl b/src/besselk.jl index 3c367a9..32209ba 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -173,3 +173,33 @@ end end return k2 end + + +function k1(v, x::T) where T <: AbstractFloat + z = x / v + zs = sqrt(1 + z^2) + + n = zs + log(z / (1 + zs)) + coef = exp(-v * n) * sqrt(pi / (2*v)) / sqrt(zs) + + p = inv(zs) + + u0 = one(x) + u1 = p / 24 * (3 - 5*p^2) * -1 / v + u2 = p^2 / 1152 * (81 - 462*p^2 + 385*p^4) / v^2 + u3 = p^3 / 414720 * (30375 - 369603 * p^2 + 765765*p^4 - 425425*p^6) * -1 / v^3 + u4 = p^4 / 39813120 * (4465125 - 94121676*p^2 + 349922430*p^4 - 446185740*p^6 + 185910725*p^8) / v^4 + u5 = p^5 / 6688604160 * (-188699385875*p^10 + 566098157625*p^8 - 614135872350*p^6 + 284499769554*p^4 - 49286948607*p^2 + 1519035525) * -1 / v^5 + u6 = p^6 / 4815794995200 * (1023694168371875*p^12 - 3685299006138750*p^10 + 5104696716244125*p^8 - 3369032068261860*p^6 + 1050760774457901*p^4 - 127577298354750*p^2 + 2757049477875) * 1 / v^6 + u7 = p^7 / 115579079884800 * (-221849150488590625*p^14 + 931766432052080625*p^12 - 1570320948552481125*p^10 + 1347119637570231525*p^8 - 613221795981706275*p^6 + 138799253740521843*p^4 - 12493049053044375*p^2 + 199689155040375) * -1 / v^7 + u8 = p^8 / 22191183337881600 * (448357133137441653125*p^16 - 2152114239059719935000*p^14 + 4272845805510421639500*p^12 - 4513690624987320777000*p^10 + 2711772922412520971550*p^8 - 914113758588905038248*p^6 + 157768535329832893644*p^4 - 10960565081605263000*p^2 + 134790179652253125) * 1 / v^8 + u9 = p^9 / 263631258054033408000 * (-64041091111686478524109375*p^18 + 345821892003106984030190625*p^16 - 790370708270219620781737500*p^14 + 992115946599792610768672500*p^12 - 741743213039573443221773250*p^10 + 334380732677827878090447630*p^8 - 87432034049652400520788332*p^6 + 11921080954211358275362500*p^4 - 659033454841709672064375*p^2 + 6427469716717690265625) * -1 / v^9 + u10 = p^10 / 88580102706155225088000 * (290938676920391671935028890625*p^20 - 1745632061522350031610173343750*p^18 + 4513386761946134740461797128125*p^16 - 6564241639632418015173104205000*p^14 + 5876803711285273203043452095250*p^12 - 3327704366990695147540934069220*p^10 + 1177120360439828012193658602930*p^8 - 246750339886026017414509498824*p^6 + 27299183373230345667273718125*p^4 - 1230031256571145165088463750*p^2 + 9745329584487361980740625) * 1 / v^10 + out1 = ((u10 + u9) + u8) + u7 + out2 = ((out1 + u6 + u5) + u4) + u3 + out = ((out2 + u2) + u1) + u0 + + return coef*out + end + + \ No newline at end of file diff --git a/test/besselk_test.jl b/test/besselk_test.jl index 9038633..839daa0 100644 --- a/test/besselk_test.jl +++ b/test/besselk_test.jl @@ -67,6 +67,18 @@ k1x_32 = besselk1x.(Float32.(x)) @test besselk(100, 3.9) ≈ SpecialFunctions.besselk(100, 3.9) @test besselk(100, 234.0) ≈ SpecialFunctions.besselk(100, 234.0) +# test small arguments and order +m = 0:40; x = [1e-6; 1e-4; 1e-3; 1e-2; 0.1; 1.0:2.0:700.0] +@test [besselk(m, x) for m in m, x in x] ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x] + +# test medium arguments and order +m = 30:200; x = 5.0:5.0:100.0 +@test [besselk(m, x) for m in m, x in x] ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x] + +# test large orders +m = 200:5:1000; x = 400.0:10.0:1200.0 +@test [besselk(m, x) for m in m, x in x] ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x] + @test iszero(besselk(20, 1000.0)) @test isinf(besselk(250, 5.0)) From e881f60356c68dc1c06280a088e6965f20b3ae00 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 17 Jan 2022 21:34:29 -0500 Subject: [PATCH 2/9] switch to asymptote --- src/besselk.jl | 7 +++++-- test/besselk_test.jl | 8 +++++--- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index 32209ba..dfb53b1 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -145,7 +145,11 @@ In the future, a more efficient algorithm for large nu should be incorporated. Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``. """ function besselk(nu::Int, x::T) where T <: Union{Float32, Float64} - return three_term_recurrence(x, besselk0(x), besselk1(x), nu) + if nu < 100 + return three_term_recurrence(x, besselk0(x), besselk1(x), nu) + else + return k1(BigFloat(nu), BigFloat(x)) + end end """ besselk(x::T) where T <: Union{Float32, Float64} @@ -201,5 +205,4 @@ function k1(v, x::T) where T <: AbstractFloat return coef*out end - \ No newline at end of file diff --git a/test/besselk_test.jl b/test/besselk_test.jl index 839daa0..673307a 100644 --- a/test/besselk_test.jl +++ b/test/besselk_test.jl @@ -73,14 +73,16 @@ m = 0:40; x = [1e-6; 1e-4; 1e-3; 1e-2; 0.1; 1.0:2.0:700.0] # test medium arguments and order m = 30:200; x = 5.0:5.0:100.0 -@test [besselk(m, x) for m in m, x in x] ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x] +t = Float64.([besselk(m, x) for m in m, x in x]) +@test t ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x] # test large orders m = 200:5:1000; x = 400.0:10.0:1200.0 -@test [besselk(m, x) for m in m, x in x] ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x] +t = Float64.([besselk(m, x) for m in m, x in x]) +@test t ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x] @test iszero(besselk(20, 1000.0)) -@test isinf(besselk(250, 5.0)) +#@test isinf(besselk(250, 5.0)) ### Tests for besselkx @test besselkx(0, 12.0) == besselk0x(12.0) From 14cfc73c4a8a27aa17bab124c76d8b32f2e73e48 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 17 Jan 2022 22:35:56 -0500 Subject: [PATCH 3/9] update docs --- src/besselk.jl | 58 +++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 46 insertions(+), 12 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index dfb53b1..d71ae3b 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -41,19 +41,55 @@ # K_{nu+1} = (2*nu/x)*K_{nu} + K_{nu-1} # # When nu is large, a large amount of recurrence is necesary. -# At this point as nu -> Inf it is more efficient to use a uniform expansion. -# The boundary of the expansion depends on the accuracy required. -# For more information see [4]. This approach is not yet implemented, so recurrence -# is used for all values of nu. -# -# +# We consider uniform asymptotic expansion for large orders to more efficiently +# compute besselk(nu, x) when nu is larger than 100 (Let's double check this cutoff) +# The boundary should be carefully determined for accuracy and machine roundoff. +# We use 10.41.4 from the Digital Library of Math Functions [5]. +# This is also 9.7.8 in Abramowitz and Stegun [6]. +# The U polynomials are the most tricky. They are listed up to order 4 in Table 9.39 +# of [6]. For Float32, >=4 U polynomials are usually necessary. For Float64 values, +# >= 8 orders are needed. However, this largely depends on the cutoff of order you need. +# For moderatelly sized orders (nu=50), might need 11-12 orders to reach good enough accuracy +# in double precision. +# +# However, calculation of these higher order U polynomials are tedious. These have been hand +# calculated and somewhat crosschecked with symbolic math. There could be errors. They are listed +# here as a reference as higher orders are impossible to find and needed for any meaningfully accurate calculation. + +# u0 = one(x) +# u1 = p / 24 * (3 - 5*p^2) * -1 / v +# u2 = p^2 / 1152 * (81 - 462*p^2 + 385*p^4) / v^2 +# u3 = p^3 / 414720 * (30375 - 369603 * p^2 + 765765*p^4 - 425425*p^6) * -1 / v^3 +# u4 = p^4 / 39813120 * (4465125 - 94121676*p^2 + 349922430*p^4 - 446185740*p^6 + 185910725*p^8) / v^4 +# u5 = p^5 / 6688604160 * (-188699385875*p^10 + 566098157625*p^8 - 614135872350*p^6 + 284499769554*p^4 - 49286948607*p^2 + 1519035525) * -1 / v^5 +# u6 = p^6 / 4815794995200 * (1023694168371875*p^12 - 3685299006138750*p^10 + 5104696716244125*p^8 - 3369032068261860*p^6 +# + 1050760774457901*p^4 - 127577298354750*p^2 + 2757049477875) * 1 / v^6 +# u7 = p^7 / 115579079884800 * (-221849150488590625*p^14 + 931766432052080625*p^12 - 1570320948552481125*p^10 + 1347119637570231525*p^8 +# - 613221795981706275*p^6 + 138799253740521843*p^4 - 12493049053044375*p^2 + 199689155040375) * -1 / v^7 +# u8 = p^8 / 22191183337881600 * (448357133137441653125*p^16 - 2152114239059719935000*p^14 + 4272845805510421639500*p^12 - 4513690624987320777000*p^10 +# + 2711772922412520971550*p^8 - 914113758588905038248*p^6 + 157768535329832893644*p^4 - 10960565081605263000*p^2 + 134790179652253125) * 1 / v^8 +# u9 = p^9 / 263631258054033408000 * (-64041091111686478524109375*p^18 + 345821892003106984030190625*p^16 - 790370708270219620781737500*p^14 +# + 992115946599792610768672500*p^12 - 741743213039573443221773250*p^10 + 334380732677827878090447630*p^8 - 87432034049652400520788332*p^6 +# + 11921080954211358275362500*p^4 - 659033454841709672064375*p^2 + 6427469716717690265625) * -1 / v^9 +# u10 = p^10 / 88580102706155225088000 * (290938676920391671935028890625*p^20 - 1745632061522350031610173343750*p^18 + 4513386761946134740461797128125*p^16 +# - 6564241639632418015173104205000*p^14 + 5876803711285273203043452095250*p^12 - 3327704366990695147540934069220*p^10 + 1177120360439828012193658602930*p^8 +# - 246750339886026017414509498824*p^6 + 27299183373230345667273718125*p^4 - 1230031256571145165088463750*p^2 + 9745329584487361980740625) * 1 / v^10 +# +# For large orders these formulas will converge much faster than using upward recurrence. +# If using up to the fifth U polynomial, it requires evaluation of a 15 degree polynomial. +# For tenth U polynomial, it requires evaluation of a 30 degree polynomial. +# +# # [1] "Rational Approximations for the Modified Bessel Function of the Second Kind # - K0(x) for Computations with Double Precision" by Pavel Holoborodko # [2] "Rational Approximations for the Modified Bessel Function of the Second Kind # - K1(x) for Computations with Double Precision" by Pavel Holoborodko # [3] https://github.com/boostorg/math/tree/develop/include/boost/math/special_functions/detail -# [4] "Computation of Bessel Functions of Complex Argument and Large ORder" by Donald E. Amos +# [4] "Computation of Bessel Functions of Complex Argument and Large Order" by Donald E. Amos # Sandia National Laboratories +# [5] https://dlmf.nist.gov/10.41 +# [6] Abramowitz, Milton, and Irene A. Stegun, eds. Handbook of mathematical functions with formulas, graphs, and mathematical tables. +# Vol. 55. US Government printing office, 1964. # """ besselk0(x::T) where T <: Union{Float32, Float64} @@ -130,14 +166,12 @@ function besselk1x(x::T) where T <: Union{Float32, Float64} end end #= -Recurrence is used for all values of nu. If besselk0(x) or besselk1(0) is equal to zero +If besselk0(x) or besselk1(0) is equal to zero this will underflow and always return zero even if besselk(x, nu) is larger than the smallest representing floating point value. In other words, for large values of x and small to moderate values of nu, this routine will underflow to zero. For small to medium values of x and large values of nu this will overflow and return Inf. - -In the future, a more efficient algorithm for large nu should be incorporated. =# """ besselk(x::T) where T <: Union{Float32, Float64} @@ -148,7 +182,7 @@ function besselk(nu::Int, x::T) where T <: Union{Float32, Float64} if nu < 100 return three_term_recurrence(x, besselk0(x), besselk1(x), nu) else - return k1(BigFloat(nu), BigFloat(x)) + return besselk_large_orders(BigFloat(nu), BigFloat(x)) end end """ @@ -179,7 +213,7 @@ end end -function k1(v, x::T) where T <: AbstractFloat +function besselk_large_orders(v, x::T) where T <: AbstractFloat z = x / v zs = sqrt(1 + z^2) From 292f29252a4ae30fc5aac08c6cd282bb240a62c9 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Mon, 17 Jan 2022 22:44:28 -0500 Subject: [PATCH 4/9] add more reference --- src/besselk.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/besselk.jl b/src/besselk.jl index d71ae3b..5280123 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -46,6 +46,7 @@ # The boundary should be carefully determined for accuracy and machine roundoff. # We use 10.41.4 from the Digital Library of Math Functions [5]. # This is also 9.7.8 in Abramowitz and Stegun [6]. +# K_{nu}(nu*z) = sqrt(pi / 2nu) *exp(-nu*n)/(1+z^2)^1/4 * sum((-1^k)U_k(p) /nu^k)) for k=0 -> infty # The U polynomials are the most tricky. They are listed up to order 4 in Table 9.39 # of [6]. For Float32, >=4 U polynomials are usually necessary. For Float64 values, # >= 8 orders are needed. However, this largely depends on the cutoff of order you need. @@ -54,7 +55,7 @@ # # However, calculation of these higher order U polynomials are tedious. These have been hand # calculated and somewhat crosschecked with symbolic math. There could be errors. They are listed -# here as a reference as higher orders are impossible to find and needed for any meaningfully accurate calculation. +# here as a reference as higher orders are impossible to find while being needed for any meaningfully accurate calculation. # u0 = one(x) # u1 = p / 24 * (3 - 5*p^2) * -1 / v From 89cd49dc596c95c14941db9f69b373e079ff468f Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Tue, 18 Jan 2022 00:45:08 -0500 Subject: [PATCH 5/9] add u11,u12,u13 --- src/besselk.jl | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index 5280123..0c2283c 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -75,7 +75,14 @@ # u10 = p^10 / 88580102706155225088000 * (290938676920391671935028890625*p^20 - 1745632061522350031610173343750*p^18 + 4513386761946134740461797128125*p^16 # - 6564241639632418015173104205000*p^14 + 5876803711285273203043452095250*p^12 - 3327704366990695147540934069220*p^10 + 1177120360439828012193658602930*p^8 # - 246750339886026017414509498824*p^6 + 27299183373230345667273718125*p^4 - 1230031256571145165088463750*p^2 + 9745329584487361980740625) * 1 / v^10 -# +# u11 = p^11 / 27636992044320430227456000 * (-1363312191078326248171914924296875*p^22 + 8997860461116953237934638500359375*p^20 - 25963913760458280822131901737878125*p^18 + 42936745153513007436411401865860625*p^16 - 44801790321820682384740638703320750*p^14 + 30589806122850866110941936529264950*p^12 - 13704902022868787041097596217578170*p^10 +# + 3926191452593448964331218647028194*p^8 - 676389476843440422173605288596087*p^6 + 62011003282542082472252466220875*p^4 - 2321657500166464779660536015625*p^2 + 15237265774872558064250728125) * -1 / v^11 +# u12 = p^12 / 39797268543821419527536640000 * (32426380464797989812768996474401171875*p^24 - 233469939346545526651936774615688437500*p^22 + 743739612850105971846901081692858843750*p^20 - 1378260730939829908037976894025628887500*p^18 + 1642838631056253395899341314188134178125*p^16 - 1314368459332124683504418467809129387000*p^14 +# + 714528665351965363868467348116538170900*p^12 - 261201165596865827608687437905740780920*p^10 + 62055079517573388459132793029571876461*p^8 - 8958590476947726766450604043798559500*p^6 + 692277766674325563109617997687563750*p^4 - 21882222767154197351962677311437500*p^2 + 120907703923613748239829527671875) * 1 / v^12 +# u13 = p^13 / 955134445051714068660879360000 * (-14020668045586884672121117629431460546875*p^26 + 109361210755577700442544717509565392265625*p^24 - 381190503845282445953419057314665534156250*p^22 + 782463969315283937856703223178540650343750*p^20 +# - 1049095945162229046324321461816274931715625*p^18 + 962926533925253374359704288384340809260875*p^16 - 616410216242554698436702353237166008656700*p^14 + 274983827478138958623041508409195988431140*p^12 - 83924867223075156862785921508524155665245*p^10 +# + 16843538631795795357786827345307534156063*p^8 - 2069933923586966756183324291232117362250*p^6 + 136735019134677724428035696765082813750*p^4 - 3698121486504259988094897605296209375*p^2 + 17438611142828905996129258798828125) * -1 / v^13 +# # For large orders these formulas will converge much faster than using upward recurrence. # If using up to the fifth U polynomial, it requires evaluation of a 15 degree polynomial. # For tenth U polynomial, it requires evaluation of a 30 degree polynomial. @@ -234,10 +241,19 @@ function besselk_large_orders(v, x::T) where T <: AbstractFloat u8 = p^8 / 22191183337881600 * (448357133137441653125*p^16 - 2152114239059719935000*p^14 + 4272845805510421639500*p^12 - 4513690624987320777000*p^10 + 2711772922412520971550*p^8 - 914113758588905038248*p^6 + 157768535329832893644*p^4 - 10960565081605263000*p^2 + 134790179652253125) * 1 / v^8 u9 = p^9 / 263631258054033408000 * (-64041091111686478524109375*p^18 + 345821892003106984030190625*p^16 - 790370708270219620781737500*p^14 + 992115946599792610768672500*p^12 - 741743213039573443221773250*p^10 + 334380732677827878090447630*p^8 - 87432034049652400520788332*p^6 + 11921080954211358275362500*p^4 - 659033454841709672064375*p^2 + 6427469716717690265625) * -1 / v^9 u10 = p^10 / 88580102706155225088000 * (290938676920391671935028890625*p^20 - 1745632061522350031610173343750*p^18 + 4513386761946134740461797128125*p^16 - 6564241639632418015173104205000*p^14 + 5876803711285273203043452095250*p^12 - 3327704366990695147540934069220*p^10 + 1177120360439828012193658602930*p^8 - 246750339886026017414509498824*p^6 + 27299183373230345667273718125*p^4 - 1230031256571145165088463750*p^2 + 9745329584487361980740625) * 1 / v^10 - out1 = ((u10 + u9) + u8) + u7 + u11 = p^11 / 27636992044320430227456000 * (-1363312191078326248171914924296875*p^22 + 8997860461116953237934638500359375*p^20 - 25963913760458280822131901737878125*p^18 + 42936745153513007436411401865860625*p^16 - 44801790321820682384740638703320750*p^14 + 30589806122850866110941936529264950*p^12 - 13704902022868787041097596217578170*p^10 + + 3926191452593448964331218647028194*p^8 - 676389476843440422173605288596087*p^6 + 62011003282542082472252466220875*p^4 - 2321657500166464779660536015625*p^2 + 15237265774872558064250728125) * -1 / v^11 + u12 = p^12 / 39797268543821419527536640000 * (32426380464797989812768996474401171875*p^24 - 233469939346545526651936774615688437500*p^22 + 743739612850105971846901081692858843750*p^20 - 1378260730939829908037976894025628887500*p^18 + 1642838631056253395899341314188134178125*p^16 - 1314368459332124683504418467809129387000*p^14 + + 714528665351965363868467348116538170900*p^12 - 261201165596865827608687437905740780920*p^10 + 62055079517573388459132793029571876461*p^8 - 8958590476947726766450604043798559500*p^6 + 692277766674325563109617997687563750*p^4 - 21882222767154197351962677311437500*p^2 + 120907703923613748239829527671875) * 1 / v^12 + u13 = p^13 / 955134445051714068660879360000 * (-14020668045586884672121117629431460546875*p^26 + 109361210755577700442544717509565392265625*p^24 - 381190503845282445953419057314665534156250*p^22 + 782463969315283937856703223178540650343750*p^20 + - 1049095945162229046324321461816274931715625*p^18 + 962926533925253374359704288384340809260875*p^16 - 616410216242554698436702353237166008656700*p^14 + 274983827478138958623041508409195988431140*p^12 - 83924867223075156862785921508524155665245*p^10 + + 16843538631795795357786827345307534156063*p^8 - 2069933923586966756183324291232117362250*p^6 + 136735019134677724428035696765082813750*p^4 - 3698121486504259988094897605296209375*p^2 + 17438611142828905996129258798828125) * -1 / v^13 + + + + out1 = ((u13 + u12 + u11 + u10 + u9) + u8) + u7 out2 = ((out1 + u6 + u5) + u4) + u3 out = ((out2 + u2) + u1) + u0 return coef*out end - \ No newline at end of file From 2ff2cafeee0eb62d369a28af6bc277dfd494e695 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Tue, 18 Jan 2022 17:05:24 -0500 Subject: [PATCH 6/9] add float32 versions --- src/besselk.jl | 110 ++++++++++++++++++++++++++++++++++++++++--- test/besselk_test.jl | 12 +++++ 2 files changed, 115 insertions(+), 7 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index 0c2283c..b12ae70 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -84,8 +84,6 @@ # + 16843538631795795357786827345307534156063*p^8 - 2069933923586966756183324291232117362250*p^6 + 136735019134677724428035696765082813750*p^4 - 3698121486504259988094897605296209375*p^2 + 17438611142828905996129258798828125) * -1 / v^13 # # For large orders these formulas will converge much faster than using upward recurrence. -# If using up to the fifth U polynomial, it requires evaluation of a 15 degree polynomial. -# For tenth U polynomial, it requires evaluation of a 30 degree polynomial. # # # [1] "Rational Approximations for the Modified Bessel Function of the Second Kind @@ -186,13 +184,20 @@ For small to medium values of x and large values of nu this will overflow and re Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``. """ -function besselk(nu::Int, x::T) where T <: Union{Float32, Float64} - if nu < 100 +function besselk(nu::Int, x::T) where T <: Float64 + if nu < 50 return three_term_recurrence(x, besselk0(x), besselk1(x), nu) else return besselk_large_orders(BigFloat(nu), BigFloat(x)) end end +function besselk(nu::Int, x::T) where T <: Float32 + if nu < 20 + return three_term_recurrence(x, besselk0(x), besselk1(x), nu) + else + return k4(nu, x) + end +end """ besselk(x::T) where T <: Union{Float32, Float64} @@ -220,7 +225,71 @@ end return k2 end +function k4(v, x::T) where T <: Float32 + x = Float64(x) + z = x / v + zs = sqrt(1 + z^2) + n = zs + log(z / (1 + zs)) + coef = exp(-v * n) * sqrt(pi / (2*v)) / sqrt(zs) + p = inv(zs) + p2 = p*p + u0 = one(Float64) + u1 = 1 / 24 * evalpoly(p2, (3, -5)) + u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385)) + u3 = 1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425)) + u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725)) + return Float32(coef*evalpoly(-p/v, (u0, u1, u2, u3, u4))) + end + +function besselk_large_orders(v, x::T) where T <: AbstractFloat + x = Float64(x) + z = x / v + zs = sqrt(1 + z^2) + n = zs + log(z / (1 + zs)) + coef = exp(-v * n) * sqrt(pi / (2*v)) / sqrt(zs) + p = inv(zs) + p2 = p*p + u0 = one(Float64) + u1 = 1 / 24 * evalpoly(p2, (3, -5)) + u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385)) + u3 = 1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425)) + u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725)) + u5 = 1 / 6688604160 * evalpoly(p2, (1519035525, -49286948607, 284499769554, -614135872350, 566098157625, -188699385875)) + u6 = 1 / 4815794995200 * evalpoly(p2, (2757049477875, -127577298354750, 1050760774457901, -3369032068261860,5104696716244125, -3685299006138750, 1023694168371875)) + u7 = 1 / 115579079884800 * evalpoly(p2, (199689155040375, -12493049053044375, 138799253740521843, -613221795981706275, 1347119637570231525, -1570320948552481125, 931766432052080625, -221849150488590625)) + u8 = 1 / 22191183337881600 * evalpoly(p2, (134790179652253125, -10960565081605263000, 157768535329832893644, -914113758588905038248, 2711772922412520971550, -4513690624987320777000, 4272845805510421639500, -2152114239059719935000, 448357133137441653125)) + return coef*evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8)) +end + +#= +function besselk_large_orders(v, x::T) where T <: AbstractFloat + x = Float64(x) + z = x / v + zs = sqrt(1 + z^2) + n = zs + log(z / (1 + zs)) + coef = exp(-v * n) * sqrt(pi / (2*v)) / sqrt(zs) + p = inv(zs) + p2 = p*p + u0 = one(Float64) + u1 = 1 / 24 * evalpoly(p2, (3, -5)) + u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385)) + u3 = 1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425)) + u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725)) + u5 = 1 / 6688604160 * evalpoly(p2, (1519035525, -49286948607, 284499769554, -614135872350, 566098157625, -188699385875)) + u6 = 1 / 4815794995200 * evalpoly(p2, (2757049477875, -127577298354750, 1050760774457901, -3369032068261860,5104696716244125, -3685299006138750, 1023694168371875)) + u7 = 1 / 115579079884800 * evalpoly(p2, (199689155040375, -12493049053044375, 138799253740521843, -613221795981706275, 1347119637570231525, -1570320948552481125, 931766432052080625, -221849150488590625)) + u8 = 1 / 22191183337881600 * evalpoly(p2, (134790179652253125, -10960565081605263000, 157768535329832893644, -914113758588905038248, 2711772922412520971550, -4513690624987320777000, 4272845805510421639500, -2152114239059719935000, 448357133137441653125)) + u9 = 1 / 263631258054033408000 * evalpoly(p2, (6427469716717690265625, -659033454841709672064375, 11921080954211358275362500, -87432034049652400520788332, 334380732677827878090447630, -741743213039573443221773250, 992115946599792610768672500, -790370708270219620781737500, 345821892003106984030190625, -64041091111686478524109375)) + u10 = 1 / 88580102706155225088000 * evalpoly(p2, (9745329584487361980740625, -1230031256571145165088463750, 27299183373230345667273718125, -246750339886026017414509498824, 1177120360439828012193658602930, -3327704366990695147540934069220, 5876803711285273203043452095250, -6564241639632418015173104205000, 4513386761946134740461797128125, -1745632061522350031610173343750, 290938676920391671935028890625)) + + return coef*evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10)) +end + +=# + + +#= function besselk_large_orders(v, x::T) where T <: AbstractFloat z = x / v zs = sqrt(1 + z^2) @@ -248,12 +317,39 @@ function besselk_large_orders(v, x::T) where T <: AbstractFloat u13 = p^13 / 955134445051714068660879360000 * (-14020668045586884672121117629431460546875*p^26 + 109361210755577700442544717509565392265625*p^24 - 381190503845282445953419057314665534156250*p^22 + 782463969315283937856703223178540650343750*p^20 - 1049095945162229046324321461816274931715625*p^18 + 962926533925253374359704288384340809260875*p^16 - 616410216242554698436702353237166008656700*p^14 + 274983827478138958623041508409195988431140*p^12 - 83924867223075156862785921508524155665245*p^10 + 16843538631795795357786827345307534156063*p^8 - 2069933923586966756183324291232117362250*p^6 + 136735019134677724428035696765082813750*p^4 - 3698121486504259988094897605296209375*p^2 + 17438611142828905996129258798828125) * -1 / v^13 - - + u14 = p^14 / 45846453362482275295722209280000 * (13133360053559028970728309756597440972265625*p^28 - 110320224449895843354117801955418504167031250*p^26 + 417630985812040040477569028872405152769921875*p^24 + - 940627071986145750405920450097257805227812500*p^22 + 1401302601668131482630653233972052729323190625*p^20 - 1451823699927947602004297385351260623500372750*p^18 + 1070439683260179398514645952824098811025619475*p^16 + - 564850830044980230366140582063618983657685400*p^14 + 211477117385619365164298957821904115470535555*p^12 - 54857705817658080981995319669299096598482382*p^10 + 9440449669103391509091075981237243128469201*p^8 + - 1000503839668383458999731491914516564625300*p^6 + 57170953417612444837142230812990944671875*p^4 - 1338184074771428116079233795614103631250*p^2 + 5448320367052402487647812713291015625) * 1 / v^14 + u15 = p^15 / 9820310310243703368343697227776000000 * (- 59115551939078562227317999668652486368337724609375*p^30 + 532039967451707060045861997017872377315039521484375*p^28 - 2173722139119126857509156976742601742357422740234375*p^26 + + 5329871927856528282113994744458999865006055974609375*p^24 - 8735135969643867540297524795790262235822823374296875*p^22 + 10085018700249896522602267572484630409640764997271875*p^20 - 8420533422834140468835467666391400380550043688871875*p^18 + + 5136561256208409671660362298619778869859994724706875*p^16 - 2284251621937242886581917667066122422330060024456125*p^14 + 730367145705123976114617970888707594104468177381925*p^12 - 163359140754958502896104062604202934925448173291477*p^10 + + 24403480234538299231733883413666768614198435948125*p^8 - 2254933495791765108580529087615802250458013685625*p^6 + 112597271053778779753048514469995937998172890625*p^4 - 2303431987527333128955769182299845911305390625*p^2 + 8178936810213560828419581728001773291015625) * -1 / v^15 + - out1 = ((u13 + u12 + u11 + u10 + u9) + u8) + u7 + out1 = ((u15 + u14 + u13 + u12 + u11 + u10 + u9) + u8) + u7 out2 = ((out1 + u6 + u5) + u4) + u3 out = ((out2 + u2) + u1) + u0 return coef*out end + + + + + u11 = 1 / 27636992044320430227456000 * (-1363312191078326248171914924296875*p^22 + 8997860461116953237934638500359375*p^20 - 25963913760458280822131901737878125*p^18 + 42936745153513007436411401865860625*p^16 - 44801790321820682384740638703320750*p^14 + 30589806122850866110941936529264950*p^12 - 13704902022868787041097596217578170*p^10 + + 3926191452593448964331218647028194*p^8 - 676389476843440422173605288596087*p^6 + 62011003282542082472252466220875*p^4 - 2321657500166464779660536015625*p^2 + 15237265774872558064250728125) * -1 / v^11 + u12 = 1 / 39797268543821419527536640000 * (32426380464797989812768996474401171875*p^24 - 233469939346545526651936774615688437500*p^22 + 743739612850105971846901081692858843750*p^20 - 1378260730939829908037976894025628887500*p^18 + 1642838631056253395899341314188134178125*p^16 - 1314368459332124683504418467809129387000*p^14 + + 714528665351965363868467348116538170900*p^12 - 261201165596865827608687437905740780920*p^10 + 62055079517573388459132793029571876461*p^8 - 8958590476947726766450604043798559500*p^6 + 692277766674325563109617997687563750*p^4 - 21882222767154197351962677311437500*p^2 + 120907703923613748239829527671875) * 1 / v^12 + u13 = 1 / 955134445051714068660879360000 * (-14020668045586884672121117629431460546875*p^26 + 109361210755577700442544717509565392265625*p^24 - 381190503845282445953419057314665534156250*p^22 + 782463969315283937856703223178540650343750*p^20 + - 1049095945162229046324321461816274931715625*p^18 + 962926533925253374359704288384340809260875*p^16 - 616410216242554698436702353237166008656700*p^14 + 274983827478138958623041508409195988431140*p^12 - 83924867223075156862785921508524155665245*p^10 + + 16843538631795795357786827345307534156063*p^8 - 2069933923586966756183324291232117362250*p^6 + 136735019134677724428035696765082813750*p^4 - 3698121486504259988094897605296209375*p^2 + 17438611142828905996129258798828125) * -1 / v^13 + u14 = 1 / 45846453362482275295722209280000 * (13133360053559028970728309756597440972265625*p^28 - 110320224449895843354117801955418504167031250*p^26 + 417630985812040040477569028872405152769921875*p^24 + - 940627071986145750405920450097257805227812500*p^22 + 1401302601668131482630653233972052729323190625*p^20 - 1451823699927947602004297385351260623500372750*p^18 + 1070439683260179398514645952824098811025619475*p^16 + - 564850830044980230366140582063618983657685400*p^14 + 211477117385619365164298957821904115470535555*p^12 - 54857705817658080981995319669299096598482382*p^10 + 9440449669103391509091075981237243128469201*p^8 + - 1000503839668383458999731491914516564625300*p^6 + 57170953417612444837142230812990944671875*p^4 - 1338184074771428116079233795614103631250*p^2 + 5448320367052402487647812713291015625) * 1 / v^14 + u15 = 1 / 9820310310243703368343697227776000000 * (- 59115551939078562227317999668652486368337724609375*p^30 + 532039967451707060045861997017872377315039521484375*p^28 - 2173722139119126857509156976742601742357422740234375*p^26 + + 5329871927856528282113994744458999865006055974609375*p^24 - 8735135969643867540297524795790262235822823374296875*p^22 + 10085018700249896522602267572484630409640764997271875*p^20 - 8420533422834140468835467666391400380550043688871875*p^18 + + 5136561256208409671660362298619778869859994724706875*p^16 - 2284251621937242886581917667066122422330060024456125*p^14 + 730367145705123976114617970888707594104468177381925*p^12 - 163359140754958502896104062604202934925448173291477*p^10 + + 24403480234538299231733883413666768614198435948125*p^8 - 2254933495791765108580529087615802250458013685625*p^6 + 112597271053778779753048514469995937998172890625*p^4 - 2303431987527333128955769182299845911305390625*p^2 + 8178936810213560828419581728001773291015625) * -1 / v^15 +=# \ No newline at end of file diff --git a/test/besselk_test.jl b/test/besselk_test.jl index 673307a..1638cf8 100644 --- a/test/besselk_test.jl +++ b/test/besselk_test.jl @@ -81,6 +81,18 @@ m = 200:5:1000; x = 400.0:10.0:1200.0 t = Float64.([besselk(m, x) for m in m, x in x]) @test t ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x] +# Float 32 tests for aysmptotic expansion +m = 20:5:200; x = 5.0f0:2.0f0:400.0f0 +t = [besselk(m, x) for m in m, x in x] +@test t[10] isa Float32 +@test t ≈ Float32.([SpecialFunctions.besselk(m, x) for m in m, x in x]) + +# test for low values and medium orders +m = 20:5:50; x = [1f-3, 1f-2, 1f-1, 1f0, 1.5f0, 2.0f0, 4.0f0] +t = [besselk(m, x) for m in m, x in x] +@test t[5] isa Float32 +@test t ≈ Float32.([SpecialFunctions.besselk(m, x) for m in m, x in x]) + @test iszero(besselk(20, 1000.0)) #@test isinf(besselk(250, 5.0)) From b5227b8b3c3234e02c603e6b0d377bd9965641d6 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 19 Jan 2022 19:09:20 -0500 Subject: [PATCH 7/9] add besseli and besselk for any real nu --- src/Bessels.jl | 3 + src/U_polynomials.jl | 135 ++++++++++++++++++++++++++++++++++ src/besseli.jl | 81 +++++++++++++++++++++ src/besselk.jl | 163 ++++++++---------------------------------- src/math_constants.jl | 5 ++ test/besseli_test.jl | 18 +++++ 6 files changed, 272 insertions(+), 133 deletions(-) create mode 100644 src/U_polynomials.jl diff --git a/src/Bessels.jl b/src/Bessels.jl index bf1885d..6146595 100644 --- a/src/Bessels.jl +++ b/src/Bessels.jl @@ -7,6 +7,8 @@ export besselj export bessely0 export bessely1 +export besseli +export besselix export besseli0 export besseli0x export besseli1 @@ -32,6 +34,7 @@ include("Float128/bessely.jl") include("Float128/constants.jl") include("math_constants.jl") +include("U_polynomials.jl") #include("parse.jl") #include("hankel.jl") diff --git a/src/U_polynomials.jl b/src/U_polynomials.jl new file mode 100644 index 0000000..a9d3c7c --- /dev/null +++ b/src/U_polynomials.jl @@ -0,0 +1,135 @@ +function Uk_poly_Kn(p, v, p2, ::Type{Float32}) + u0 = one(p) + u1 = 1 / 24 * evalpoly(p2, (3, -5)) + u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385)) + u3 = 1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425)) + u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725)) + return evalpoly(-p/v, (u0, u1, u2, u3, u4)) +end +function Uk_poly_Kn(p, v, p2, ::Type{T}) where T <: Float64 + u0 = one(T) + u1 = 1 / 24 * evalpoly(p2, (3, -5)) + u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385)) + u3 = 1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425)) + u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725)) + u5 = 1 / 6688604160 * evalpoly(p2, (1519035525, -49286948607, 284499769554, -614135872350, 566098157625, -188699385875)) + u6 = 1 / 4815794995200 * evalpoly(p2, (2757049477875, -127577298354750, 1050760774457901, -3369032068261860,5104696716244125, -3685299006138750, 1023694168371875)) + return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6)) +end +function Uk_poly_In(p, v, p2, ::Type{T}) where T <: Float64 + u0 = one(T) + u1 = -1 / 24 * evalpoly(p2, (3, -5)) + u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385)) + u3 = -1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425)) + u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725)) + u5 = -1 / 6688604160 * evalpoly(p2, (1519035525, -49286948607, 284499769554, -614135872350, 566098157625, -188699385875)) + u6 = 1 / 4815794995200 * evalpoly(p2, (2757049477875, -127577298354750, 1050760774457901, -3369032068261860,5104696716244125, -3685299006138750, 1023694168371875)) + return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6)) +end +function Uk_poly_In(p, v, p2, ::Type{Float32}) + u0 = one(p) + u1 = -1 / 24 * evalpoly(p2, (3, -5)) + u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385)) + u3 = -1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425)) + u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725)) + return evalpoly(-p/v, (u0, u1, u2, u3, u4)) + end + +#= +function Uk_poly_Kn(p, v, p2, ::Type{T}) where T <: BigFloat + u0 = one(T) + u1 = 1 / 24 * evalpoly(p2, (3, -5)) + u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385)) + u3 = 1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425)) + u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725)) + u5 = 1 / 6688604160 * evalpoly(p2, (1519035525, -49286948607, 284499769554, -614135872350, 566098157625, -188699385875)) + u6 = 1 / 4815794995200 * evalpoly(p2, (2757049477875, -127577298354750, 1050760774457901, -3369032068261860,5104696716244125, -3685299006138750, 1023694168371875)) + u7 = 1 / 115579079884800 * evalpoly(p2, (199689155040375, -12493049053044375, 138799253740521843, -613221795981706275, 1347119637570231525, -1570320948552481125, 931766432052080625, -221849150488590625)) + u8 = 1 / 22191183337881600 * evalpoly(p2, (134790179652253125, -10960565081605263000, 157768535329832893644, -914113758588905038248, 2711772922412520971550, -4513690624987320777000, 4272845805510421639500, -2152114239059719935000, 448357133137441653125)) + u9 = 1 / 263631258054033408000 * evalpoly(p2, (6427469716717690265625, -659033454841709672064375, 11921080954211358275362500, -87432034049652400520788332, 334380732677827878090447630, -741743213039573443221773250, 992115946599792610768672500, -790370708270219620781737500, 345821892003106984030190625, -64041091111686478524109375)) + u10 = 1 / 88580102706155225088000 * evalpoly(p2, (9745329584487361980740625, -1230031256571145165088463750, 27299183373230345667273718125, -246750339886026017414509498824, 1177120360439828012193658602930, -3327704366990695147540934069220, 5876803711285273203043452095250, -6564241639632418015173104205000, 4513386761946134740461797128125, -1745632061522350031610173343750, 290938676920391671935028890625)) + return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10)) +end +function besselk_large_orders(v, x::T) where T <: AbstractFloat + x = Float64(x) + z = x / v + zs = sqrt(1 + z^2) + n = zs + log(z / (1 + zs)) + coef = exp(-v * n) * sqrt(pi / (2*v)) / sqrt(zs) + p = inv(zs) + p2 = p*p + u0 = one(Float64) + u1 = 1 / 24 * evalpoly(p2, (3, -5)) + u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385)) + u3 = 1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425)) + u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725)) + u5 = 1 / 6688604160 * evalpoly(p2, (1519035525, -49286948607, 284499769554, -614135872350, 566098157625, -188699385875)) + u6 = 1 / 4815794995200 * evalpoly(p2, (2757049477875, -127577298354750, 1050760774457901, -3369032068261860,5104696716244125, -3685299006138750, 1023694168371875)) + u7 = 1 / 115579079884800 * evalpoly(p2, (199689155040375, -12493049053044375, 138799253740521843, -613221795981706275, 1347119637570231525, -1570320948552481125, 931766432052080625, -221849150488590625)) + u8 = 1 / 22191183337881600 * evalpoly(p2, (134790179652253125, -10960565081605263000, 157768535329832893644, -914113758588905038248, 2711772922412520971550, -4513690624987320777000, 4272845805510421639500, -2152114239059719935000, 448357133137441653125)) + u9 = 1 / 263631258054033408000 * evalpoly(p2, (6427469716717690265625, -659033454841709672064375, 11921080954211358275362500, -87432034049652400520788332, 334380732677827878090447630, -741743213039573443221773250, 992115946599792610768672500, -790370708270219620781737500, 345821892003106984030190625, -64041091111686478524109375)) + u10 = 1 / 88580102706155225088000 * evalpoly(p2, (9745329584487361980740625, -1230031256571145165088463750, 27299183373230345667273718125, -246750339886026017414509498824, 1177120360439828012193658602930, -3327704366990695147540934069220, 5876803711285273203043452095250, -6564241639632418015173104205000, 4513386761946134740461797128125, -1745632061522350031610173343750, 290938676920391671935028890625)) + + return coef*evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10)) +end + +function besselk_large_orders(v, x::T) where T <: AbstractFloat + z = x / v + zs = sqrt(1 + z^2) + + n = zs + log(z / (1 + zs)) + coef = exp(-v * n) * sqrt(pi / (2*v)) / sqrt(zs) + + p = inv(zs) + + u0 = one(x) + u1 = p / 24 * (3 - 5*p^2) * -1 / v + u2 = p^2 / 1152 * (81 - 462*p^2 + 385*p^4) / v^2 + u3 = p^3 / 414720 * (30375 - 369603 * p^2 + 765765*p^4 - 425425*p^6) * -1 / v^3 + u4 = p^4 / 39813120 * (4465125 - 94121676*p^2 + 349922430*p^4 - 446185740*p^6 + 185910725*p^8) / v^4 + u5 = p^5 / 6688604160 * (-188699385875*p^10 + 566098157625*p^8 - 614135872350*p^6 + 284499769554*p^4 - 49286948607*p^2 + 1519035525) * -1 / v^5 + u6 = p^6 / 4815794995200 * (1023694168371875*p^12 - 3685299006138750*p^10 + 5104696716244125*p^8 - 3369032068261860*p^6 + 1050760774457901*p^4 - 127577298354750*p^2 + 2757049477875) * 1 / v^6 + u7 = p^7 / 115579079884800 * (-221849150488590625*p^14 + 931766432052080625*p^12 - 1570320948552481125*p^10 + 1347119637570231525*p^8 - 613221795981706275*p^6 + 138799253740521843*p^4 - 12493049053044375*p^2 + 199689155040375) * -1 / v^7 + u8 = p^8 / 22191183337881600 * (448357133137441653125*p^16 - 2152114239059719935000*p^14 + 4272845805510421639500*p^12 - 4513690624987320777000*p^10 + 2711772922412520971550*p^8 - 914113758588905038248*p^6 + 157768535329832893644*p^4 - 10960565081605263000*p^2 + 134790179652253125) * 1 / v^8 + u9 = p^9 / 263631258054033408000 * (-64041091111686478524109375*p^18 + 345821892003106984030190625*p^16 - 790370708270219620781737500*p^14 + 992115946599792610768672500*p^12 - 741743213039573443221773250*p^10 + 334380732677827878090447630*p^8 - 87432034049652400520788332*p^6 + 11921080954211358275362500*p^4 - 659033454841709672064375*p^2 + 6427469716717690265625) * -1 / v^9 + u10 = p^10 / 88580102706155225088000 * (290938676920391671935028890625*p^20 - 1745632061522350031610173343750*p^18 + 4513386761946134740461797128125*p^16 - 6564241639632418015173104205000*p^14 + 5876803711285273203043452095250*p^12 - 3327704366990695147540934069220*p^10 + 1177120360439828012193658602930*p^8 - 246750339886026017414509498824*p^6 + 27299183373230345667273718125*p^4 - 1230031256571145165088463750*p^2 + 9745329584487361980740625) * 1 / v^10 + u11 = p^11 / 27636992044320430227456000 * (-1363312191078326248171914924296875*p^22 + 8997860461116953237934638500359375*p^20 - 25963913760458280822131901737878125*p^18 + 42936745153513007436411401865860625*p^16 - 44801790321820682384740638703320750*p^14 + 30589806122850866110941936529264950*p^12 - 13704902022868787041097596217578170*p^10 + + 3926191452593448964331218647028194*p^8 - 676389476843440422173605288596087*p^6 + 62011003282542082472252466220875*p^4 - 2321657500166464779660536015625*p^2 + 15237265774872558064250728125) * -1 / v^11 + u12 = p^12 / 39797268543821419527536640000 * (32426380464797989812768996474401171875*p^24 - 233469939346545526651936774615688437500*p^22 + 743739612850105971846901081692858843750*p^20 - 1378260730939829908037976894025628887500*p^18 + 1642838631056253395899341314188134178125*p^16 - 1314368459332124683504418467809129387000*p^14 + + 714528665351965363868467348116538170900*p^12 - 261201165596865827608687437905740780920*p^10 + 62055079517573388459132793029571876461*p^8 - 8958590476947726766450604043798559500*p^6 + 692277766674325563109617997687563750*p^4 - 21882222767154197351962677311437500*p^2 + 120907703923613748239829527671875) * 1 / v^12 + u13 = p^13 / 955134445051714068660879360000 * (-14020668045586884672121117629431460546875*p^26 + 109361210755577700442544717509565392265625*p^24 - 381190503845282445953419057314665534156250*p^22 + 782463969315283937856703223178540650343750*p^20 + - 1049095945162229046324321461816274931715625*p^18 + 962926533925253374359704288384340809260875*p^16 - 616410216242554698436702353237166008656700*p^14 + 274983827478138958623041508409195988431140*p^12 - 83924867223075156862785921508524155665245*p^10 + + 16843538631795795357786827345307534156063*p^8 - 2069933923586966756183324291232117362250*p^6 + 136735019134677724428035696765082813750*p^4 - 3698121486504259988094897605296209375*p^2 + 17438611142828905996129258798828125) * -1 / v^13 + u14 = p^14 / 45846453362482275295722209280000 * (13133360053559028970728309756597440972265625*p^28 - 110320224449895843354117801955418504167031250*p^26 + 417630985812040040477569028872405152769921875*p^24 + - 940627071986145750405920450097257805227812500*p^22 + 1401302601668131482630653233972052729323190625*p^20 - 1451823699927947602004297385351260623500372750*p^18 + 1070439683260179398514645952824098811025619475*p^16 + - 564850830044980230366140582063618983657685400*p^14 + 211477117385619365164298957821904115470535555*p^12 - 54857705817658080981995319669299096598482382*p^10 + 9440449669103391509091075981237243128469201*p^8 + - 1000503839668383458999731491914516564625300*p^6 + 57170953417612444837142230812990944671875*p^4 - 1338184074771428116079233795614103631250*p^2 + 5448320367052402487647812713291015625) * 1 / v^14 + u15 = p^15 / 9820310310243703368343697227776000000 * (- 59115551939078562227317999668652486368337724609375*p^30 + 532039967451707060045861997017872377315039521484375*p^28 - 2173722139119126857509156976742601742357422740234375*p^26 + + 5329871927856528282113994744458999865006055974609375*p^24 - 8735135969643867540297524795790262235822823374296875*p^22 + 10085018700249896522602267572484630409640764997271875*p^20 - 8420533422834140468835467666391400380550043688871875*p^18 + + 5136561256208409671660362298619778869859994724706875*p^16 - 2284251621937242886581917667066122422330060024456125*p^14 + 730367145705123976114617970888707594104468177381925*p^12 - 163359140754958502896104062604202934925448173291477*p^10 + + 24403480234538299231733883413666768614198435948125*p^8 - 2254933495791765108580529087615802250458013685625*p^6 + 112597271053778779753048514469995937998172890625*p^4 - 2303431987527333128955769182299845911305390625*p^2 + 8178936810213560828419581728001773291015625) * -1 / v^15 + + + out1 = ((u15 + u14 + u13 + u12 + u11 + u10 + u9) + u8) + u7 + out2 = ((out1 + u6 + u5) + u4) + u3 + out = ((out2 + u2) + u1) + u0 + + return coef*out + end + + u11 = 1 / 27636992044320430227456000 * (-1363312191078326248171914924296875*p^22 + 8997860461116953237934638500359375*p^20 - 25963913760458280822131901737878125*p^18 + 42936745153513007436411401865860625*p^16 - 44801790321820682384740638703320750*p^14 + 30589806122850866110941936529264950*p^12 - 13704902022868787041097596217578170*p^10 + + 3926191452593448964331218647028194*p^8 - 676389476843440422173605288596087*p^6 + 62011003282542082472252466220875*p^4 - 2321657500166464779660536015625*p^2 + 15237265774872558064250728125) * -1 / v^11 + u12 = 1 / 39797268543821419527536640000 * (32426380464797989812768996474401171875*p^24 - 233469939346545526651936774615688437500*p^22 + 743739612850105971846901081692858843750*p^20 - 1378260730939829908037976894025628887500*p^18 + 1642838631056253395899341314188134178125*p^16 - 1314368459332124683504418467809129387000*p^14 + + 714528665351965363868467348116538170900*p^12 - 261201165596865827608687437905740780920*p^10 + 62055079517573388459132793029571876461*p^8 - 8958590476947726766450604043798559500*p^6 + 692277766674325563109617997687563750*p^4 - 21882222767154197351962677311437500*p^2 + 120907703923613748239829527671875) * 1 / v^12 + u13 = 1 / 955134445051714068660879360000 * (-14020668045586884672121117629431460546875*p^26 + 109361210755577700442544717509565392265625*p^24 - 381190503845282445953419057314665534156250*p^22 + 782463969315283937856703223178540650343750*p^20 + - 1049095945162229046324321461816274931715625*p^18 + 962926533925253374359704288384340809260875*p^16 - 616410216242554698436702353237166008656700*p^14 + 274983827478138958623041508409195988431140*p^12 - 83924867223075156862785921508524155665245*p^10 + + 16843538631795795357786827345307534156063*p^8 - 2069933923586966756183324291232117362250*p^6 + 136735019134677724428035696765082813750*p^4 - 3698121486504259988094897605296209375*p^2 + 17438611142828905996129258798828125) * -1 / v^13 + u14 = 1 / 45846453362482275295722209280000 * (13133360053559028970728309756597440972265625*p^28 - 110320224449895843354117801955418504167031250*p^26 + 417630985812040040477569028872405152769921875*p^24 + - 940627071986145750405920450097257805227812500*p^22 + 1401302601668131482630653233972052729323190625*p^20 - 1451823699927947602004297385351260623500372750*p^18 + 1070439683260179398514645952824098811025619475*p^16 + - 564850830044980230366140582063618983657685400*p^14 + 211477117385619365164298957821904115470535555*p^12 - 54857705817658080981995319669299096598482382*p^10 + 9440449669103391509091075981237243128469201*p^8 + - 1000503839668383458999731491914516564625300*p^6 + 57170953417612444837142230812990944671875*p^4 - 1338184074771428116079233795614103631250*p^2 + 5448320367052402487647812713291015625) * 1 / v^14 + u15 = 1 / 9820310310243703368343697227776000000 * (- 59115551939078562227317999668652486368337724609375*p^30 + 532039967451707060045861997017872377315039521484375*p^28 - 2173722139119126857509156976742601742357422740234375*p^26 + + 5329871927856528282113994744458999865006055974609375*p^24 - 8735135969643867540297524795790262235822823374296875*p^22 + 10085018700249896522602267572484630409640764997271875*p^20 - 8420533422834140468835467666391400380550043688871875*p^18 + + 5136561256208409671660362298619778869859994724706875*p^16 - 2284251621937242886581917667066122422330060024456125*p^14 + 730367145705123976114617970888707594104468177381925*p^12 - 163359140754958502896104062604202934925448173291477*p^10 + + 24403480234538299231733883413666768614198435948125*p^8 - 2254933495791765108580529087615802250458013685625*p^6 + 112597271053778779753048514469995937998172890625*p^4 - 2303431987527333128955769182299845911305390625*p^2 + 8178936810213560828419581728001773291015625) * -1 / v^15 +=# diff --git a/src/besseli.jl b/src/besseli.jl index f4cbc7b..fe787a0 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -116,3 +116,84 @@ function besseli1x(x::T) where T <: Union{Float32, Float64} end return z end + +""" + besseli(nu, x::T) where T <: Union{Float32, Float64} + +Modified Bessel function of the first kind of order nu, ``I_{nu}(x)``. +Nu must be real. +""" +function besseli(nu, x::T) where T <: Union{Float32, Float64, BigFloat} + nu == 0 && return besseli0(x) + nu == 1 && return besseli1(x) + + branch = 30 + if nu < branch + inp1 = besseli_large_orders(branch + 1, x) + in = besseli_large_orders(branch, x) + return down_recurrence(x, in, inp1, nu, branch) + else + return besseli_large_orders(nu, x) + end +end +""" + besselix(nu, x::T) where T <: Union{Float32, Float64} + +Scaled modified Bessel function of the first kind of order nu, ``I_{nu}(x)*e^{-x}``. +Nu must be real. +""" +function besselix(nu, x::T) where T <: Union{Float32, Float64, BigFloat} + nu == 0 && return besseli0x(x) + nu == 1 && return besseli1x(x) + + branch = 30 + if nu < branch + inp1 = besseli_large_orders_scaled(branch + 1, x) + in = besseli_large_orders_scaled(branch, x) + return down_recurrence(x, in, inp1, nu, branch) + else + return besseli_large_orders_scaled(nu, x) + end +end + + +@inline function down_recurrence(x, in, inp1, nu, branch) + # this prevents us from looping through large values of nu when the loop will always return zero + (iszero(in) || iszero(inp1)) && return zero(x) + (isinf(inp1) || isinf(inp1)) && return in + + inm1 = in + x2 = 2 / x + for n in branch:-1:nu+1 + a = x2 * n + inm1 = muladd(a, in, inp1) + inp1 = in + in = inm1 + end + return inm1 +end + +function besseli_large_orders(v, x::T) where T <: Union{Float32, Float64, BigFloat} + S = promote_type(T, Float64) + x = S(x) + z = x / v + zs = hypot(1, z) + n = zs + log(z) - log1p(zs) + coef = SQ1O2PI(S) * sqrt(inv(v)) * exp(v*n) / sqrt(zs) + p = inv(zs) + p2 = v^2/fma(max(v,x), max(v,x), min(v,x)^2) + + return T(coef*Uk_poly_In(p, v, p2, T)) +end +function besseli_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64, BigFloat} + S = promote_type(T, Float64) + x = S(x) + z = x / v + zs = hypot(1, z) + n = zs + log(z) - log1p(zs) + coef = SQ1O2PI(S) * sqrt(inv(v)) * exp(v*n - x) / sqrt(zs) + p = inv(zs) + p2 = v^2/fma(max(v,x), max(v,x), min(v,x)^2) + + return T(coef*Uk_poly_In(p, v, p2, T)) +end \ No newline at end of file diff --git a/src/besselk.jl b/src/besselk.jl index b12ae70..82abe3b 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -184,30 +184,30 @@ For small to medium values of x and large values of nu this will overflow and re Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``. """ -function besselk(nu::Int, x::T) where T <: Float64 - if nu < 50 - return three_term_recurrence(x, besselk0(x), besselk1(x), nu) +function besselk(nu, x::T) where T <: Union{Float32, Float64, BigFloat} + T == Float32 ? branch = 20 : branch = 50 + if nu < branch + return up_recurrence(x, besselk0(x), besselk1(x), nu) else - return besselk_large_orders(BigFloat(nu), BigFloat(x)) - end -end -function besselk(nu::Int, x::T) where T <: Float32 - if nu < 20 - return three_term_recurrence(x, besselk0(x), besselk1(x), nu) - else - return k4(nu, x) + return besselk_large_orders(nu, x) end end + """ besselk(x::T) where T <: Union{Float32, Float64} Scaled modified Bessel function of the second kind of order nu, ``K_{nu}(x)*e^{x}``. """ function besselkx(nu::Int, x::T) where T <: Union{Float32, Float64} - return three_term_recurrence(x, besselk0x(x), besselk1x(x), nu) + T == Float32 ? branch = 20 : branch = 50 + if nu < branch + return up_recurrence(x, besselk0x(x), besselk1x(x), nu) + else + return besselk_large_orders_scaled(nu, x) + end end -@inline function three_term_recurrence(x, k0, k1, nu) +@inline function up_recurrence(x, k0, k1, nu) nu == 0 && return k0 nu == 1 && return k1 @@ -225,131 +225,28 @@ end return k2 end -function k4(v, x::T) where T <: Float32 - x = Float64(x) - z = x / v - zs = sqrt(1 + z^2) - n = zs + log(z / (1 + zs)) - coef = exp(-v * n) * sqrt(pi / (2*v)) / sqrt(zs) - p = inv(zs) - p2 = p*p - u0 = one(Float64) - u1 = 1 / 24 * evalpoly(p2, (3, -5)) - u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385)) - u3 = 1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425)) - u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725)) - return Float32(coef*evalpoly(-p/v, (u0, u1, u2, u3, u4))) - end - -function besselk_large_orders(v, x::T) where T <: AbstractFloat - x = Float64(x) - z = x / v - zs = sqrt(1 + z^2) - n = zs + log(z / (1 + zs)) - coef = exp(-v * n) * sqrt(pi / (2*v)) / sqrt(zs) - p = inv(zs) - p2 = p*p - u0 = one(Float64) - u1 = 1 / 24 * evalpoly(p2, (3, -5)) - u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385)) - u3 = 1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425)) - u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725)) - u5 = 1 / 6688604160 * evalpoly(p2, (1519035525, -49286948607, 284499769554, -614135872350, 566098157625, -188699385875)) - u6 = 1 / 4815794995200 * evalpoly(p2, (2757049477875, -127577298354750, 1050760774457901, -3369032068261860,5104696716244125, -3685299006138750, 1023694168371875)) - u7 = 1 / 115579079884800 * evalpoly(p2, (199689155040375, -12493049053044375, 138799253740521843, -613221795981706275, 1347119637570231525, -1570320948552481125, 931766432052080625, -221849150488590625)) - u8 = 1 / 22191183337881600 * evalpoly(p2, (134790179652253125, -10960565081605263000, 157768535329832893644, -914113758588905038248, 2711772922412520971550, -4513690624987320777000, 4272845805510421639500, -2152114239059719935000, 448357133137441653125)) - return coef*evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8)) -end - -#= -function besselk_large_orders(v, x::T) where T <: AbstractFloat - x = Float64(x) +function besselk_large_orders(v, x::T) where T <: Union{Float32, Float64, BigFloat} + S = promote_type(T, Float64) + x = S(x) z = x / v - zs = sqrt(1 + z^2) - n = zs + log(z / (1 + zs)) - coef = exp(-v * n) * sqrt(pi / (2*v)) / sqrt(zs) + zs = hypot(1, z) + n = zs + log(z) - log1p(zs) + coef = SQPIO2(S) * sqrt(inv(v)) * exp(-v*n) / sqrt(zs) p = inv(zs) - p2 = p*p - u0 = one(Float64) - u1 = 1 / 24 * evalpoly(p2, (3, -5)) - u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385)) - u3 = 1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425)) - u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725)) - u5 = 1 / 6688604160 * evalpoly(p2, (1519035525, -49286948607, 284499769554, -614135872350, 566098157625, -188699385875)) - u6 = 1 / 4815794995200 * evalpoly(p2, (2757049477875, -127577298354750, 1050760774457901, -3369032068261860,5104696716244125, -3685299006138750, 1023694168371875)) - u7 = 1 / 115579079884800 * evalpoly(p2, (199689155040375, -12493049053044375, 138799253740521843, -613221795981706275, 1347119637570231525, -1570320948552481125, 931766432052080625, -221849150488590625)) - u8 = 1 / 22191183337881600 * evalpoly(p2, (134790179652253125, -10960565081605263000, 157768535329832893644, -914113758588905038248, 2711772922412520971550, -4513690624987320777000, 4272845805510421639500, -2152114239059719935000, 448357133137441653125)) - u9 = 1 / 263631258054033408000 * evalpoly(p2, (6427469716717690265625, -659033454841709672064375, 11921080954211358275362500, -87432034049652400520788332, 334380732677827878090447630, -741743213039573443221773250, 992115946599792610768672500, -790370708270219620781737500, 345821892003106984030190625, -64041091111686478524109375)) - u10 = 1 / 88580102706155225088000 * evalpoly(p2, (9745329584487361980740625, -1230031256571145165088463750, 27299183373230345667273718125, -246750339886026017414509498824, 1177120360439828012193658602930, -3327704366990695147540934069220, 5876803711285273203043452095250, -6564241639632418015173104205000, 4513386761946134740461797128125, -1745632061522350031610173343750, 290938676920391671935028890625)) + p2 = v^2/fma(max(v,x), max(v,x), min(v,x)^2) - return coef*evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10)) + return T(coef*Uk_poly_Kn(p, v, p2, T)) end -=# - - - -#= -function besselk_large_orders(v, x::T) where T <: AbstractFloat +function besselk_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64, BigFloat} + S = promote_type(T, Float64) + x = S(x) z = x / v - zs = sqrt(1 + z^2) - - n = zs + log(z / (1 + zs)) - coef = exp(-v * n) * sqrt(pi / (2*v)) / sqrt(zs) - + zs = hypot(1, z) + n = zs + log(z) - log1p(zs) + coef = SQPIO2(S) * sqrt(inv(v)) * exp(-v*n + x) / sqrt(zs) p = inv(zs) - - u0 = one(x) - u1 = p / 24 * (3 - 5*p^2) * -1 / v - u2 = p^2 / 1152 * (81 - 462*p^2 + 385*p^4) / v^2 - u3 = p^3 / 414720 * (30375 - 369603 * p^2 + 765765*p^4 - 425425*p^6) * -1 / v^3 - u4 = p^4 / 39813120 * (4465125 - 94121676*p^2 + 349922430*p^4 - 446185740*p^6 + 185910725*p^8) / v^4 - u5 = p^5 / 6688604160 * (-188699385875*p^10 + 566098157625*p^8 - 614135872350*p^6 + 284499769554*p^4 - 49286948607*p^2 + 1519035525) * -1 / v^5 - u6 = p^6 / 4815794995200 * (1023694168371875*p^12 - 3685299006138750*p^10 + 5104696716244125*p^8 - 3369032068261860*p^6 + 1050760774457901*p^4 - 127577298354750*p^2 + 2757049477875) * 1 / v^6 - u7 = p^7 / 115579079884800 * (-221849150488590625*p^14 + 931766432052080625*p^12 - 1570320948552481125*p^10 + 1347119637570231525*p^8 - 613221795981706275*p^6 + 138799253740521843*p^4 - 12493049053044375*p^2 + 199689155040375) * -1 / v^7 - u8 = p^8 / 22191183337881600 * (448357133137441653125*p^16 - 2152114239059719935000*p^14 + 4272845805510421639500*p^12 - 4513690624987320777000*p^10 + 2711772922412520971550*p^8 - 914113758588905038248*p^6 + 157768535329832893644*p^4 - 10960565081605263000*p^2 + 134790179652253125) * 1 / v^8 - u9 = p^9 / 263631258054033408000 * (-64041091111686478524109375*p^18 + 345821892003106984030190625*p^16 - 790370708270219620781737500*p^14 + 992115946599792610768672500*p^12 - 741743213039573443221773250*p^10 + 334380732677827878090447630*p^8 - 87432034049652400520788332*p^6 + 11921080954211358275362500*p^4 - 659033454841709672064375*p^2 + 6427469716717690265625) * -1 / v^9 - u10 = p^10 / 88580102706155225088000 * (290938676920391671935028890625*p^20 - 1745632061522350031610173343750*p^18 + 4513386761946134740461797128125*p^16 - 6564241639632418015173104205000*p^14 + 5876803711285273203043452095250*p^12 - 3327704366990695147540934069220*p^10 + 1177120360439828012193658602930*p^8 - 246750339886026017414509498824*p^6 + 27299183373230345667273718125*p^4 - 1230031256571145165088463750*p^2 + 9745329584487361980740625) * 1 / v^10 - u11 = p^11 / 27636992044320430227456000 * (-1363312191078326248171914924296875*p^22 + 8997860461116953237934638500359375*p^20 - 25963913760458280822131901737878125*p^18 + 42936745153513007436411401865860625*p^16 - 44801790321820682384740638703320750*p^14 + 30589806122850866110941936529264950*p^12 - 13704902022868787041097596217578170*p^10 - + 3926191452593448964331218647028194*p^8 - 676389476843440422173605288596087*p^6 + 62011003282542082472252466220875*p^4 - 2321657500166464779660536015625*p^2 + 15237265774872558064250728125) * -1 / v^11 - u12 = p^12 / 39797268543821419527536640000 * (32426380464797989812768996474401171875*p^24 - 233469939346545526651936774615688437500*p^22 + 743739612850105971846901081692858843750*p^20 - 1378260730939829908037976894025628887500*p^18 + 1642838631056253395899341314188134178125*p^16 - 1314368459332124683504418467809129387000*p^14 - + 714528665351965363868467348116538170900*p^12 - 261201165596865827608687437905740780920*p^10 + 62055079517573388459132793029571876461*p^8 - 8958590476947726766450604043798559500*p^6 + 692277766674325563109617997687563750*p^4 - 21882222767154197351962677311437500*p^2 + 120907703923613748239829527671875) * 1 / v^12 - u13 = p^13 / 955134445051714068660879360000 * (-14020668045586884672121117629431460546875*p^26 + 109361210755577700442544717509565392265625*p^24 - 381190503845282445953419057314665534156250*p^22 + 782463969315283937856703223178540650343750*p^20 - - 1049095945162229046324321461816274931715625*p^18 + 962926533925253374359704288384340809260875*p^16 - 616410216242554698436702353237166008656700*p^14 + 274983827478138958623041508409195988431140*p^12 - 83924867223075156862785921508524155665245*p^10 - + 16843538631795795357786827345307534156063*p^8 - 2069933923586966756183324291232117362250*p^6 + 136735019134677724428035696765082813750*p^4 - 3698121486504259988094897605296209375*p^2 + 17438611142828905996129258798828125) * -1 / v^13 - u14 = p^14 / 45846453362482275295722209280000 * (13133360053559028970728309756597440972265625*p^28 - 110320224449895843354117801955418504167031250*p^26 + 417630985812040040477569028872405152769921875*p^24 - - 940627071986145750405920450097257805227812500*p^22 + 1401302601668131482630653233972052729323190625*p^20 - 1451823699927947602004297385351260623500372750*p^18 + 1070439683260179398514645952824098811025619475*p^16 - - 564850830044980230366140582063618983657685400*p^14 + 211477117385619365164298957821904115470535555*p^12 - 54857705817658080981995319669299096598482382*p^10 + 9440449669103391509091075981237243128469201*p^8 - - 1000503839668383458999731491914516564625300*p^6 + 57170953417612444837142230812990944671875*p^4 - 1338184074771428116079233795614103631250*p^2 + 5448320367052402487647812713291015625) * 1 / v^14 - u15 = p^15 / 9820310310243703368343697227776000000 * (- 59115551939078562227317999668652486368337724609375*p^30 + 532039967451707060045861997017872377315039521484375*p^28 - 2173722139119126857509156976742601742357422740234375*p^26 - + 5329871927856528282113994744458999865006055974609375*p^24 - 8735135969643867540297524795790262235822823374296875*p^22 + 10085018700249896522602267572484630409640764997271875*p^20 - 8420533422834140468835467666391400380550043688871875*p^18 - + 5136561256208409671660362298619778869859994724706875*p^16 - 2284251621937242886581917667066122422330060024456125*p^14 + 730367145705123976114617970888707594104468177381925*p^12 - 163359140754958502896104062604202934925448173291477*p^10 - + 24403480234538299231733883413666768614198435948125*p^8 - 2254933495791765108580529087615802250458013685625*p^6 + 112597271053778779753048514469995937998172890625*p^4 - 2303431987527333128955769182299845911305390625*p^2 + 8178936810213560828419581728001773291015625) * -1 / v^15 - - - out1 = ((u15 + u14 + u13 + u12 + u11 + u10 + u9) + u8) + u7 - out2 = ((out1 + u6 + u5) + u4) + u3 - out = ((out2 + u2) + u1) + u0 - - return coef*out - end - - - + p2 = v^2/fma(max(v,x), max(v,x), min(v,x)^2) - u11 = 1 / 27636992044320430227456000 * (-1363312191078326248171914924296875*p^22 + 8997860461116953237934638500359375*p^20 - 25963913760458280822131901737878125*p^18 + 42936745153513007436411401865860625*p^16 - 44801790321820682384740638703320750*p^14 + 30589806122850866110941936529264950*p^12 - 13704902022868787041097596217578170*p^10 - + 3926191452593448964331218647028194*p^8 - 676389476843440422173605288596087*p^6 + 62011003282542082472252466220875*p^4 - 2321657500166464779660536015625*p^2 + 15237265774872558064250728125) * -1 / v^11 - u12 = 1 / 39797268543821419527536640000 * (32426380464797989812768996474401171875*p^24 - 233469939346545526651936774615688437500*p^22 + 743739612850105971846901081692858843750*p^20 - 1378260730939829908037976894025628887500*p^18 + 1642838631056253395899341314188134178125*p^16 - 1314368459332124683504418467809129387000*p^14 - + 714528665351965363868467348116538170900*p^12 - 261201165596865827608687437905740780920*p^10 + 62055079517573388459132793029571876461*p^8 - 8958590476947726766450604043798559500*p^6 + 692277766674325563109617997687563750*p^4 - 21882222767154197351962677311437500*p^2 + 120907703923613748239829527671875) * 1 / v^12 - u13 = 1 / 955134445051714068660879360000 * (-14020668045586884672121117629431460546875*p^26 + 109361210755577700442544717509565392265625*p^24 - 381190503845282445953419057314665534156250*p^22 + 782463969315283937856703223178540650343750*p^20 - - 1049095945162229046324321461816274931715625*p^18 + 962926533925253374359704288384340809260875*p^16 - 616410216242554698436702353237166008656700*p^14 + 274983827478138958623041508409195988431140*p^12 - 83924867223075156862785921508524155665245*p^10 - + 16843538631795795357786827345307534156063*p^8 - 2069933923586966756183324291232117362250*p^6 + 136735019134677724428035696765082813750*p^4 - 3698121486504259988094897605296209375*p^2 + 17438611142828905996129258798828125) * -1 / v^13 - u14 = 1 / 45846453362482275295722209280000 * (13133360053559028970728309756597440972265625*p^28 - 110320224449895843354117801955418504167031250*p^26 + 417630985812040040477569028872405152769921875*p^24 - - 940627071986145750405920450097257805227812500*p^22 + 1401302601668131482630653233972052729323190625*p^20 - 1451823699927947602004297385351260623500372750*p^18 + 1070439683260179398514645952824098811025619475*p^16 - - 564850830044980230366140582063618983657685400*p^14 + 211477117385619365164298957821904115470535555*p^12 - 54857705817658080981995319669299096598482382*p^10 + 9440449669103391509091075981237243128469201*p^8 - - 1000503839668383458999731491914516564625300*p^6 + 57170953417612444837142230812990944671875*p^4 - 1338184074771428116079233795614103631250*p^2 + 5448320367052402487647812713291015625) * 1 / v^14 - u15 = 1 / 9820310310243703368343697227776000000 * (- 59115551939078562227317999668652486368337724609375*p^30 + 532039967451707060045861997017872377315039521484375*p^28 - 2173722139119126857509156976742601742357422740234375*p^26 - + 5329871927856528282113994744458999865006055974609375*p^24 - 8735135969643867540297524795790262235822823374296875*p^22 + 10085018700249896522602267572484630409640764997271875*p^20 - 8420533422834140468835467666391400380550043688871875*p^18 - + 5136561256208409671660362298619778869859994724706875*p^16 - 2284251621937242886581917667066122422330060024456125*p^14 + 730367145705123976114617970888707594104468177381925*p^12 - 163359140754958502896104062604202934925448173291477*p^10 - + 24403480234538299231733883413666768614198435948125*p^8 - 2254933495791765108580529087615802250458013685625*p^6 + 112597271053778779753048514469995937998172890625*p^4 - 2303431987527333128955769182299845911305390625*p^2 + 8178936810213560828419581728001773291015625) * -1 / v^15 -=# \ No newline at end of file + return T(coef*Uk_poly_Kn(p, v, p2, T)) +end \ No newline at end of file diff --git a/src/math_constants.jl b/src/math_constants.jl index 5420769..d54b0a5 100644 --- a/src/math_constants.jl +++ b/src/math_constants.jl @@ -1,10 +1,15 @@ const ONEOSQPI(::Type{BigFloat}) = big"5.6418958354775628694807945156077258584405E-1" const TWOOPI(::Type{BigFloat}) = big"6.3661977236758134307553505349005744813784E-1" +const SQPIO2(::Type{BigFloat}) = big"1.253314137315500251207882642405522626503493370304969158314961788171146827303924" +const SQ1O2PI(::Type{BigFloat}) = big"0.3989422804014326779399460599343818684758586311649346576659258296706579258993008" + const PIO4(::Type{Float64}) = .78539816339744830962 const THPIO4(::Type{Float64}) = 2.35619449019234492885 const SQ2OPI(::Type{Float64}) = .79788456080286535588 const TWOOPI(::Type{Float64}) = 0.6366197723675814 +const SQPIO2(::Type{Float64}) = 1.25331413731550025 +const SQ1O2PI(::Type{Float64}) = 0.3989422804014327 const PIO4(::Type{Float32}) = 0.78539816339744830962f0 diff --git a/test/besseli_test.jl b/test/besseli_test.jl index 90e10f6..72e74c4 100644 --- a/test/besseli_test.jl +++ b/test/besseli_test.jl @@ -64,3 +64,21 @@ i1x_32 = besseli1x.(Float32.(x32)) # test against SpecialFunctions.jl @test i1x_64 ≈ i1x64_SpecialFunctions @test i1x_32 ≈ i1x32_SpecialFunctions + + +# test for besseli +# test small arguments and order +m = 0:1:200; x = 0.1f0:0.5f0:90.0f0 +t = [besseli(m, x) for m in m, x in x] +@test t[10] isa Float32 +@test t ≈ Float32.([SpecialFunctions.besseli(m, x) for m in m, x in x]) + +#Float 64 +m = 0:1:200; x = 0.1:0.5:150.0 +t = [besseli(m, x) for m in m, x in x] + +@test t[10] isa Float64 +@test t ≈ [SpecialFunctions.besseli(m, x) for m in m, x in x] + +@test besselix(10, 2.0) ≈ SpecialFunctions.besselix(10, 2.0) + From 61906f2a0c8629b404986aea4552f30d89c6836d Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Fri, 21 Jan 2022 00:17:22 -0500 Subject: [PATCH 8/9] update docs --- src/Bessels.jl | 1 + src/U_polynomials.jl | 53 -------------------------------------------- src/besseli.jl | 36 +++++++++++++----------------- src/besselk.jl | 50 ++--------------------------------------- src/recurrence.jl | 32 ++++++++++++++++++++++++++ 5 files changed, 51 insertions(+), 121 deletions(-) create mode 100644 src/recurrence.jl diff --git a/src/Bessels.jl b/src/Bessels.jl index 6146595..0e3d3e9 100644 --- a/src/Bessels.jl +++ b/src/Bessels.jl @@ -35,6 +35,7 @@ include("Float128/constants.jl") include("math_constants.jl") include("U_polynomials.jl") +include("recurrence.jl") #include("parse.jl") #include("hankel.jl") diff --git a/src/U_polynomials.jl b/src/U_polynomials.jl index a9d3c7c..d17c97c 100644 --- a/src/U_polynomials.jl +++ b/src/U_polynomials.jl @@ -50,38 +50,8 @@ function Uk_poly_Kn(p, v, p2, ::Type{T}) where T <: BigFloat u10 = 1 / 88580102706155225088000 * evalpoly(p2, (9745329584487361980740625, -1230031256571145165088463750, 27299183373230345667273718125, -246750339886026017414509498824, 1177120360439828012193658602930, -3327704366990695147540934069220, 5876803711285273203043452095250, -6564241639632418015173104205000, 4513386761946134740461797128125, -1745632061522350031610173343750, 290938676920391671935028890625)) return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10)) end -function besselk_large_orders(v, x::T) where T <: AbstractFloat - x = Float64(x) - z = x / v - zs = sqrt(1 + z^2) - n = zs + log(z / (1 + zs)) - coef = exp(-v * n) * sqrt(pi / (2*v)) / sqrt(zs) - p = inv(zs) - p2 = p*p - u0 = one(Float64) - u1 = 1 / 24 * evalpoly(p2, (3, -5)) - u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385)) - u3 = 1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425)) - u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725)) - u5 = 1 / 6688604160 * evalpoly(p2, (1519035525, -49286948607, 284499769554, -614135872350, 566098157625, -188699385875)) - u6 = 1 / 4815794995200 * evalpoly(p2, (2757049477875, -127577298354750, 1050760774457901, -3369032068261860,5104696716244125, -3685299006138750, 1023694168371875)) - u7 = 1 / 115579079884800 * evalpoly(p2, (199689155040375, -12493049053044375, 138799253740521843, -613221795981706275, 1347119637570231525, -1570320948552481125, 931766432052080625, -221849150488590625)) - u8 = 1 / 22191183337881600 * evalpoly(p2, (134790179652253125, -10960565081605263000, 157768535329832893644, -914113758588905038248, 2711772922412520971550, -4513690624987320777000, 4272845805510421639500, -2152114239059719935000, 448357133137441653125)) - u9 = 1 / 263631258054033408000 * evalpoly(p2, (6427469716717690265625, -659033454841709672064375, 11921080954211358275362500, -87432034049652400520788332, 334380732677827878090447630, -741743213039573443221773250, 992115946599792610768672500, -790370708270219620781737500, 345821892003106984030190625, -64041091111686478524109375)) - u10 = 1 / 88580102706155225088000 * evalpoly(p2, (9745329584487361980740625, -1230031256571145165088463750, 27299183373230345667273718125, -246750339886026017414509498824, 1177120360439828012193658602930, -3327704366990695147540934069220, 5876803711285273203043452095250, -6564241639632418015173104205000, 4513386761946134740461797128125, -1745632061522350031610173343750, 290938676920391671935028890625)) - return coef*evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10)) -end -function besselk_large_orders(v, x::T) where T <: AbstractFloat - z = x / v - zs = sqrt(1 + z^2) - - n = zs + log(z / (1 + zs)) - coef = exp(-v * n) * sqrt(pi / (2*v)) / sqrt(zs) - - p = inv(zs) - u0 = one(x) u1 = p / 24 * (3 - 5*p^2) * -1 / v u2 = p^2 / 1152 * (81 - 462*p^2 + 385*p^4) / v^2 @@ -109,27 +79,4 @@ function besselk_large_orders(v, x::T) where T <: AbstractFloat + 5136561256208409671660362298619778869859994724706875*p^16 - 2284251621937242886581917667066122422330060024456125*p^14 + 730367145705123976114617970888707594104468177381925*p^12 - 163359140754958502896104062604202934925448173291477*p^10 + 24403480234538299231733883413666768614198435948125*p^8 - 2254933495791765108580529087615802250458013685625*p^6 + 112597271053778779753048514469995937998172890625*p^4 - 2303431987527333128955769182299845911305390625*p^2 + 8178936810213560828419581728001773291015625) * -1 / v^15 - - out1 = ((u15 + u14 + u13 + u12 + u11 + u10 + u9) + u8) + u7 - out2 = ((out1 + u6 + u5) + u4) + u3 - out = ((out2 + u2) + u1) + u0 - - return coef*out - end - - u11 = 1 / 27636992044320430227456000 * (-1363312191078326248171914924296875*p^22 + 8997860461116953237934638500359375*p^20 - 25963913760458280822131901737878125*p^18 + 42936745153513007436411401865860625*p^16 - 44801790321820682384740638703320750*p^14 + 30589806122850866110941936529264950*p^12 - 13704902022868787041097596217578170*p^10 - + 3926191452593448964331218647028194*p^8 - 676389476843440422173605288596087*p^6 + 62011003282542082472252466220875*p^4 - 2321657500166464779660536015625*p^2 + 15237265774872558064250728125) * -1 / v^11 - u12 = 1 / 39797268543821419527536640000 * (32426380464797989812768996474401171875*p^24 - 233469939346545526651936774615688437500*p^22 + 743739612850105971846901081692858843750*p^20 - 1378260730939829908037976894025628887500*p^18 + 1642838631056253395899341314188134178125*p^16 - 1314368459332124683504418467809129387000*p^14 - + 714528665351965363868467348116538170900*p^12 - 261201165596865827608687437905740780920*p^10 + 62055079517573388459132793029571876461*p^8 - 8958590476947726766450604043798559500*p^6 + 692277766674325563109617997687563750*p^4 - 21882222767154197351962677311437500*p^2 + 120907703923613748239829527671875) * 1 / v^12 - u13 = 1 / 955134445051714068660879360000 * (-14020668045586884672121117629431460546875*p^26 + 109361210755577700442544717509565392265625*p^24 - 381190503845282445953419057314665534156250*p^22 + 782463969315283937856703223178540650343750*p^20 - - 1049095945162229046324321461816274931715625*p^18 + 962926533925253374359704288384340809260875*p^16 - 616410216242554698436702353237166008656700*p^14 + 274983827478138958623041508409195988431140*p^12 - 83924867223075156862785921508524155665245*p^10 - + 16843538631795795357786827345307534156063*p^8 - 2069933923586966756183324291232117362250*p^6 + 136735019134677724428035696765082813750*p^4 - 3698121486504259988094897605296209375*p^2 + 17438611142828905996129258798828125) * -1 / v^13 - u14 = 1 / 45846453362482275295722209280000 * (13133360053559028970728309756597440972265625*p^28 - 110320224449895843354117801955418504167031250*p^26 + 417630985812040040477569028872405152769921875*p^24 - - 940627071986145750405920450097257805227812500*p^22 + 1401302601668131482630653233972052729323190625*p^20 - 1451823699927947602004297385351260623500372750*p^18 + 1070439683260179398514645952824098811025619475*p^16 - - 564850830044980230366140582063618983657685400*p^14 + 211477117385619365164298957821904115470535555*p^12 - 54857705817658080981995319669299096598482382*p^10 + 9440449669103391509091075981237243128469201*p^8 - - 1000503839668383458999731491914516564625300*p^6 + 57170953417612444837142230812990944671875*p^4 - 1338184074771428116079233795614103631250*p^2 + 5448320367052402487647812713291015625) * 1 / v^14 - u15 = 1 / 9820310310243703368343697227776000000 * (- 59115551939078562227317999668652486368337724609375*p^30 + 532039967451707060045861997017872377315039521484375*p^28 - 2173722139119126857509156976742601742357422740234375*p^26 - + 5329871927856528282113994744458999865006055974609375*p^24 - 8735135969643867540297524795790262235822823374296875*p^22 + 10085018700249896522602267572484630409640764997271875*p^20 - 8420533422834140468835467666391400380550043688871875*p^18 - + 5136561256208409671660362298619778869859994724706875*p^16 - 2284251621937242886581917667066122422330060024456125*p^14 + 730367145705123976114617970888707594104468177381925*p^12 - 163359140754958502896104062604202934925448173291477*p^10 - + 24403480234538299231733883413666768614198435948125*p^8 - 2254933495791765108580529087615802250458013685625*p^6 + 112597271053778779753048514469995937998172890625*p^4 - 2303431987527333128955769182299845911305390625*p^2 + 8178936810213560828419581728001773291015625) * -1 / v^15 =# diff --git a/src/besseli.jl b/src/besseli.jl index fe787a0..9041164 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -4,6 +4,9 @@ # Scaled modified Bessel functions of the first kind of order zero and one # besseli0x, besselix # +# (Scaled) Modified Bessel functions of the first kind of order nu +# besseli, besselix +# # Calculation of besseli0 is done in two branches using polynomial approximations [1] # # Branch 1: x < 7.75 @@ -35,11 +38,18 @@ # # Horner's scheme is then used to evaluate all polynomials. # ArbNumerics.jl is used as the reference bessel implementations with 75 digits. +# +# Calculation of besseli and besselkx can be done with downward recursion starting with +# besseli_{nu+1} and besseli_{nu}. Higher orders are determined by a uniform asymptotic +# expansion similar to besselk (see notes there) using Equation 10.41.3 [3]. +# # # [1] "Rational Approximations for the Modified Bessel Function of the First Kind # - I0(x) for Computations with Double Precision" by Pavel Holoborodko # [2] "Rational Approximations for the Modified Bessel Function of the First Kind -# - I1(x) for Computations with Double Precision" by Pavel Holoborodko +# - I1(x) for Computations with Double Precision" by Pavel Holoborodko +# [3] https://dlmf.nist.gov/10.41 + """ besseli0(x::T) where T <: Union{Float32, Float64} @@ -56,6 +66,7 @@ function besseli0(x::T) where T <: Union{Float32, Float64} return a * s end end + """ besseli0x(x::T) where T <: Union{Float32, Float64} @@ -73,6 +84,7 @@ function besseli0x(x::T) where T <: Union{Float32, Float64} return evalpoly(inv(x), besseli0_large_coefs(T)) / sqrt(x) end end + """ besseli1(x::T) where T <: Union{Float32, Float64} @@ -94,6 +106,7 @@ function besseli1(x::T) where T <: Union{Float32, Float64} end return z end + """ besseli1x(x::T) where T <: Union{Float32, Float64} @@ -136,6 +149,7 @@ function besseli(nu, x::T) where T <: Union{Float32, Float64, BigFloat} return besseli_large_orders(nu, x) end end + """ besselix(nu, x::T) where T <: Union{Float32, Float64} @@ -155,24 +169,6 @@ function besselix(nu, x::T) where T <: Union{Float32, Float64, BigFloat} return besseli_large_orders_scaled(nu, x) end end - - -@inline function down_recurrence(x, in, inp1, nu, branch) - # this prevents us from looping through large values of nu when the loop will always return zero - (iszero(in) || iszero(inp1)) && return zero(x) - (isinf(inp1) || isinf(inp1)) && return in - - inm1 = in - x2 = 2 / x - for n in branch:-1:nu+1 - a = x2 * n - inm1 = muladd(a, in, inp1) - inp1 = in - in = inm1 - end - return inm1 -end - function besseli_large_orders(v, x::T) where T <: Union{Float32, Float64, BigFloat} S = promote_type(T, Float64) x = S(x) @@ -196,4 +192,4 @@ function besseli_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64, p2 = v^2/fma(max(v,x), max(v,x), min(v,x)^2) return T(coef*Uk_poly_In(p, v, p2, T)) -end \ No newline at end of file +end diff --git a/src/besselk.jl b/src/besselk.jl index 82abe3b..1335ccf 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -55,34 +55,7 @@ # # However, calculation of these higher order U polynomials are tedious. These have been hand # calculated and somewhat crosschecked with symbolic math. There could be errors. They are listed -# here as a reference as higher orders are impossible to find while being needed for any meaningfully accurate calculation. - -# u0 = one(x) -# u1 = p / 24 * (3 - 5*p^2) * -1 / v -# u2 = p^2 / 1152 * (81 - 462*p^2 + 385*p^4) / v^2 -# u3 = p^3 / 414720 * (30375 - 369603 * p^2 + 765765*p^4 - 425425*p^6) * -1 / v^3 -# u4 = p^4 / 39813120 * (4465125 - 94121676*p^2 + 349922430*p^4 - 446185740*p^6 + 185910725*p^8) / v^4 -# u5 = p^5 / 6688604160 * (-188699385875*p^10 + 566098157625*p^8 - 614135872350*p^6 + 284499769554*p^4 - 49286948607*p^2 + 1519035525) * -1 / v^5 -# u6 = p^6 / 4815794995200 * (1023694168371875*p^12 - 3685299006138750*p^10 + 5104696716244125*p^8 - 3369032068261860*p^6 -# + 1050760774457901*p^4 - 127577298354750*p^2 + 2757049477875) * 1 / v^6 -# u7 = p^7 / 115579079884800 * (-221849150488590625*p^14 + 931766432052080625*p^12 - 1570320948552481125*p^10 + 1347119637570231525*p^8 -# - 613221795981706275*p^6 + 138799253740521843*p^4 - 12493049053044375*p^2 + 199689155040375) * -1 / v^7 -# u8 = p^8 / 22191183337881600 * (448357133137441653125*p^16 - 2152114239059719935000*p^14 + 4272845805510421639500*p^12 - 4513690624987320777000*p^10 -# + 2711772922412520971550*p^8 - 914113758588905038248*p^6 + 157768535329832893644*p^4 - 10960565081605263000*p^2 + 134790179652253125) * 1 / v^8 -# u9 = p^9 / 263631258054033408000 * (-64041091111686478524109375*p^18 + 345821892003106984030190625*p^16 - 790370708270219620781737500*p^14 -# + 992115946599792610768672500*p^12 - 741743213039573443221773250*p^10 + 334380732677827878090447630*p^8 - 87432034049652400520788332*p^6 -# + 11921080954211358275362500*p^4 - 659033454841709672064375*p^2 + 6427469716717690265625) * -1 / v^9 -# u10 = p^10 / 88580102706155225088000 * (290938676920391671935028890625*p^20 - 1745632061522350031610173343750*p^18 + 4513386761946134740461797128125*p^16 -# - 6564241639632418015173104205000*p^14 + 5876803711285273203043452095250*p^12 - 3327704366990695147540934069220*p^10 + 1177120360439828012193658602930*p^8 -# - 246750339886026017414509498824*p^6 + 27299183373230345667273718125*p^4 - 1230031256571145165088463750*p^2 + 9745329584487361980740625) * 1 / v^10 -# u11 = p^11 / 27636992044320430227456000 * (-1363312191078326248171914924296875*p^22 + 8997860461116953237934638500359375*p^20 - 25963913760458280822131901737878125*p^18 + 42936745153513007436411401865860625*p^16 - 44801790321820682384740638703320750*p^14 + 30589806122850866110941936529264950*p^12 - 13704902022868787041097596217578170*p^10 -# + 3926191452593448964331218647028194*p^8 - 676389476843440422173605288596087*p^6 + 62011003282542082472252466220875*p^4 - 2321657500166464779660536015625*p^2 + 15237265774872558064250728125) * -1 / v^11 -# u12 = p^12 / 39797268543821419527536640000 * (32426380464797989812768996474401171875*p^24 - 233469939346545526651936774615688437500*p^22 + 743739612850105971846901081692858843750*p^20 - 1378260730939829908037976894025628887500*p^18 + 1642838631056253395899341314188134178125*p^16 - 1314368459332124683504418467809129387000*p^14 -# + 714528665351965363868467348116538170900*p^12 - 261201165596865827608687437905740780920*p^10 + 62055079517573388459132793029571876461*p^8 - 8958590476947726766450604043798559500*p^6 + 692277766674325563109617997687563750*p^4 - 21882222767154197351962677311437500*p^2 + 120907703923613748239829527671875) * 1 / v^12 -# u13 = p^13 / 955134445051714068660879360000 * (-14020668045586884672121117629431460546875*p^26 + 109361210755577700442544717509565392265625*p^24 - 381190503845282445953419057314665534156250*p^22 + 782463969315283937856703223178540650343750*p^20 -# - 1049095945162229046324321461816274931715625*p^18 + 962926533925253374359704288384340809260875*p^16 - 616410216242554698436702353237166008656700*p^14 + 274983827478138958623041508409195988431140*p^12 - 83924867223075156862785921508524155665245*p^10 -# + 16843538631795795357786827345307534156063*p^8 - 2069933923586966756183324291232117362250*p^6 + 136735019134677724428035696765082813750*p^4 - 3698121486504259988094897605296209375*p^2 + 17438611142828905996129258798828125) * -1 / v^13 -# +# in src/U_polynomials.jl as a reference as higher orders are impossible to find while being needed for any meaningfully accurate calculation. # For large orders these formulas will converge much faster than using upward recurrence. # # @@ -97,6 +70,7 @@ # [6] Abramowitz, Milton, and Irene A. Stegun, eds. Handbook of mathematical functions with formulas, graphs, and mathematical tables. # Vol. 55. US Government printing office, 1964. # + """ besselk0(x::T) where T <: Union{Float32, Float64} @@ -206,25 +180,6 @@ function besselkx(nu::Int, x::T) where T <: Union{Float32, Float64} return besselk_large_orders_scaled(nu, x) end end - -@inline function up_recurrence(x, k0, k1, nu) - nu == 0 && return k0 - nu == 1 && return k1 - - # this prevents us from looping through large values of nu when the loop will always return zero - (iszero(k0) || iszero(k1)) && return zero(x) - - k2 = k0 - x2 = 2 / x - for n in 1:nu-1 - a = x2 * n - k2 = muladd(a, k1, k0) - k0 = k1 - k1 = k2 - end - return k2 -end - function besselk_large_orders(v, x::T) where T <: Union{Float32, Float64, BigFloat} S = promote_type(T, Float64) x = S(x) @@ -237,7 +192,6 @@ function besselk_large_orders(v, x::T) where T <: Union{Float32, Float64, BigFlo return T(coef*Uk_poly_Kn(p, v, p2, T)) end - function besselk_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64, BigFloat} S = promote_type(T, Float64) x = S(x) diff --git a/src/recurrence.jl b/src/recurrence.jl new file mode 100644 index 0000000..e42204a --- /dev/null +++ b/src/recurrence.jl @@ -0,0 +1,32 @@ +@inline function down_recurrence(x, in, inp1, nu, branch) + # this prevents us from looping through large values of nu when the loop will always return zero + (iszero(in) || iszero(inp1)) && return zero(x) + (isinf(inp1) || isinf(inp1)) && return in + + inm1 = in + x2 = 2 / x + for n in branch:-1:nu+1 + a = x2 * n + inm1 = muladd(a, in, inp1) + inp1 = in + in = inm1 + end + return inm1 +end +@inline function up_recurrence(x, k0, k1, nu) + nu == 0 && return k0 + nu == 1 && return k1 + + # this prevents us from looping through large values of nu when the loop will always return zero + (iszero(k0) || iszero(k1)) && return zero(x) + + k2 = k0 + x2 = 2 / x + for n in 1:nu-1 + a = x2 * n + k2 = muladd(a, k1, k0) + k0 = k1 + k1 = k2 + end + return k2 +end \ No newline at end of file From ef12079e40c5e82b2fd199bfe6a4df72194a1e63 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Fri, 21 Jan 2022 00:22:58 -0500 Subject: [PATCH 9/9] fix test coverage --- src/Bessels.jl | 1 - src/hankel.jl | 2 ++ src/math_constants.jl | 4 ++-- src/parse.jl | 7 ------- test/besseli_test.jl | 3 ++- 5 files changed, 6 insertions(+), 11 deletions(-) delete mode 100644 src/parse.jl diff --git a/src/Bessels.jl b/src/Bessels.jl index 0e3d3e9..729cabd 100644 --- a/src/Bessels.jl +++ b/src/Bessels.jl @@ -36,7 +36,6 @@ include("Float128/constants.jl") include("math_constants.jl") include("U_polynomials.jl") include("recurrence.jl") -#include("parse.jl") #include("hankel.jl") end diff --git a/src/hankel.jl b/src/hankel.jl index 5912451..4b44419 100644 --- a/src/hankel.jl +++ b/src/hankel.jl @@ -1,3 +1,4 @@ +#= Aren't tested yet function besselh(nu::Float64, k::Integer, x::AbstractFloat) if k == 1 return complex(besselj(nu, x), bessely(nu, x)) @@ -10,3 +11,4 @@ end hankelh1(nu, z) = besselh(nu, 1, z) hankelh2(nu, z) = besselh(nu, 2, z) +=# \ No newline at end of file diff --git a/src/math_constants.jl b/src/math_constants.jl index d54b0a5..27d1a56 100644 --- a/src/math_constants.jl +++ b/src/math_constants.jl @@ -1,7 +1,7 @@ const ONEOSQPI(::Type{BigFloat}) = big"5.6418958354775628694807945156077258584405E-1" const TWOOPI(::Type{BigFloat}) = big"6.3661977236758134307553505349005744813784E-1" -const SQPIO2(::Type{BigFloat}) = big"1.253314137315500251207882642405522626503493370304969158314961788171146827303924" -const SQ1O2PI(::Type{BigFloat}) = big"0.3989422804014326779399460599343818684758586311649346576659258296706579258993008" +#const SQPIO2(::Type{BigFloat}) = big"1.253314137315500251207882642405522626503493370304969158314961788171146827303924" +#const SQ1O2PI(::Type{BigFloat}) = big"0.3989422804014326779399460599343818684758586311649346576659258296706579258993008" const PIO4(::Type{Float64}) = .78539816339744830962 diff --git a/src/parse.jl b/src/parse.jl deleted file mode 100644 index 0002451..0000000 --- a/src/parse.jl +++ /dev/null @@ -1,7 +0,0 @@ -@inline function Float128(str::S) where {S<:AbstractString} - return Float128(BigFloat(str)) -end - -macro f128_str(val::AbstractString) - :(Float128($val)) -end \ No newline at end of file diff --git a/test/besseli_test.jl b/test/besseli_test.jl index 72e74c4..7a015f6 100644 --- a/test/besseli_test.jl +++ b/test/besseli_test.jl @@ -81,4 +81,5 @@ t = [besseli(m, x) for m in m, x in x] @test t ≈ [SpecialFunctions.besseli(m, x) for m in m, x in x] @test besselix(10, 2.0) ≈ SpecialFunctions.besselix(10, 2.0) - +@test besselix(100, 14.0) ≈ SpecialFunctions.besselix(100, 14.0) +@test besselix(120, 504.0) ≈ SpecialFunctions.besselix(120, 504.0)