[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Axiom-math] Writing packages
From: |
Fabio S. |
Subject: |
[Axiom-math] Writing packages |
Date: |
Fri, 12 May 2006 18:31:36 +0200 (CEST) |
Hi all,
I am trying to start to write packages, but the package I wrote does not
compile. Yet I started from a function which compiles cleanly and works as
expected. I include at the end both the function and the package.
Now some detail. To start, I have chosen to write a function which
computes the minimal polynomial of a matrix (if I am not wrong, there is
not such a function in the library, isn't it?)
The function that I wrote seems to work and compiles. The algorithm I used
is quite dumb: if m is an nxn matrix, then it tries to solve the linear
systems
m^i+x_1*m^(i-1)*...*x_i=0 for i=1,..,n
The minimal i for which this system is solvable is the degree of the
minimal polynomial and the x_i are its coefficients.
The fact that there must be a solution for i<=n follows from the
Hamilton-Cayley theorem, of course.
As I said, the following function works and compiles. The following
package, on the other hand, doesn't compile. Yet I just add what is
necessary to convert the function into a package (I have read thoroughly
chapter 11 of the book and also I looked at other packages in the source
tree, in particular eigen.spad).
Any help?
BTW, I also have another question on this topic: if I compute the
characteristic polynomial of a matrix and then I evaluate the polynomial
at the matrix, I would expect to have the zero matrix as result. This
does not happens: try
m:=matrix([[1,2],[3,4]])
p:=characteristicPolynomial(m)
p(m)
Why?
Thanks
Fabio
---------------------- This is the function
R ==> INT
F ==> Fraction R
P ==> Polynomial F
M ==> Matrix(F)
NNI ==> NonNegativeInteger
minpoly: (M,Symbol) -> P
minpoly(m,x) ==
dimm : NNI := nrows m ++ dimension of m
dimm ~= ncols m => error "The matrix is not square"
minpol : P := 0
dimm = 1 => return(minpol : P := x - m(1,1))
powerm : Matrix INT := new(dimm,dimm,0) ++ powers of m
for i in 1..dimm repeat powerm(i,i) := 1 ++ firts power of m:
m^0 = I
eqtns : List Vector F := [[] for i in 1..dimm*dimm] ++ Matrix of
coefficients
++ for i in 1..dimm solve the linear system
++ m^i+x_1*m^(i-1)*...*x_i=0
++ or, equivalently,
++ x_1*m^(i-1)*...*x_i=-m^i
for i in 1..dimm repeat
for j in 1..dimm*dimm repeat
ji := (j-1) quo dimm + 1
jj := (j-1) rem dimm + 1
eqtns(j) := append(eqtns(j),[powerm(ji,jj)])
powerm := powerm*m
sols : List F := [] ++ Vector of solutions
for j in 1..dimm repeat sols := append(sols,-1*row(powerm,j))
coeffs := solve(eqtns,sols)
if coeffs.particular case "failed" then iterate
++ The first time the system has a solution, it means that
++ we have a (monic) polynomial of least degree which is
++ 0 on m : then leave the loop i and reconstruct the polynomial
++ from the solution of the system.
break
for j in 1..i repeat
minpol := minpol+(coeffs.particular)(j)*x^((j-1)::NNI)
minpol := minpol+x^i
------------------- This is the package
)abbrev package MMPOL MatrixMinimalPolynomialPackage
MatrixMinimalPolynomialPackage(R) : Exports == Implementation where
R : IntegralDomain
F ==> Fraction R
P ==> Polynomial F
M ==> Matrix(F)
NNI ==> NonNegativeInteger
Exports == with
minpoly : (M,Symbol) -> P
++ minimalPolynomial(m,var) returns the minimal polynomial
++ of the matrix m using the symbol var as the main variabe.
Implementation == add
minpoly: (M,Symbol) -> P
minpoly(m,x) ==
dimm : NNI := nrows m ++ dimension of m
dimm ~= ncols m => error "The matrix is not square"
minpol : P := 0
dimm = 1 => return(minpol : P := x - m(1,1))
powerm : Matrix INT := new(dimm,dimm,0) ++ powers of m
for i in 1..dimm repeat powerm(i,i) := 1 ++ firt power of m:
m^0 = I
eqtns : List Vector F := [[] for i in 1..dimm*dimm] ++ Matrix of
coefficients
++ for i in 1..dimm solve the linear system
++ m^i+x_1*m^(i-1)*...*x_i=0
++ or, equivalently,
++ x_1*m^(i-1)*...*x_i=-m^i
for i in 1..dimm repeat
for j in 1..dimm*dimm repeat
ji := (j-1) quo dimm + 1
jj := (j-1) rem dimm + 1
eqtns(j) := append(eqtns(j),[powerm(ji,jj)])
powerm := powerm*m
sols : List F := [] ++ Vector of
solutions
for j in 1..dimm repeat sols := append(sols,-1*row(powerm,j))
coeffs := solve(eqtns,sols)
if coeffs.particular case "failed" then iterate
++ The first time the system has a solution, it means that
++ we have a (monic) polynomial of least degree which is
++ 0 on m : then leave the loop i and reconstruct the polynomial
++ from the solution of the system.
break
for j in 1..i repeat
minpol := minpol+(coeffs.particular)(j)*x^((j-1)::NNI)
minpol := minpol+x^i
- [Axiom-math] Writing packages,
Fabio S. <=