diff options
Diffstat (limited to 'src/algebra/pade.spad.pamphlet')
| -rw-r--r-- | src/algebra/pade.spad.pamphlet | 247 |
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} |