[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Axiom-math] Rational Interpolation
From: |
Martin Rubey |
Subject: |
Re: [Axiom-math] Rational Interpolation |
Date: |
Sat, 13 Aug 2005 21:26:29 +0200 |
In fact, I just implemented yet another algorithm, this time by Beckermann and
Labahn. It seems that this one is even faster. In fact, a *lot* faster.
If you are interested, here comes a preliminary version.
fffg.spad
Description: Binary data
usage:
x: List POLY INT := [i for i in 1..30]; _
y: List POLY INT := [random(10) for i in 1..#x]; _
C:SquareMatrix(#x, POLY INT) := diagonalMatrix x; _
f:Vector SUP POLY INT := [-1, fun]; _
c2:(NNI, Vector SUP POLY INT, Vector SUP POLY INT) -> POLY INT
c2(k, f, M) == eval(dot(f,M)(x.(k+1)),fun=y.(k+1))
r := fffg(C,c2,f,append([1 for i in 1..14],[2 for i in 15..#x]));
computes a certain list of matrix, the interpolant of degree 14/15 then is
r.1.(1,1)/r.1.(2,1)
Another advantage is that this algorithm also computes other recurrence
relations as needed in my guessing package :-)
The corresponding article is available online: I implemented only the Algorithm
from Table 1, the more general algorithm from Table 2 is not yet
done. Furthermore, I did not optimize in any way. In particular, the
representation of the elements of the vector space, the representation of the
special multiplication rule and the computation of c_k are handled in a very
unsatisfactory way. I'm very surprised that the algorithm is even in this
setting so fast. Well, maybe I'm missing something. (Or maybe there is an
efficiency bug in the implementation of Egecioglu-Koc-Gemignanis algorithm...)
"Fraction-free Computation of Matrix Rational Interpolants and Matrix GCDs"
http://math.univ-lille1.fr/~bbecker//ano/pub/1997/ano375.ps
Other publications:
http://math.univ-lille1.fr/~bbecker//bb_pub.html
The following might also be useful:
http://math.univ-lille1.fr/~bbecker//ano/pub/2002/ano445.ps
If you feel like programming, how about implementing a domain for holonomic
(= d-finite = P-recursive) functions?
You might want to read
http://www.risc.uni-linz.ac.at/research/combinat/software/GeneratingFunctions/pub/mallinger96.pdf
to this end.
Martin