Changes http://wiki.axiom-developer.org/FractionFreeFastGaussianElimination/diff
--
??changed:
-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}
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{fffg.spad}
+\author{Martin Rubey}
+\maketitle
+\begin{abstract}
+ The packages defined in this file provide fast fraction free rational
+ interpolation algorithms.
+\end{abstract}
+\tableofcontents
+\section{package FAMR2 FiniteAbelianMonoidRingFunctions2}
+<<package FAMR2 FiniteAbelianMonoidRingFunctions2>>=
+)abbrev package FAMR2 FiniteAbelianMonoidRingFunctions2
+++ Author: Martin Rubey
+++ Description:
+++ This package provides a mapping function for
\spadtype{FiniteAbelianMonoidRing}
+FiniteAbelianMonoidRingFunctions2(E: OrderedAbelianMonoid,
+ R1: Ring,
+ A1: FiniteAbelianMonoidRing(R1, E),
+ R2: Ring,
+ A2: FiniteAbelianMonoidRing(R2, E))
+ : Exports == Implementation where
+ Exports == with
+
+ map: (R1 -> R2, A1) -> A2
+
+ Implementation == add
+
+-- we assume that 0 is mapped onto 0
+ map(f: R1 -> R2, a: A1): A2 ==
+ if zero? a then 0$A2
+ else
+ monomial(f leadingCoefficient a, degree a)$A2 + map(f, reductum a)
+
+@
+
+\section{package FFFG FractionFreeFastGaussian}
+<<package FFFG FractionFreeFastGaussian>>=
)abbrev package FFFG FractionFreeFastGaussian
??changed:
)abbrev package FFFG FractionFreeFastGaussian
-FractionFreeFastGaussian2(D): Exports == Implementation where
- SUP ==> SparseUnivariatePolynomial
-
+++ Author: Martin Rubey
+++ Description:
+++ This package implements the interpolation algorithm proposed in Beckermann,
+++ Bernhard and Labahn, George, Fraction-free computation of matrix rational
+++ interpolants and matrix GCDs, SIAM Journal on Matrix Analysis and
+++ Applications 22.
+FractionFreeFastGaussian(D, V): Exports == Implementation where
D: Join(IntegralDomain, GcdDomain)
??changed:
D: Join(IntegralDomain, GcdDomain)
- V ==> SUP D
-
- cFunction ==> (NonNegativeInteger, Vector V) -> D
+ V: AbelianMonoidRing(D, NonNegativeInteger)
+
+ SUP ==> SparseUnivariatePolynomial
+
+ cFunction ==> (NonNegativeInteger, Vector SUP D) -> D
+
+ CoeffAction ==> (NonNegativeInteger, NonNegativeInteger, V) -> D
+-- coeffAction(k, l, f) is the coefficient of x^k in z^l f(x)
??changed:
- fffg: (List D, cFunction, NonNegativeInteger, List NonNegativeInteger) -> _
- List Matrix V
+ fffg: (List D, cFunction, List NonNegativeInteger) -> Matrix SUP D
+ ++ \spad{fffg} is the general algorithm as proposed by Beckermann and
+ ++ Labahn.
??changed:
- 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
+ interpolate: (List D, List D, NonNegativeInteger) -> Fraction SUP D
+
+ interpolate: (List Fraction D, List Fraction D, NonNegativeInteger)
+ -> Fraction SUP D
+
+ generalInterpolation: (List D, CoeffAction,
+ Vector V, List NonNegativeInteger) -> Matrix SUP D
+ ++ \spad{generalInterpolation(l, CA, f, eta)} performs Hermite-Pade
+ ++ approximation using the given action CA of polynomials on the elements
+ ++ of f. The result is guaranteed to be correct up to order
+ ++ |eta|-1. Given that eta is a "normal" point, the degrees on the
+ ++ diagonal are given by eta. The degrees of column i are in this case
+ ++ eta + e.i - [1,1,...,1], where the degree of zero is -1.
+
+ generalInterpolation: (List D, CoeffAction,
+ Vector V, NonNegativeInteger, NonNegativeInteger)
+ -> Stream Matrix SUP D
+ ++ \spad{generalInterpolation(l, CA, f, totalDegree, maxDegree)} applies
+ ++ \spad{generalInterpolation(l, CA, f, eta)} for all possible \spad{eta}
+ ++ with maximal entry \spad{maxDegree} and sum of entries
+ ++ \spad{totalDegree}.
+
+ generalCoefficient: (CoeffAction, Vector V,
+ NonNegativeInteger, Vector SUP D) -> D
+ ++ \spad{generalCoefficient(action, f, k, p)} gives the coefficient of
+ ++ x^k in p(z)\dot f(x), where the action of z^l on a polynomial in x is
+ ++ given by action.
+
+ ShiftAction: (NonNegativeInteger, NonNegativeInteger, V) -> D
+ ++ \spad{ShiftAction(k, l, g)} gives the coefficient of x^k in z^l g(x),
+ ++ where \spad{z*(a+b*x+c*x^2+d*x^3+...) = (b*x+2*c*x^2+3*d*x^3+...)}. In
+ ++ terms of sequences, z*u(n)=n*u(n).
+
+ ShiftC: NonNegativeInteger -> List D
+ ++ \spad{ShiftC} gives the coefficients c_{k,k} in the expansion <x^k> z
+ ++ g(x) = sum_{i=0}^k c_{k,i} <x^i> g(x), where z acts on g(x) by
+ ++ shifting. In fact, the result is [0,1,2,...]
+
+ DiffAction: (NonNegativeInteger, NonNegativeInteger, V) -> D
+ ++ \spad{DiffAction(k, l, g)} gives the coefficient of x^k in z^l g(x),
+ ++ where z*(a+b*x+c*x^2+d*x^3+...) = (a*x+b*x^2+c*x^3+...).
+
+ DiffC: NonNegativeInteger -> List D
+ ++ \spad{DiffC} gives the coefficients c_{k,k} in the expansion <x^k> z
+ ++ g(x) = sum_{i=0}^k c_{k,i} <x^i> g(x), where z acts on g(x) by
+ ++ shifting. In fact, the result is [0,0,0,...]
+
+ qShiftAction: (D, NonNegativeInteger, NonNegativeInteger, V) -> D
+ ++ \spad{qShiftAction(q, k, l, g)} gives the coefficient of x^k in z^l
+ ++ g(x), where z*(a+b*x+c*x^2+d*x^3+...) =
+ ++ (a+q*b*x+q^2*c*x^2+q^3*d*x^3+...). In terms of sequences,
+ ++ z*u(n)=q^n*u(n).
+
+ qShiftC: (D, NonNegativeInteger) -> List D
+ ++ \spad{qShiftC} gives the coefficients c_{k,k} in the expansion <x^k> z
+ ++ g(x) = sum_{i=0}^k c_{k,i} <x^i> g(x), where z acts on g(x) by
+ ++ shifting. In fact, the result is [1,q,q^2,...]
??changed:
- interpolate(x: List Fraction D, y: List Fraction D, _
- d: NonNegativeInteger) ==
+-------------------------------------------------------------------------------
+-- Shift Operator
+-------------------------------------------------------------------------------
+
+-- ShiftAction(k, l, f) is the CoeffAction appropriate for the shift operator.
+
+ ShiftAction(k: NonNegativeInteger, l: NonNegativeInteger, f: V): D ==
+ k**l*coefficient(f, k)
+
+
+ ShiftC(total: NonNegativeInteger): List D ==
+ [i::D for i in 0..total-1]
+
+-------------------------------------------------------------------------------
+-- q-Shift Operator
+-------------------------------------------------------------------------------
+
+-- q-ShiftAction(k, l, f) is the CoeffAction appropriate for the shift
operator.
+
+ qShiftAction(q: D, k: NonNegativeInteger, l: NonNegativeInteger, f: V): D
==
+ q**(k*l)*coefficient(f, k)
+
+
+ qShiftC(q: D, total: NonNegativeInteger): List D ==
+ [q**i for i in 0..total-1]
+
+-------------------------------------------------------------------------------
+-- Differentiation Operator
+-------------------------------------------------------------------------------
+
+-- DiffAction(k, l, f) is the CoeffAction appropriate for the differentiation
+-- operator.
+
+ DiffAction(k: NonNegativeInteger, l: NonNegativeInteger, f: V): D ==
+ coefficient(f, (k-l)::NonNegativeInteger)
+
+
+ DiffC(total: NonNegativeInteger): List D ==
+ [0 for i in 1..total]
+
+-------------------------------------------------------------------------------
+-- general - suitable for functions f
+-------------------------------------------------------------------------------
+
+-- get the coefficient of z^k in the scalar product of p and f, the action
+-- being defined by coeffAction
+
+ generalCoefficient(coeffAction: CoeffAction, f: Vector V,
+ k: NonNegativeInteger, p: Vector SUP D): D ==
+ res: D := 0
+ for i in 1..#f repeat
+
+-- Defining a and b and summing only over those coefficients that might be
+-- nonzero makes a huge difference in speed
+ a := f.i
+ b := p.i
+ for l in minimumDegree b..degree b repeat
+ if not zero? coefficient(b, l)
+ then res := res + coefficient(b, l) * coeffAction(k, l, a)
+ res
+
+
+ generalInterpolation(C: List D, coeffAction: CoeffAction,
+ f: Vector V,
+ eta: List NonNegativeInteger): Matrix SUP D ==
+
+ c: cFunction := generalCoefficient(coeffAction, f,
+ (#1-1)::NonNegativeInteger, #2)
+ fffg(C, c, eta)
+
+
+
+-- returns the lexicographically next vector with non-negative components
+-- smaller than p with the same sum as v
+ nextVector!(p: NonNegativeInteger, v: List NonNegativeInteger)
+ : Union("failed", List NonNegativeInteger) ==
+ n := #v
+ pos := position(#1 < p, v)
+-- pos should never be zero, since we assume that we have at least one
+-- non-maximal, non-zero element.
+ if pos = 1 then
+ sum: Integer := v.1
+ for i in 2..n repeat
+ if v.i < p and sum > 0 then
+ v.i := v.i + 1
+ sum := sum - 1
+ for j in 1..i-1 repeat
+ if sum > p then
+ v.j := p
+ sum := sum - p
+ else
+ v.j := sum::NonNegativeInteger
+ sum := 0
+ return v
+ else sum := sum + v.i
+ return "failed"
+ else
+ v.pos := v.pos + 1
+ v.(pos-1) := (v.(pos-1) - 1)::NonNegativeInteger
+
+ v
+
+ vectorStream(p: NonNegativeInteger, v: List NonNegativeInteger)
+ : Stream List NonNegativeInteger == delay
+ next := nextVector!(p, copy v)
+ (next case "failed") => empty()$Stream(List NonNegativeInteger)
+ cons(next, vectorStream(p, next))
+
+ vectorStream2(p: NonNegativeInteger, v: List NonNegativeInteger)
+ : Stream List NonNegativeInteger == delay
+ next := nextVector!(p, copy v)
+ (next case "failed") => empty()$Stream(List NonNegativeInteger)
+ next2 := nextVector!(p, copy next)
+ (next2 case "failed") => cons(next, empty())
+ cons(next2, vectorStream2(p, next2))
+
+ generalInterpolation(C: List D, coeffAction: CoeffAction,
+ f: Vector V,
+ totalDegree: NonNegativeInteger,
+ maxDegree: NonNegativeInteger)
+ : Stream Matrix SUP D ==
+ sum: Integer := totalDegree
+ entry: Integer
+ eta: List NonNegativeInteger
+ := [(if sum < maxDegree _
+ then (entry := sum; sum := 0) _
+ else (entry := maxDegree; sum := sum - maxDegree); _
+ entry::NonNegativeInteger) for i in 1..#f]
+
+ (sum > 0) => empty()$Stream(Matrix SUP D)
+
+ if #f = 2 then
+ map(generalInterpolation(C, coeffAction, f, #1),
+ cons(eta, vectorStream2(maxDegree, eta)))
+ $StreamFunctions2(List NonNegativeInteger,
+ Matrix SUP D)
+ else
+ map(generalInterpolation(C, coeffAction, f, #1),
+ cons(eta, vectorStream(maxDegree, eta)))
+ $StreamFunctions2(List NonNegativeInteger,
+ Matrix SUP D)
+
+-------------------------------------------------------------------------------
+-- rational interpolation
+-------------------------------------------------------------------------------
+
+ interpolate(x: List Fraction D, y: List Fraction D, d: NonNegativeInteger)
+ : Fraction SUP D ==
gx := splitDenominator(x)$InnerCommonDenominator(D, Fraction D, _
??changed:
- interpolate(x: List D, y: List D, d: NonNegativeInteger) ==
+ interpolate(x: List D, y: List D, d: NonNegativeInteger): Fraction SUP D ==
-- berechne Interpolante mit Graden d und N-d-1
??changed:
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)
+ M := fffg(x, c, eta)
+
+ if zero?(M.(2,1)) then M.(1,2)/M.(2,2)
+ else M.(1,1)/M.(2,1)
??changed:
-
--- 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!
-[3 more lines...]
+-------------------------------------------------------------------------------
+-- fffg
+-------------------------------------------------------------------------------
+
+-- a major part of the time is spent here
+ recurrence(M: Matrix SUP D, lambda: NonNegativeInteger, m:
NonNegativeInteger,
+ r: Vector D, d: D, z: SUP D, Ck: D, p: Vector D): Matrix SUP D
==
+
+ rLambda := qelt(r, lambda)
+ polyf := rLambda * (z - Ck::SUP D)
+
+ for i in 1..m repeat
+ Milambda := qelt(M, i, lambda)
+ newMilambda := polyf * Milambda
+
+
+-- update columns ~= lambda and calculate their sum
+ for l in 1..m | l ~= lambda repeat
+ rl := qelt(r, l)
+ Mil := ((qelt(M, i, l) * rLambda - Milambda * rl) exquo d)::SUP D
+ qsetelt!(M, i, l, Mil)
+
+ pl := qelt(p, l)
+ newMilambda := newMilambda - pl * Mil
+
+-- update column lambda
+
+ qsetelt!(M, i, lambda, (newMilambda exquo d)::SUP D)
+
+ M
+
+
+ fffg(C: List D, c: cFunction, eta: List NonNegativeInteger): Matrix SUP D
==
+
+-- eta is the vector of degrees. We compute M with degrees eta+e_i-1, i=1..m
+ z: SUP D := monomial(1, 1)
+ m: NonNegativeInteger := #eta
+ M: Matrix SUP D := scalarMatrix(m, 1)
+ d: D := 1
K: NonNegativeInteger := reduce(_+, eta)
--removed:
p: Vector D := zero(m)
- Mold: Matrix V
- Mnew: Matrix V
- Col: Vector V := zero(m)
Lambda: List Integer
??changed:
-- k = sigma+1
- Mold := M.first
-
- for l in 1..m repeat
- r.l := c(k, column(Mold, l))
+
+ for l in 1..m repeat r.l := c(k, column(M, l))
??changed:
- 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
-[18 more lines...]
+-- if Lambda is empty, then M, d and etak remain unchanged. Otherwise, we look
+-- for the next closest para-normal point.
+
+ (empty? Lambda) => "iterate"
+
+ lambdaMax := reduce(max, Lambda)
+ lambda := 1
+ while eta.lambda-etak.lambda < lambdaMax or r.lambda = 0 repeat
+ lambda := lambda + 1
+
+-- Calculate leading coefficients
+
+ for l in 1..m | l ~= lambda repeat
+ if etak.l > 0 then
+ p.l := coefficient(M.(l, lambda), (etak.l-1)::NonNegativeInteger)
+ else
+ p.l := 0
+
+-- increase order and adjust degree constraints
+
+ M := recurrence(M, lambda, m, r, d, z, C.k, p)
+
+ d := r.lambda
+ etak.lambda := etak.lambda + 1
??changed:
M
-\end{spad}
+
+@
+%$
+
+\section{package FFFGF FractionFreeFastGaussianFractions}
+<<package FFFGF FractionFreeFastGaussianFractions>>=
+)abbrev package FFFGF FractionFreeFastGaussianFractions
+++ Author: Martin Rubey
+++ Description:
+++ This package lifts the interpolation functions from
+++ \spadtype{FractionFreeFastGaussian} to fractions.
+FractionFreeFastGaussianFractions(D, V, VF): Exports == Implementation where
+ D: Join(IntegralDomain, GcdDomain)
+ V: FiniteAbelianMonoidRing(D, NonNegativeInteger)
+ VF: FiniteAbelianMonoidRing(Fraction D, NonNegativeInteger)
+
+ F ==> Fraction D
+
+ SUP ==> SparseUnivariatePolynomial
+
+ FFFG ==> FractionFreeFastGaussian
+
+ FAMR2 ==> FiniteAbelianMonoidRingFunctions2
+
+ cFunction ==> (NonNegativeInteger, Vector SUP D) -> D
+
+ CoeffAction ==> (NonNegativeInteger, NonNegativeInteger, V) -> D
+-- coeffAction(k, l, f) is the coefficient of x^k in z^l f(x)
+
+ Exports == with
+
+ generalInterpolation: (List D, CoeffAction, Vector VF, List
NonNegativeInteger)
+ -> Matrix SUP D
+ ++ \spad{generalInterpolation(l, CA, f, eta)} performs Hermite-Pade
+ ++ approximation using the given action CA of polynomials on the elements
+ ++ of f. The result is guaranteed to be correct up to order
+ ++ |eta|-1. Given that eta is a "normal" point, the degrees on the
+ ++ diagonal are given by eta. The degrees of column i are in this case
+ ++ eta + e.i - [1,1,...,1], where the degree of zero is -1.
+
+ generalInterpolation: (List D, CoeffAction,
+ Vector VF, NonNegativeInteger, NonNegativeInteger)
+ -> Stream Matrix SUP D
+ ++ \spad{generalInterpolation(l, CA, f, totalDegree, maxDegree)} applies
+ ++ generalInterpolation(l, CA, f, eta) for all possible eta with maximal
+ ++ entry maxDegree and sum of entries totalDegree
+
+ Implementation == add
+
+ multiplyRows!(v: Vector D, M: Matrix SUP D): Matrix SUP D ==
+ n := #v
+ for i in 1..n repeat
+ for j in 1..n repeat
+ M.(i,j) := v.i*M.(i,j)
+
+ M
+
+ generalInterpolation(C: List D, coeffAction: CoeffAction,
+ f: Vector VF, eta: List NonNegativeInteger): Matrix
SUP D ==
+ n := #f
+ g: Vector V := new(n, 0)
+ den: Vector D := new(n, 0)
+
+ for i in 1..n repeat
+ c := coefficients(f.i)
+ den.i := commonDenominator(c)$CommonDenominator(D, F, List F)
+ g.i := map(retract(#1*den.i)@D, f.i)
+ $FAMR2(NonNegativeInteger, Fraction D, VF, D, V)
+
+ M := generalInterpolation(C, coeffAction, g, eta)$FFFG(D, V)
+
+-- The following is necessary since I'm multiplying each row with a factor, not
+-- each column. Possibly I could factor out gcd den, but I'm not sure whether
+-- this is efficient.
+
+ multiplyRows!(den, M)
+
+ generalInterpolation(C: List D, coeffAction: CoeffAction,
+ f: Vector VF,
+ totalDegree: NonNegativeInteger,
+ maxDegree: NonNegativeInteger)
+ : Stream Matrix SUP D ==
+
+ n := #f
+ g: Vector V := new(n, 0)
+ den: Vector D := new(n, 0)
+
+ for i in 1..n repeat
+ c := coefficients(f.i)
+ den.i := commonDenominator(c)$CommonDenominator(D, F, List F)
+ g.i := map(retract(#1*den.i)@D, f.i)
+ $FAMR2(NonNegativeInteger, Fraction D, VF, D, V)
+
+ c: cFunction := generalCoefficient(coeffAction, g,
+ (#1-1)::NonNegativeInteger,
#2)$FFFG(D, V)
+
+
+ MS: Stream Matrix SUP D := generalInterpolation(C, coeffAction, g,
totalDegree, maxDegree)$FFFG(D, V)
+
+-- The following is necessary since I'm multiplying each row with a factor, not
+-- each column. Possibly I could factor out gcd den, but I'm not sure whether
+-- this is efficient.
+
+ map(multiplyRows!(den, #1), MS)$Stream(Matrix SUP D)
+
+@
+
+
+\section{package NEWTON NewtonInterpolation}
+<<package NEWTON NewtonInterpolation>>=
+)abbrev package NEWTON NewtonInterpolation
+++ Description:
+++ This package exports Newton interpolation for the special case where the
+++ result is known to be in the original integral domain
+NewtonInterpolation F: Exports == Implementation where
+ F: IntegralDomain
+ Exports == with
+ newton: List F -> SparseUnivariatePolynomial F
+
+ Implementation == add
+
+ differences(yl: List F): List F ==
+ [y2-y1 for y1 in yl for y2 in rest yl]
+
+ z:SparseUnivariatePolynomial(F) := monomial(1,1)
+
+-- we assume x=[1,2,3,...,n]
+ newtonAux(k: F, fact: F, yl: List F): SparseUnivariatePolynomial(F) ==
+ if empty? rest yl
+ then ((yl.1) exquo fact)::F::SparseUnivariatePolynomial(F)
+ else ((yl.1) exquo fact)::F::SparseUnivariatePolynomial(F)
+ + (z-k::SparseUnivariatePolynomial(F)) _
+ * newtonAux(k+1$F, fact*k, differences yl)
+
+
+ newton yl == newtonAux(1$F, 1$F, yl)
+
+@
+%$
+
+\section{License}
+<<license>>=
+-- Copyright (C) 2006 Martin Rubey <[EMAIL PROTECTED]>
+--
+-- This program is free software; you can redistribute it and/or
+-- modify it under the terms of the GNU General Public License as
+-- published by the Free Software Foundation; either version 2 of
+-- the License, or (at your option) any later version.
+--
+-- This program is distributed in the hope that it will be
+-- useful, but WITHOUT ANY WARRANTY; without even the implied
+-- warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+-- PURPOSE. See the GNU General Public License for more details.
+
+@
+
+<<*>>=
+<<license>>
+
+<<package FAMR2 FiniteAbelianMonoidRingFunctions2>>
+<<package FFFG FractionFreeFastGaussian>>
+<<package FFFGF FractionFreeFastGaussianFractions>>
+<<package NEWTON NewtonInterpolation>>
+@
+\end{document}
+
--
forwarded from http://wiki.axiom-developer.org/[EMAIL PROTECTED]
