pole hopping for 0F0

src/drummond.jl

@@ -12,16 +12,56 @@ function pFqdrummond(::Tuple{}, ::Tuple{}, z::T; kmax::Int = KMAX) where T
         return one(T)
     end
     μlo = one(z)
-    μhi = inv(2-z)
     Tlo = one(z)
-    Thi = Tlo + 2z*μhi
-    μhi, μlo = inv(3-z + z*μhi), μhi
-    Thi, Tlo = ((3-z)*Thi + z*(Tlo+z)*μlo)*μhi, Thi
-    k = 2
-    while k < kmax && errcheck(Tlo, Thi, 8eps(real(T)))
+    k = 0
+    if iszero(2-z)
+        μhi = inv(2-z)
+        Thi = Tlo + 2z*μhi
+        k += 1
+        μhi, μlo = inv(k+2-z + k*z*μhi), μhi
+        Thi, Tlo = 6-z + Tlo, Thi
+        k += 1
         μhi, μlo = inv(k+2-z + k*z*μhi), μhi
-        Thi, Tlo = ((k+2-z)*Thi + k*z*Tlo*μlo)*μhi, Thi
+        Thi, Tlo = ((k+2-z)*Thi + k*z*2)*μhi, Thi
         k += 1
+    else
+        μhi = inv(2-z)
+        Thi = Tlo + 2z*μhi
+        k += 1
+        if iszero(3-z + z*μhi)
+            μhi, μlo = inv(3-z + z*μhi), μhi
+            cst = (3-z)*Thi + z*(Tlo+z)*μlo
+            Thi, Tlo = cst*μhi, Thi
+            k += 1
+            μhi, μlo = inv(k+2-z + k*z*μhi), μhi
+            Thi, Tlo = (k+2-z)*cst/(k*z) + Tlo, Thi
+            k += 1
+            μhi, μlo = inv(k+2-z + k*z*μhi), μhi
+            Thi, Tlo = Thi + (k*cst/T(k-1))*μhi, Thi
+            k += 1
+        else
+            μhi, μlo = inv(3-z + z*μhi), μhi
+            Thi, Tlo = ((3-z)*Thi + z*(Tlo+z)*μlo)*μhi, Thi
+            k += 1
+        end
+    end
+    while k < kmax && errcheck(Tlo, Thi, 8eps(real(T)))
+        if iszero(k+2-z + k*z*μhi)
+            μhi, μlo = inv(k+2-z + k*z*μhi), μhi
+            cst = ((k+2-z)*Thi + k*z*Tlo*μlo)
+            Thi, Tlo = cst*μhi, Thi
+            k += 1
+            μhi, μlo = inv(k+2-z + k*z*μhi), μhi
+            Thi, Tlo = (k+2-z)*cst/(k*z) + Tlo, Thi
+            k += 1
+            μhi, μlo = inv(k+2-z + k*z*μhi), μhi
+            Thi, Tlo = Thi + (k*cst/T(k-1))*μhi, Thi
+            k += 1
+        else
+            μhi, μlo = inv(k+2-z + k*z*μhi), μhi
+            Thi, Tlo = ((k+2-z)*Thi + k*z*Tlo*μlo)*μhi, Thi
+            k += 1
+        end
     end
     k < kmax || @warn "Rational approximation to "*pFq2string(Val(0), Val(0))*" reached the maximum type of ("*string(kmax, ", ", kmax)*")."
     return isfinite(Thi) ? Thi : Tlo

src/weniger.jl

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

test/runtests.jl

@@ -362,6 +362,17 @@ end
             @test pFqdrummond((), (), z) ≈ CS(pFq((), (), CT(z))) atol=atol rtol=rtol
             @test pFqweniger((), (), z) ≈ CS(pFq((), (), CT(z))) atol=atol rtol=rtol
         end
+        @test pFqweniger((), (), S(2)) ≈ S(pFq((), (), T(2))) atol=atol rtol=rtol
+        @test pFqdrummond((), (), S(2)) ≈ S(pFq((), (), T(2))) atol=atol rtol=rtol
+        # real root of z³-6z²+18z-24
+        z = big"2.62581681895846671601188893376528333127944414121175230343105061856653653362177"
+        @test pFqdrummond((), (), S(z)) ≈ S(pFq((), (), T(z))) atol=atol rtol=rtol
+        # real root of z³-12z²+60z-120
+        z = big"4.64437070925217118582294142140806396986328876579916801835768782380859352462557"
+        @test pFqweniger((), (), S(z)) ≈ S(pFq((), (), T(z))) atol=atol rtol=rtol
+        # real root of z⁵-30z⁴+420z³-3360z²+15120z-30240
+        z = big"7.29347719065928651947033927231889084041392903533230450877470222528006165013418"
+        @test pFqweniger((), (), S(z)) ≈ S(pFq((), (), T(z))) atol=atol rtol=rtol
     end
 end