path: root/src/algebra/pade.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra pade.spad}
\author{Barry Trager, William Burge, Martin Hassner,Stephen M. Watt}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package PADEPAC PadeApproximantPackage}
<<package PADEPAC PadeApproximantPackage>>=
)abbrev package PADEPAC PadeApproximantPackage
++ This package computes reliable Pad&ea. approximants using
++ a generalized Viskovatov continued fraction algorithm.
++ Authors: Trager,Burge, Hassner & Watt.
++ Date Created: April 1987
++ Date Last Updated: 12 April 1990
++ Keywords: Pade, series
++ Examples:
++ References:
++   "Pade Approximants, Part I: Basic Theory", Baker & Graves-Morris.

PadeApproximantPackage(R: Field, x:Symbol, pt:R): Exports == Implementation where
   PS ==> UnivariateTaylorSeries(R,x,pt)
   UP ==> UnivariatePolynomial(x,R)
   QF ==> Fraction UP
   CF  ==> ContinuedFraction UP
   NNI ==> NonNegativeInteger
   Exports ==> with
     pade:   (NNI,NNI,PS,PS) -> Union(QF,"failed")
      ++ pade(nd,dd,ns,ds) computes the approximant as a quotient of polynomials
      ++ (if it exists) for arguments
      ++ nd (numerator degree of approximant),
      ++ dd (denominator degree of approximant),
      ++ ns (numerator series of function), and
      ++ ds (denominator series of function).
     pade:   (NNI,NNI,PS) -> Union(QF,"failed")
      ++ pade(nd,dd,s)
      ++ computes the quotient of polynomials
      ++ (if it exists) with numerator degree at
      ++ most nd and denominator degree at most dd
      ++ which matches the series s to order \spad{nd + dd}.

   Implementation ==> add
     n,m : NNI
     u,v : PS
     pa := PadeApproximants(R,PS,UP)
     pade(n,m,u,v) ==
       ans:=pade(n,m,u,v)$pa
       ans case "failed" => ans
       pt = 0 => ans
       num := numer(ans::QF)
       den := denom(ans::QF)
       xpt : UP := monomial(1,1)-monomial(pt,0)
       num := num(xpt)
       den := den(xpt)
       num/den
     pade(n,m,u) == pade(n,m,u,1)

@
\section{package PADE PadeApproximants}
<<package PADE PadeApproximants>>=
)abbrev package PADE PadeApproximants
++ This package computes reliable Pad&ea. approximants using
++ a generalized Viskovatov continued fraction algorithm.
++ Authors: Burge, Hassner & Watt.
++ Date Created: April 1987
++ Date Last Updated: 12 April 1990
++ Keywords: Pade, series
++ Examples:
++ References:
++   "Pade Approximants, Part I: Basic Theory", Baker & Graves-Morris.
PadeApproximants(R,PS,UP): Exports == Implementation where
  R:  Field -- IntegralDomain
  PS: UnivariateTaylorSeriesCategory R
  UP: UnivariatePolynomialCategory R
 
  NNI ==> NonNegativeInteger
  QF  ==> Fraction UP
  CF  ==> ContinuedFraction UP
 
  Exports ==> with
    pade:   (NNI,NNI,PS,PS) -> Union(QF,"failed")
      ++ pade(nd,dd,ns,ds)
      ++ computes the approximant as a quotient of polynomials
      ++ (if it exists) for arguments
      ++ nd (numerator degree of approximant),
      ++ dd (denominator degree of approximant),
      ++ ns (numerator series of function), and
      ++ ds (denominator series of function).
    padecf: (NNI,NNI,PS,PS) -> Union(CF, "failed")
      ++ padecf(nd,dd,ns,ds)
      ++ computes the approximant as a continued fraction of
      ++ polynomials (if it exists) for arguments
      ++ nd (numerator degree of approximant),
      ++ dd (denominator degree of approximant),
      ++ ns (numerator series of function), and
      ++ ds (denominator series of function).
 
  Implementation ==> add
    -- The approximant is represented as
    --   p0 + x**a1/(p1 + x**a2/(...))
 
    PadeRep ==> Record(ais: List UP, degs: List NNI) -- #ais= #degs
    PadeU   ==> Union(PadeRep, "failed")             -- #ais= #degs+1
 
    constInner(up:UP):PadeU == [[up], []]
 
    truncPoly(p:UP,n:NNI):UP ==
      while n < degree p repeat p := reductum p
      p
 
    truncSeries(s:PS,n:NNI):UP ==
      p: UP := 0
      for i in 0..n repeat p := p + monomial(coefficient(s,i),i)
      p
 
    -- Assumes s starts with a<n>*x**n + ... and divides out x**n.
    divOutDegree(s:PS,n:NNI):PS ==
      for i in 1..n repeat s := quoByVar s
      s
 
    padeNormalize: (NNI,NNI,PS,PS) -> PadeU
    padeInner:     (NNI,NNI,PS,PS) -> PadeU
 
    pade(l,m,gps,dps) ==
      (ad := padeNormalize(l,m,gps,dps)) case "failed" => "failed"
      plist := ad.ais; dlist := ad.degs
      approx := first(plist) :: QF
      for d in dlist for p in rest plist repeat
        approx := p::QF + (monomial(1,d)$UP :: QF)/approx
      approx
 
    padecf(l,m,gps,dps) ==
      (ad := padeNormalize(l,m,gps,dps)) case "failed" => "failed"
      alist := reverse(ad.ais)
      blist := [monomial(1,d)$UP for d in reverse ad.degs]
      continuedFraction(first(alist),_
                          blist::Stream UP,(rest alist) :: Stream UP)
 
    padeNormalize(l,m,gps,dps) ==
      zero? dps => "failed"
      zero? gps => constInner 0
      -- Normalize so numerator or denominator has constant term.
      ldeg:= min(order dps,order gps)
      if ldeg > 0 then
        dps := divOutDegree(dps,ldeg)
        gps := divOutDegree(gps,ldeg)
      padeInner(l,m,gps,dps)
 
    padeInner(l, m, gps, dps) ==
      zero? coefficient(gps,0) and zero? coefficient(dps,0) =>
        error "Pade' problem not normalized."
      plist: List UP := nil()
      alist: List NNI := nil()
      -- Ensure denom has constant term.
      if zero? coefficient(dps,0) then
        -- g/d = 0 + z**0/(d/g)
        (gps,dps) := (dps,gps)
        (l,m)     := (m,l)
        plist := concat(0,plist)
        alist := concat(0,alist)
      -- Ensure l >= m, maintaining coef(dps,0)~=0.
      if l < m then
        --   (a<n>*x**n + a<n+1>*x**n+1 + ...)/b
        -- = x**n/b + (a<n> + a<n+1>*x + ...)/b
        alpha := order gps
        if alpha > l then return "failed"
        gps := divOutDegree(gps, alpha)
        (l,m) := (m,(l-alpha) :: NNI)
        (gps,dps) := (dps,gps)
        plist := concat(0,plist)
        alist := concat(alpha,alist)
      degbd: NNI := l + m + 1
      g := truncSeries(gps,degbd)
      d := truncSeries(dps,degbd)
      for j in 0.. repeat
        -- Normalize d so constant coefs cancel. (B&G-M is wrong)
        d0 := coefficient(d,0)
        d := (1/d0) * d; g := (1/d0) * g
        p : UP := 0; s := g
        if negative?(l-m+1) then error "Internal pade error"
        degbd := (l-m+1) :: NNI
        for k in 1..degbd repeat
          pk := coefficient(s,0)
          p  := p + monomial(pk,(k-1) :: NNI)
          s  := s - pk*d
          s  := (s exquo monomial(1,1)) :: UP
        plist := concat(p,plist)
        s = 0 => return [plist,alist]
        alpha := minimumDegree(s) + degbd
        alpha > l + m => return [plist,alist]
        alpha > l     => return "failed"
        alist := concat(alpha,alist)
        h := (s exquo monomial(1,minimumDegree s)) :: UP
        degbd := (l + m - alpha) :: NNI
        g := truncPoly(d,degbd)
        d := truncPoly(h,degbd)
        (l,m) := (m,(l-alpha) :: NNI)

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>

<<package PADEPAC PadeApproximantPackage>>
<<package PADE PadeApproximants>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}