Karim Belabas on Tue, 29 Jul 2008 01:52:42 +0200

[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

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

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

```