/usr/share/axiom-20170501/src/algebra/PADE.spad is in axiom-source 20170501-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 | )abbrev package PADE PadeApproximants
++ Authors: Burge, Hassner & Watt.
++ Date Created: April 1987
++ Date Last Updated: 12 April 1990
++ References:
++ George A. Baker and Peter Graves-Morris
++ ``Pade Approximants''
++ Cambridge University Press, March 1996 ISBN 9870521450072
++ Description:
++ This package computes reliable Pad&ea. approximants using
++ a generalized Viskovatov continued fraction algorithm.
PadeApproximants(R,PS,UP) : SIG == CODE where
R : Field -- IntegralDomain
PS : UnivariateTaylorSeriesCategory R
UP : UnivariatePolynomialCategory R
NNI ==> NonNegativeInteger
QF ==> Fraction UP
CF ==> ContinuedFraction UP
SIG ==> 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).
CODE ==> 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)
|