Initial population.

1 files changed, 247 insertions, 0 deletions

diff --git a/src/algebra/pade.spad.pamphlet b/src/algebra/pade.spad.pamphlet
new file mode 100644
index 00000000..de1a71b2
--- /dev/null
+++ b/src/algebra/pade.spad.pamphlet
@@ -0,0 +1,247 @@
+\documentclass{article}
+\usepackage{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 l-m+1 < 0 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}