Changes http://wiki.axiom-developer.org/FractionFreeFastGaussianElimination/diff
--
This is an implementation of the algorithm of Bernhard Beckermann and George
Labahn as described in::
Beckermann, Bernhard(F-LILL-LA); Labahn, George(3-WTRL-C)
Fraction-free computation of matrix rational interpolants and matrix GCDs.
(English. English summary)
SIAM J. Matrix Anal. Appl. 22 (2000), no. 1, 114--144 (electronic).
\begin{spad}
)abbrev package FFFG2 FractionFreeFastGaussian2
FractionFreeFastGaussian2(D): Exports == Implementation where
SUP ==> SparseUnivariatePolynomial
D: Join(IntegralDomain, GcdDomain)
V ==> SUP D
cFunction ==> (NonNegativeInteger, Vector V) -> D
Exports == with
fffg: (List D, cFunction, NonNegativeInteger, List NonNegativeInteger) -> _
List Matrix V
-- c(k, M) computes $c_k(M)$ as in Equation (2). Note that the information
-- about $f$ is encoded in $c$.
interpolate: (List D, List D, NonNegativeInteger) -> Fraction V
interpolate: (List Fraction D, List Fraction D, NonNegativeInteger) -> _
Fraction V
HermitePade: (Vector V, List NonNegativeInteger) -> Matrix V
Implementation ==> add
interpolate(x: List Fraction D, y: List Fraction D, _
d: NonNegativeInteger) ==
gx := splitDenominator(x)$InnerCommonDenominator(D, Fraction D, _
List D, _
List Fraction D)
gy := splitDenominator(y)$InnerCommonDenominator(D, Fraction D, _
List D, _
List Fraction D)
r := interpolate(gx.num, gy.num, d)
elt(numer r, monomial(gx.den,1))/(gy.den*elt(denom r, monomial(gx.den,1)))
interpolate(x: List D, y: List D, d: NonNegativeInteger) ==
-- berechne Interpolante mit Graden d und N-d-1
if (N := #x) ~= #y then
error "interpolate: number of points and values must match"
if N <= d then
error "interpolate: numerator degree must be smaller than number of
data points"
c: cFunction := y.#1 * elt(#2.2, x.#1) - elt(#2.1, x.#1)
eta: List NonNegativeInteger := [d, (N-d)::NonNegativeInteger]
M := fffg(x, c, 2, eta)
if zero?(M.1.(2,1)) then M.1.(1,2)/M.1.(2,2)
else M.1.(1,1)/M.1.(2,1)
-- wegen Lemma 5.3 können M.1.(2,1) und M.1.(2,2) nicht beide verschwinden,
-- denn eta_sigma ist immer sigma-normal (Theorem 7.2) und damit auch
-- paranormal (Dfn 4.2)
-- wegen Lemma 5.1 ist M.1.(*,2) eine Lösung des Interpolationsproblems, wenn
-- M.1.(2,1) verschwindet.
-- get the coefficient of z^k in the scalar product
coeff(k: NonNegativeInteger, f: Vector V, g: Vector V): D ==
res: D := 0
for i in 1..#f repeat
for l in 0..k repeat
res := res + coefficient(f.i, l) _
*coefficient(g.i, (k-l)::NonNegativeInteger)
res
HermitePade(f, eta) ==
c: cFunction := coeff((#1-1)::NonNegativeInteger, f, #2)
C: List D := [0 for i in 1..reduce(_+, eta)]
fffg(C, c, #f, eta).1
fffg(C, c, m, eta) ==
-- eta ist der Gradvektor!
-- berechne M mit Graden eta+e_i-1, i=1..m
z: V := monomial(1,1)
M: List Matrix V := [scalarMatrix(m,1)]
d: List D := [1]
K: NonNegativeInteger := reduce(_+, eta)
etak: Vector NonNegativeInteger := zero(m)
r: Vector D := zero(m)
p: Vector D := zero(m)
Mold: Matrix V
Mnew: Matrix V
Col: Vector V := zero(m)
Lambda: List Integer
lambdaMax: Integer
lambda: NonNegativeInteger
for k in 1..K repeat
-- k = sigma+1
Mold := M.first
for l in 1..m repeat
r.l := c(k, column(Mold, l))
Lambda := [eta.l-etak.l for l in 1..m | r.l ~= 0]
if empty?(Lambda) then
output("Lambda leer")$OutputPackage
M := cons(M.first, M)
d := cons(d.first, d)
-- etak bleibt gleich
else
lambdaMax := reduce(max, Lambda)
lambda := 1
while eta.lambda-etak.lambda < lambdaMax or r.lambda = 0 repeat
lambda := lambda + 1
for l in 1..m | l ~= lambda repeat
if etak.l > 0 then
p.l := coefficient(Mold.(l, lambda), _
(etak.l-1)::NonNegativeInteger)
else
p.l := 0
Mnew := new(m, m, 0)
for l in 1..m | l ~= lambda repeat
setColumn!(Mnew, l, map(((#1*r.lambda-#2*r.l) exquo d.first)::V, _
column(Mold, l), column(Mold, lambda)))
Col := (r.lambda * (z-(C.k)::V)) * column(Mold, lambda)
for l in 1..m | l ~= lambda repeat
Col := Col - (p.l::V) * column(Mnew, l)
Col := map((#1 exquo d.first)::V, Col)
setColumn!(Mnew, lambda, Col)
M := cons(Mnew, M)
d := cons(r.lambda, d)
etak.lambda := etak.lambda + 1
M
\end{spad}
--
forwarded from http://wiki.axiom-developer.org/[EMAIL PROTECTED]
