Writing packages

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