[arndt@jjj.de: bugfix for charpoly with finite fields]

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

----- Forwarded message from Joerg Arndt <arndt@jjj.de> -----

From: Joerg Arndt <arndt@jjj.de>
To: Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr>
Subject: bugfix for charpoly with finite fields

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).

The routines use the Hessenberg-based method, straight from Cohen's
book.  Might be a good idea to forward to the list(s).

Note the timings!

cheers,  jj


P.S.: URL is http://www.jjj.de/pari/charpoly2.inc.gp




\\% Compute characteristic polynomial for matrices over finite fields
\\ Author: Joerg Arndt
\\ online at http://www.jjj.de/pari/
\\ version: 2008-July-27 (13:11)

/*

This is a bugfix for pari/gp:
charpoly(matid(3)*Mod(1,2))
  *** charpoly: impossible inverse modulo: Mod(0, 2).

Also quite fast for small characteristic and dense matrices:
with a dense 256 x 256 matrix and p==2 we have

  charpoly(M);  \\ workaround 1: computation over integers
    ***   last result computed in 3mn, 24,765 ms.

  matdet(Mod(1,2)*('x*matid(n)-M));  \\ workaround 2: use determinant
    ***   last result computed in 11mn, 21,971 ms.

  charpoly_ff(M,2);
    ***   last result computed in 3,236 ms.

  charpoly_2(M);
    ***   last result computed in 1,937 ms.

*/


charpoly_ff(H, p=2)= \\ p must be prime, H a square matrix
{
    local(n, P, t);
    H = mathess( Mod(1,p)*H );
    n = matsize(H)[1];
    P = vector(n+1);  P[1+0] = 1;
    for (m=1, n,
        P[1+m] = ('x-H[m,m])*P[1+m-1];
        t = 1;
        for (i=1, m-1,
            t *= H[m-i+1, m-i];
            if ( 0==t, break(); );  \\ good with small characteristic
            P[1+m] -= ( t * H[m-i,m] * P[1+m-i-1] );
        );
    );
    return( lift( P[n+1] ) );
} /* ----------- */


charpoly_2(H)= \\ H must be a square matrix
{
    local(n, P);
    H = lift( mathess( Mod(1,2)*H ) );
    n = matsize(H)[1];
    P = vector(n+1);  P[1+0] = 1;
    for (m=1, n,
\\        P[1+m] = ('x-H[m,m])*P[1+m-1];
        P[1+m] = P[1+m-1] << 1;
        if ( H[m,m],  P[1+m] = bitxor(P[1+m], P[1+m-1]) );
\\        t = 1;
        for (i=1, m-1,
\\            t = H[m-i+1, m-i];
            if ( 0==H[m-i+1, m-i], break(); );  \\ good with small characteristic
\\            if ( H[m-i,m],   P[1+m] -= P[1+m-i-1] );
            if ( H[m-i,m],   P[1+m] = bitxor(P[1+m], P[1+m-i-1]) );
        );
    );
\\    return( lift( P[n+1] ) );
    return( Pol( binary( P[n+1] ) ) );
} /* ----------- */



----- End forwarded message -----

    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]
`