Merge pull request #10 from heltonmc/continued_fraction

Use continued fractions for besseli for medium sized orders 2<nu<100

Author: heltonmc

src/U_polynomials.jl

@@ -2,37 +2,29 @@ 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))
+    return evalpoly(-p/v, (u0, u1, u2))
 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))
+    return evalpoly(-p/v, (u0, u1, u2, u3, u4))
 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))
+    return evalpoly(-p/v, (u0, u1, u2, u3, u4))
 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))
+    return evalpoly(-p/v, (u0, u1, u2))
  end
 
 #=

src/besseli.jl

@@ -136,17 +136,18 @@ end
 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}
+function besseli(nu, x::T) where T <: Union{Float32, Float64}
     nu == 0 && return besseli0(x)
     nu == 1 && return besseli1(x)
-
-    branch = 60
-    if nu < branch
-        inp1 = besseli_large_orders(branch + 1, x)
-        in = besseli_large_orders(branch, x)
-        return down_recurrence(x, in, inp1, nu, branch)
+    
+    if x > maximum((T(30), nu^2 / 4))
+        return T(besseli_large_argument(nu, x))
+    elseif x <= 2 * sqrt(nu + 1)
+        return T(besseli_small_arguments(nu, x))
+    elseif nu < 100
+        return T(_besseli_continued_fractions(nu, x))
     else
-        return besseli_large_orders(nu, x)
+        return T(besseli_large_orders(nu, x))
     end
 end
 
@@ -156,20 +157,21 @@ end
 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}
+function besselix(nu, x::T) where T <: Union{Float32, Float64}
     nu == 0 && return besseli0x(x)
     nu == 1 && return besseli1x(x)
 
-    branch = 60
-    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)
+    if x > maximum((T(30), nu^2 / 4))
+        return T(besseli_large_argument_scaled(nu, x))
+    elseif x <= 2 * sqrt(nu + 1)
+        return T(besseli_small_arguments(nu, x)) * exp(-x)
+    elseif nu < 100
+        return T(_besseli_continued_fractions_scaled(nu, x))
     else
         return besseli_large_orders_scaled(nu, x)
     end
 end
-function besseli_large_orders(v, x::T) where T <: Union{Float32, Float64, BigFloat}
+function besseli_large_orders(v, x::T) where T <: Union{Float32, Float64}
     S = promote_type(T, Float64)
     x = S(x)
     z = x / v
@@ -179,9 +181,9 @@ function besseli_large_orders(v, x::T) where T <: Union{Float32, Float64, BigFlo
     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))
+    return coef*Uk_poly_In(p, v, p2, T)
 end
-function besseli_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64, BigFloat}
+function besseli_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64}
     S = promote_type(T, Float64)
     x = S(x)
     z = x / v
@@ -193,3 +195,94 @@ function besseli_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64,
 
     return T(coef*Uk_poly_In(p, v, p2, T))
 end
+function _besseli_continued_fractions(nu, x::T) where T
+    S = promote_type(T, Float64)
+    xx = S(x)
+    knu, knum1 = up_recurrence(xx, besselk0(xx), besselk1(xx), nu)
+    # if knu or knum1 is zero then besseli will likely overflow
+    (iszero(knu) || iszero(knum1)) && return throw(DomainError(x, "Overflow error"))
+    return 1 / (x * (knum1 + knu / steed(nu, x)))
+end
+function _besseli_continued_fractions_scaled(nu, x::T) where T
+    S = promote_type(T, Float64)
+    xx = S(x)
+    knu, knum1 = up_recurrence(xx, besselk0x(xx), besselk1x(xx), nu)
+    # if knu or knum1 is zero then besseli will likely overflow
+    (iszero(knu) || iszero(knum1)) && return throw(DomainError(x, "Overflow error"))
+    return 1 / (x * (knum1 + knu / steed(nu, x)))
+end
+function steed(n, x::T) where T
+    MaxIter = 1000
+    xinv = inv(x)
+    xinv2 = 2 * xinv
+    d = x / (n + n)
+    a = d
+    h = a
+    b = muladd(2, n, 2) * xinv
+    for _ in 1:MaxIter
+        d = inv(b + d)
+        a *= muladd(b, d, -1)
+        h = h + a
+        b = b + xinv2
+        abs(a / h) <= eps(T) && break
+    end
+    return h
+end
+function besseli_large_argument(v, z::T) where T
+    MaxIter = 1000
+    a = exp(z / 2)
+    coef = a / sqrt(2 * T(pi) * z)
+    fv2 = 4 * v^2
+    term = one(T)
+    res = term
+    s = -term
+    for i in 1:MaxIter
+        i = T(i)
+        offset = muladd(2, i, -1)
+        term *= T(0.125) * muladd(offset, -offset, fv2) / (z * i)
+        res = muladd(term, s, res)
+        s = -s
+        abs(term) <= eps(T) && break
+    end
+    return res * coef * a
+end
+function besseli_large_argument_scaled(v, z::T) where T
+    MaxIter = 1000
+    coef = inv(sqrt(2 * T(pi) * z))
+    fv2 = 4 * v^2
+    term = one(T)
+    res = term
+    s = -term
+    for i in 1:MaxIter
+        i = T(i)
+        offset = muladd(2, i, -1)
+        term *= T(0.125) * muladd(offset, -offset, fv2) / (z * i)
+        res = muladd(term, s, res)
+        s = -s
+        abs(term) <= eps(T) && break
+    end
+    return res * coef
+end
+
+function besseli_small_arguments(v, z::T) where T
+    S = promote_type(T, Float64)
+    x = S(z)
+    if v < 20
+        coef = (x / 2)^v / factorial(v)
+    else
+        vinv = inv(v)
+        coef = sqrt(vinv / (2 * π)) * MathConstants.e^(v * (log(x / (2 * v)) + 1)) 
+        coef *= evalpoly(vinv, (1, -1/12, 1/288,  139/51840, -571/2488320, -163879/209018880, 5246819/75246796800, 534703531/902961561600))
+    end
+
+    MaxIter = 1000
+    out = one(S)
+    zz = x^2 / 4
+    a = one(S)
+    for k in 1:MaxIter
+        a *= zz / (k * (k + v))
+        out += a
+        a <= eps(T) && break
+    end
+    return coef * out
+end

src/besselk.jl

@@ -161,7 +161,7 @@ Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``.
 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)
+        return up_recurrence(x, besselk0(x), besselk1(x), nu)[1]
     else
         return besselk_large_orders(nu, x)
     end
@@ -175,7 +175,7 @@ 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}
     T == Float32 ? branch = 20 : branch = 50
     if nu < branch
-        return up_recurrence(x, besselk0x(x), besselk1x(x), nu)
+        return up_recurrence(x, besselk0x(x), besselk1x(x), nu)[1]
     else
         return besselk_large_orders_scaled(nu, x)
     end

src/recurrence.jl

@@ -1,3 +1,5 @@
+# no longer used for besseli but could be used in future for Jn, Yn
+#=
 @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)
@@ -13,6 +15,7 @@
     end
     return inm1
 end
+=#
 @inline function up_recurrence(x, k0, k1, nu)
     nu == 0 && return k0
     nu == 1 && return k1
@@ -28,5 +31,5 @@ end
         k0 = k1
         k1 = k2
     end
-    return k2
+    return k2, k0
 end

test/besseli_test.jl

@@ -68,18 +68,26 @@ i1x_32 = besseli1x.(Float32.(x32))
 
 # 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])
+m = 0:1:200; x = 0.5f0:0.5f0:90.0f0
+@test besseli(10, 1.0f0) isa Float32
+@test besseli(2, 80.0f0) isa Float32
+@test besseli(112, 80.0f0) isa Float32
+
+for m in m, x in x
+    @test besseli(m, x) ≈ Float32(SpecialFunctions.besseli(m, x))
+end
 
 #Float 64
 m = 0:1:200; x = 0.1:0.5:150.0
+@test besseli(10, 1.0) isa Float64
+@test besseli(2, 80.0) isa Float64
+@test besseli(112, 80.0) isa Float64
 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)
-@test besselix(100, 14.0) ≈ SpecialFunctions.besselix(100, 14.0)
-@test besselix(120, 504.0) ≈ SpecialFunctions.besselix(120, 504.0)
+t = [besselix(m, x) for m in m, x in x]
+@test t[10] isa Float64
+@test t ≈ [SpecialFunctions.besselix(m, x) for m in m, x in x]
+