Changes http://wiki.axiom-developer.org/FffgSpad/diff
--
\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
++ 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)
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)
Exports == with
fffg: (List D, cFunction, List NonNegativeInteger) -> Matrix SUP D
++ \spad{fffg} is the general algorithm as proposed by Beckermann and
++ Labahn.
-- 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 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,...]
Implementation ==> add
-------------------------------------------------------------------------------
-- 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, _
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): Fraction SUP D ==
-- 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, eta)
if zero?(M.(2,1)) then M.(1,2)/M.(2,2)
else M.(1,1)/M.(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.
-------------------------------------------------------------------------------
-- 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)
etak: Vector NonNegativeInteger := zero(m)
r: Vector D := zero(m)
p: Vector D := zero(m)
Lambda: List Integer
lambdaMax: Integer
lambda: NonNegativeInteger
for k in 1..K repeat
-- k = sigma+1
for l in 1..m repeat r.l := c(k, column(M, l))
Lambda := [eta.l-etak.l for l in 1..m | r.l ~= 0]
-- 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
M
@
%$
\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]
