Re: [Axiom-math] Writing packages

[Ralf Hemmecke]

Hello Fabio,

Let me be the bad guy... ;-)
(And what I write here is not meant personal, but more sort of a general 
remark.)

If you start writing a package you should immediately stop thinking in 
SPAD and use Aldor. And even more important you should neither write 
.spad nor .as but .as.pamphlet. Anything else is considered "old" style 
and will not survive the 30 year horizon.

Yes, yes. Everyone knows that this is a pain at first. But after your 
code starts getting more complex this pain becomes a relief. So how 
complex must your code be till you think to switch to a literate 
programming style? If you already start in a non-literate fashion, you 
will find it harder to switch later. So why not start with LP from the 
beginning???

> ---------------------- This is the function
>
> R     ==> INT
> F     ==> Fraction R
> P     ==> Polynomial F

Use SparseUnivariatePolynomial.

> M     ==> Matrix(F)

Why don't you use SquareMatrix? Does a minimal polynomial make sense for 
something non-square?

> 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

I am surprised that the compiler will know to which symbol that ++ 
documentation is to be associated.

>  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

Anyway, you have not told that the function is in an .input file and the 
domain in a .spad file. Now it should become quite clear to you that if 
two programs (the interpreter for .input and the SPAD compiler for 
.spad) handle your more or less identical code pieces, the semantics is 
not the same. The reason is that the interpreter tries hard to figure 
out what you could have meant. It is dangerous to rely on that guessing. 
In fact, that is more or less quite similar to programming without 
types. The spad compiler does not guess. So you have to be a bit more 
careful to get the types right yourself, ie, you have to be very 
explicit with coercing your objects to the right type.

Please start to use Aldor, after a while you'll recognise that the error 
messages help you to figure out which coerce you have to add. There is 
no way around. If you want to use Axiom for programming, then you MUST 
thing about the types of your objects. If you don't want to do this then 
AXIOM is not for you. Or I better should say, ...then you are bound to 
use .input files and never be able to write complex programs. It's 
simply up to you whether you want to take the pain of thinking also 
about the types. Believe me, after you are through the initial hard 
part, you'll love types.

Sorry, that I did not help you by pointing you exactly to the error. I 
would have to browse through hyperdoc myself, since I am not overly 
familiar with the Axiom library. But switching to Aldor would certainly 
help you.

Ralf