\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} <>= )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} <>= )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*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*x**n + a*x**n+1 + ...)/b -- = x**n/b + (a + a*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} <>= --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. @ <<*>>= <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}