charpoly using too much stack space !

[Gonzalo Tornaria <tornaria@mail.ma.utexas.edu>]

When computing the characteristic polynomial of big matrices (not so
much, in the example 50x50), PARI uses A LOT of stack space, which I
believe is not really necesary. I've written a gp function for
computing the characteristic polynomial (matrix_charpoly), using what
I believe is the same algorithm of charpoly (computing traces, etc.)

My GP function takes almost the same time as the PARI function (which
is sensible, since most part of the time is spent multiplying
matrices), but uses way less memory.

Here is a transcript for a 50x50 matrix.

---------
GP/PARI CALCULATOR Version 2.2.5 (development CHANGES-1.578)
[...]
parisize = 500000, primelimit = 10000000
? \r charpoly.gp
? #
   timer = 1 (on)
? cp50=matrix_charpoly(m50);
time = 2,910 ms.
? allocatemem(7000000)
? cp50==charpoly(m50)
  ***   the PARI stack overflows !
  current stack size: 7000000 (6.676 Mbytes)
  [hint] you can increase GP stack with allocatemem()

? allocatemem(7500000)
? cp50==charpoly(m50)
time = 2,690 ms.
%3 = 1
---------

As you can see, for the GP function it suffices to have parisize=500.000,
whereas for the PARI function parisize=7:000.000 it's not enough. It
works with 7:500.000 and shows that the result is the same (better
be). The matrix m50 is a random matrix with small coefficients.
The difference in time is less than 10%.

The difference in stack space needed seems to be more impressive for
bigger matrices (100x100, 150x150, etc.), and the difference in time
seems less noticeable.

Previously, I was thinking that computing characteristic polynomial of
250x250 matrices was not feasible (I don't know how much memory it's
needed for that using charpoly!), but I'm doing just that now... :-)

What's surprising is that, if I compile the GP function (with GP2C),
it doesn't get any faster, but it needs more stack space (even with -g)!!!

Here is the file charpoly.gp:

/*******************************************************************************/

matrix_charpoly(m) =
{ local(n, p, X, F, pk);
  n=matsize(m)[1];
  X=Pol([1,0]);
  p=1; F=m;
  for(k=1,n,
    pk=trace(F)/k;
    p=X*p-pk; F=m*(F-pk);
  );
  if(F!=0, error("matrix_charpoly"));
  p
}

m50=matrix(50,50,i,j,random(10000));

/*******************************************************************************/

Any comments?

Thanks,
Gonzalo
-- 
GM/CS/S d? a-- C++(+++)  UL+++(++++)  P++>+++  L+++>++++ E---
W-(+) N+(++) w--- O M-(--) V-(--) PGP(--) b++(+++) G++ e>++++

... A mathematician is a machine for converting coffee into theorems.