[Axiom-math] Writing packages

[Fabio S.]

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