Re: bugfix for charpoly with finite fields

[Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr>]

Hello,

* Joerg Arndt [2008-07-27 11:39]:
> Attached a workaround/bugfix for pari/gp failing to compute the
> charpoly if the coefficients of the matrix are in a finite field.
> (Guess you are using the method involving division by i at step i,
> this is only good with characteristic zero, Cohen should mention
> this).

or dim(M) < characteristic, which covers some more cases.

> The routines use the Hessenberg-based method, straight from Cohen's
> book.  Might be a good idea to forward to the list(s).
[...]
> This is a bugfix for pari/gp:
> charpoly(matid(3)*Mod(1,2))
>   *** charpoly: impossible inverse modulo: Mod(0, 2).

  version 2.4.2 was in better shape already:

(01:37) gp > charpoly(matid(3)*Mod(1,2))
%1 = x^3 + Mod(1, 2)*x^2 + Mod(1, 2)*x + Mod(1, 2)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
(19:32) gp > ??charpoly
charpoly(A,{v = x},{flag = 3}):
[...]
   Let n be the dimension of A.

   If  flag  = 0,  the method used is the same as for computing the adjoint
matrix,   i.e. computing  the traces of the powers of A  (Le Verrier's method).
Assumes that n! is invertible; uses O(n^4) scalar operations.

   If  flag  =  1,   uses  Lagrange  interpolation which is usually the slowest
method. Assumes that n! is invertible; uses O(n^4) scalar operations.

   If  flag  =  2,  uses the Hessenberg form.   Assumes that the base ring is a
field. Uses O(n^3) scalar operations, but suffers from coefficient explosion.

   If flag = 3,  uses Berkowitz's division free algorithm,  valid over any ring
(commutative, with unit). Uses O(n^4) scalar operations.

   In practice one should use the default (Berkowitz) unless the base ring is a
field  where  coefficient explosion does not occur,  e.g. a finite field or the
reals. In which case the Hessenberg form should be preferred.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[ actually, this help text is from svn since I fixed a number of typos in the
doc: e.g. it stated that the Hessenberg algorithm used O(n^4) operations;
albeit true, this is a little off the mark ]

> Note the timings!

Hum.

1) while computing the characteristic polynomial of the Hessenberg matrix,
stopping the partial product as soon as we reach 0  is a clever idea 
in small characteristic, and costs nothing otherwise. I've incorported it.
Thanks !

2) charpoly(,,2) has been using way too much memory. I just fixed that.

(19:37) gp > M = matrix(256,256,i,j,random(2) * Mod(1,2));
(19:37) gp > charpoly(M,,2);
time = 4,296 ms.

3) just for fun, I've implemented that same Hessenberg algorithm over a small
prime field (the "Flm" pari subtype)

  install(Flm_charpoly,GL)
  install(ZM_to_Flm,GL)
  install(Flx_to_ZX,G)
  char(x, p) = Flx_to_ZX( Flm_charpoly(ZM_to_Flm(x,p), p) )

(19:37) gp > char(M, 2);
time = 388 ms. \\ ten times faster for a straightforward adaptation

4) since it had become easy, I added a chinese remaindering scheme for
matrices over Z, as charpoly(M,,4).

(19:37) gp > charpoly(lift(M),, 2);
time = 1mn, 1,315 ms.

[ going down to ~35s if I switched to a Monte Carlo algorithm: stop as
soon as result stabilizes; but I kept the rigorous algorithm ]

Thanks !

    K.B.
--
Karim Belabas, IMB (UMR 5251)  Tel: (+33) (0)5 40 00 26 17
Universite Bordeaux 1          Fax: (+33) (0)5 40 00 69 50
351, cours de la Liberation    http://www.math.u-bordeaux.fr/~belabas/
F-33405 Talence (France)       http://pari.math.u-bordeaux.fr/  [PARI/GP]
`