/* -*- gp-script -*- */ \\% Compute characteristic polynomial for matrices over finite fields \\ Author: Joerg Arndt \\ License: GPL version 3 or later \\ online at http://www.jjj.de/pari/ \\ version: 2014-October-16 (18:32) /* 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 { my(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 { my(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 of file ====