pole hopping for 1F0

src/weniger.jl

@@ -60,24 +60,71 @@ function pFqweniger(α::Tuple{T1}, ::Tuple{}, z::T2; kmax::Int = KMAX) where {T1
     if norm(z) < eps(real(T)) || norm(α) < eps(absα)
         return one(T)
     end
+    z2 = z*z
     μlo = α
-    μhi = inv(2-(α+1)*z)
     Rlo = one(T)
-    Rhi = Rlo + 2z*μhi*μlo
-    if norm(α+1) < eps(absα+1)
-        return Rhi
-    end
-    μhi *= α+1
-    k = 1
-    z2 = z*z
-    while k < kmax && errcheck(Rlo, Rhi, 8eps(real(T)))
+    k = 0
+    if iszero(2-(α+1)*z)
+        μhi = inv(2-(α+1)*z)
+        Rhi = Rlo + 2z*μhi*μlo
+        μhi *= α+1
+        k += 1
         μhi, μlo = inv((2k+1)*(2-z) - (k-α)*z2*μhi), μhi
-        Rhi, Rlo = ((2k+1)*(2-z)*Rhi - (k-α)*z2*Rlo*μlo)*μhi, Rhi
+        Rhi, Rlo = Rlo - (2k+1)*(2-z)*2*α/(z*(α+1)*(k-α)), Rhi
         if norm(α+k+1) < eps(absα+k+1)
             return Rhi
         end
         μhi *= α+k+1
         k += 1
+        μhi, μlo = inv((2k+1)*(2-z) - (k-α)*z2*μhi), μhi
+        Rhi, Rlo = Rhi + 2z*(k-α)*α*(α+2)/((1-α)*(α+1))*μhi, Rhi
+        if norm(α+k+1) < eps(absα+k+1)
+            return Rhi
+        end
+        μhi *= α+k+1
+        k += 1
+    else
+        μhi = inv(2-(α+1)*z)
+        Rhi = Rlo + 2z*μhi*μlo
+        if norm(α+1) < eps(absα+1)
+            return Rhi
+        end
+        μhi *= α+1
+        k += 1
+    end
+    while k < kmax && errcheck(Rlo, Rhi, 8eps(real(T)))
+        if iszero((2k+1)*(2-z) - (k-α)*z2*μhi)
+            μhi, μlo = inv((2k+1)*(2-z) - (k-α)*z2*μhi), μhi
+            cst = (2k+1)*(2-z)*Rhi - (k-α)*z2*Rlo*μlo
+            Rhi, Rlo = cst*μhi, Rhi
+            if norm(α+k+1) < eps(absα+k+1)
+                return Rhi
+            end
+            μhi *= α+k+1
+            k += 1
+            μhi, μlo = inv((2k+1)*(2-z) - (k-α)*z2*μhi), μhi # ✓
+            Rhi, Rlo = Rlo - (2k+1)*(2-z)*cst/((k-α)*z2*(α+k)), Rhi
+            if norm(α+k+1) < eps(absα+k+1)
+                return Rhi
+            end
+            μhi *= α+k+1
+            k += 1
+            μhi, μlo = inv((2k+1)*(2-z) - (k-α)*z2*μhi), μhi
+            Rhi, Rlo = Rhi + cst*(k-α)*(α+k)/((k-1-α)*(α+k-1))*μhi, Rhi
+            if norm(α+k+1) < eps(absα+k+1)
+                return Rhi
+            end
+            μhi *= α+k+1
+            k += 1
+        else
+            μhi, μlo = inv((2k+1)*(2-z) - (k-α)*z2*μhi), μhi
+            Rhi, Rlo = ((2k+1)*(2-z)*Rhi - (k-α)*z2*Rlo*μlo)*μhi, Rhi
+            if norm(α+k+1) < eps(absα+k+1)
+                return Rhi
+            end
+            μhi *= α+k+1
+            k += 1
+        end
     end
     k < kmax || @warn "Rational approximation to "*pFq2string(Val(1), Val{0}())*" reached the maximum type of ("*string(kmax, ", ", kmax)*")."
     return isfinite(Rhi) ? Rhi : Rlo

test/runtests.jl

@@ -387,14 +387,26 @@ end
             @test pFqdrummond((α, ), (), z) ≈ S(pFq((T(α), ), (), T(z))) atol=atol rtol=rtol
             @test pFqweniger((α, ), (), z) ≈ S(pFq((T(α), ), (), T(z))) atol=atol rtol=rtol
         end
+        #for α = S(-1.5):S(0.5):S(1.5)
+        #    z = 2/(α+1)
+        #    @test pFqdrummond((α, ), (), z) ≈ S(pFq((T(α), ), (), 2/(T(α)+1))) atol=atol rtol=rtol
+        #    @test pFqweniger((α, ), (), z) ≈ S(pFq((T(α), ), (), 2/(T(α)+1))) atol=atol rtol=rtol
+        #end
         atol *= 2
         rtol *= 2
         for α in S(-1.5):S(0.5):S(1.5), x in S(-0.75):S(0.25):S(0.25), y in S(-0.75):S(0.25):S(0.25)
             z = complex(x, y)
             @test pFqdrummond((α, ), (), z) ≈ CS(pFq((T(α), ), (), CT(z))) atol=atol rtol=rtol
             @test pFqweniger((α, ), (), z) ≈ CS(pFq((T(α), ), (), CT(z))) atol=atol rtol=rtol
         end
+        z = S(2.0)^-104
+        α = S(0)+S(1)*im
+        @test real(pFq((α, ), (), z)) ≈ real(exp(-α*log1p(-z)))
+        @test imag(pFq((α, ), (), z)) ≈ imag(exp(-α*log1p(-z))) # High relative accuracy in the imaginary part is precisely why we don't use (1-z)^-α
     end
+    z = 0.43019497816657903
+    α = 8.5
+    @test pFqweniger((α, ), (), z) ≈ pFq((α, ), (), z)
 end
 
 @testset "₀F₁" begin