Merge pull request #16 from MikaelSlevinsky/master

Basic support for Pade approximants

Author: Keno

src/Polynomials.jl

@@ -4,6 +4,7 @@ module Polynomials
 #todo: sparse polynomials?
 
 export Poly, polyval, polyint, polyder, poly, roots
+export Pade, padeval
 
 import Base: length, endof, getindex, setindex!, copy, zero, one, convert
 import Base: show, print, *, /, //, -, +, ==, divrem, rem, eltype
@@ -341,7 +342,7 @@ function roots{T}(p::Poly{T})
     return r
 end
 
-function gcd{T<:FloatingPoint, S<:FloatingPoint}(a::Poly{T}, b::Poly{S})
+function gcd{T, S}(a::Poly{T}, b::Poly{S})
     #Finds the Greatest Common Denominator of two polynomials recursively using
     #Euclid's algorithm: http://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor#Euclid.27s_algorithm
     if all(abs(b.a).<=2*eps(S))
@@ -351,4 +352,7 @@ function gcd{T<:FloatingPoint, S<:FloatingPoint}(a::Poly{T}, b::Poly{S})
         return gcd(b, r)
     end
 end
+
+include("pade.jl")
+
 end # module Poly

src/pade.jl

@@ -0,0 +1,33 @@
+immutable Pade{T<:Number,S<:Number}
+    p::Poly{T}
+    q::Poly{S}
+    var::Symbol
+    function Pade(p::Poly{T},q::Poly{S})
+        if p.var != q.var
+            error("Polynomials must have same variable")
+        end
+        new(p, q, p.var)
+    end
+end
+Pade{T<:Number,S<:Number}(p::Poly{T}, q::Poly{S}) = Pade{T,S}(p,q)
+
+function Pade{T}(c::Poly{T},m::Int,n::Int)
+    @assert m+n < length(c)
+    rold,rnew = Poly([zeros(T,m+n+1),one(T)],c.var),Poly(c.a[1:m+n+1],c.var)
+    uold,vold = Poly([one(T)],c.var),Poly([zero(T)],c.var)
+    unew,vnew = vold,uold
+    for i=1:n
+        temp0,temp1,temp2 = rnew,unew,vnew
+        q,rnew = divrem(rold,rnew)
+        unew,vnew = uold-q*unew,vold-q*vnew
+        rold,uold,vold = temp0,temp1,temp2
+
+    end
+    if vnew[0] == 0
+        d = gcd(rnew,vnew)
+        rnew /= d
+        vnew /= d
+    end
+    Pade(rnew/vnew[0],vnew/vnew[0])
+end
+padeval{T}(PQ::Pade{T},x) = polyval(PQ.p,x)./polyval(PQ.q,x)

test/tests.jl

@@ -86,3 +86,30 @@ pS3 = Poly([1, 2, 3, 4, 5], :s)
 @test_throws ErrorException pS1 * pX
 @test_throws ErrorException pS1 / pX
 @test_throws ErrorException pS1 % pX
+
+#Tests for Pade approximants
+
+println("Test for the exponential function.")
+a = Poly(1.//convert(Vector{Int},gamma(1:17)),"x")
+PQexp = Pade(a,8,8)
+@test padeval(PQexp,1.0) == exp(1.0)
+@test padeval(PQexp,-1.0) == exp(-1.0)
+
+println("Test for the sine function.")
+b = Poly(convert(Vector{BigInt},sinpi((0:16)/2)).//convert(Vector{BigInt},gamma(BigFloat("1.0"):BigFloat("17.0"))),"x")
+PQsin = Pade(b,8,7)
+@test isapprox(padeval(PQsin,1.0),sin(1.0))
+@test isapprox(padeval(PQsin,-1.0),sin(-1.0))
+
+println("Test for the cosine function.")
+c = Poly(convert(Vector{BigInt},sinpi((1:17)/2)).//convert(Vector{BigInt},gamma(BigFloat("1.0"):BigFloat("17.0"))),"x")
+PQcos = Pade(c,8,8)
+@test isapprox(padeval(PQcos,1.0),cos(1.0))
+@test isapprox(padeval(PQcos,-1.0),cos(-1.0))
+
+println("Test for the summation of a factorially divergent series.")
+d = Poly(convert(Vector{BigInt},(-1).^(0:60).*gamma(BigFloat("1.0"):BigFloat("61.0"))).//1,"x")
+PQexpint = Pade(d,30,30)
+println("The approximate sum of the divergent series is:  ",float64(padeval(PQexpint,1.0)))
+println("The approximate sum of the convergent series is: ",exp(1)*(-γ-sum([(-1).^k/k./gamma(k+1) for k=1:20])))
+@test isapprox(padeval(PQexpint,1.0) , exp(1)*(-γ-sum([(-1).^k/k./gamma(k+1) for k=1:20])))