Merge pull request #312 from JuliaMath/an/En_gamma

Fix evaluation of continued fraction for En_gamma

Author: andreasnoack

src/expint.jl

@@ -214,17 +214,19 @@ function En_cf_gamma(ν::Number, z::Number, n::Int=1000)
     ϵ = 10*eps(real(B))
     scale = sqrt(floatmax(real(A)))
 
-    iters = 1
-    j = 1
-    for i = 2:n
+    iters = 0
+    for i = 1:n
         iters += 1
 
-        term = iseven(i) ? -(i÷2 - 1 - ν)*z : z*(i÷2)
-        A, Aprev = (i - ν)*A + term * Aprev, A
-        B, Bprev = (i - ν)*B + term * Bprev, B
+        a = iseven(i) ? (i÷2)*z : -((i + 1)÷2 - ν)*z
+        b = (i + 1 - ν)
+
+        A, Aprev = b*A + a*Aprev, A
+        B, Bprev = b*B + a*Bprev, B
 
-        conv = abs(Aprev*B - A*Bprev) < ϵ*abs(B*Bprev)
-        conv && break
+        if abs(Aprev*B - A*Bprev) < ϵ*abs(Aprev*B)
+            break
+        end
 
         if max(abs(real(A)), abs(imag(A))) > scale
             A     /= scale
@@ -436,6 +438,7 @@ function _expint(ν::Number, z::Number, niter::Int=1000, ::Val{expscaled}=Val{fa
         else
             g, cf, _ = zero(z), En_cf_nogamma(ν, z, niter)...
         end
+
         g != 0 && (g *= gmult)
         cf = expscaled ? cf : En_safeexpmult(-z, cf)
         return g + cf

test/expint.jl

@@ -123,6 +123,48 @@ using Base.MathConstants
         for (ν, z, E) in zip(complex_orders, complex_inputs, sympy_complex_outputs)
             abs(ν) ≤ 50 ? (@test expint(ν, z) ⩭ E) : (@test_throws ArgumentError expint(ν, z) ⩭ E)
         end
+
+        xs = range(2.0, stop=5.0, step=0.05)
+        ys = [0.05702612399289205, 0.05308665097847697, 0.04944294741593074, 
+            0.04607036967463369, 0.04294659589019342, 0.04005137577047593, 
+            0.03736631128005022, 0.03487466387112418, 0.03256118461138553, 
+            0.03041196412241836, 0.02841429970928976, 0.02655657745053252, 
+            0.02482816734246761, 0.0232193298641301, 0.02172113255824923, 
+            0.02032537541727814, 0.01902452402744482, 0.01781164956316038, 
+            0.01668037484292225, 0.01562482575941789, 0.01463958748361088, 
+            0.01371966491744429, 0.01286044693430924, 0.01205767400216573, 
+            0.01130740883247835, 0.01060600974003776, 0.00995010643520129, 
+            0.00933657800187881, 0.00876253284236996, 0.00822529039448117, 
+            0.007722364447687571, 0.007251447903858146, 0.006810398844575473, 
+            0.006397227781646947, 0.006010085980275144, 0.005647254755746198, 
+            0.005307135654593942, 0.004988241440163916, 0.004689187810475, 
+            0.004408685783377122, 0.004145534690336333, 0.003898615725833893, 
+            0.003666886004423625, 0.003449373082020276, 0.00324516990205213, 
+            0.003053430130755408, 0.002873363849164126, 0.0027042335722962, 
+            0.002545350568691675, 0.002396071455853417, 0.002255795049301819, 
+            0.002123959444908601, 0.00200003931594166, 0.001883543407852843, 
+            0.001774012215290864, 0.001671015827136935, 0.00157415192655548, 
+            0.001483043934137846, 0.001397339283204469, 0.001316707817230031, 
+            0.001240840300175119]
+
+        @testset "En_cf_nogamma" begin
+            for (x, y) in zip(xs, ys)
+                cf, iter = SpecialFunctions.En_cf_nogamma(1 / 2, x)
+                @test cf*exp(-x) ≈ y
+            end
+        end
+        @testset "En_cf_gamma" begin
+            for (x, y) in zip(xs, ys)
+                Γ, cf, iter = SpecialFunctions.En_cf_gamma(1 / 2, x)
+                @test Γ + cf*exp(-x) ≈ y
+            end
+        end
+        @testset "En_expand_origin" begin
+            for (x, y) in zip(xs, ys)
+                res = SpecialFunctions.En_expand_origin(1 / 2, x, 1000)
+                @test res ≈ y
+            end
+        end
     end
 
     # issue #259: positive vs. negative integer order for small/larger args