Karim Belabas on Tue, 29 Jul 2008 01:43:42 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: bugfix for charpoly with finite fields |
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] `