Re: charpoly via Newton

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

* Pascal Molin [2012-10-24 13:01]:
> (sorry for writing in french)
> 
> En écrivant des sujets de tp, je découvre que calculer un polynôme
> caractéristique par les sommes de Newton n'est pas déraisonnable,
> en tous cas avec pari pour des matrices en caractéristique moyenne (0
> ou >dimension dans la fonction ci-dessous):
> 
> charpoly_newton(A) = {
>   my(d = matsize(A)[1]);
>   my(B = matid(d));
>   my(S = x*Ser(vector(d,i,B*=A;-trace(B)/i)));
>   polrecip(Pol(exp(S)));
> };
> 
> [gp-2.6 sur paridev] ============
> ? A = matrix(100,100,i,j,random(Mod(1,1009)));
> ? charpoly(A);
> time = 2,329 ms.
> ? charpoly_newton(A);
> time = 408 ms.

For some reason, charpoly() with default flag still doesn't inspect the input
object to choose the "best" available algorithm. It just defaults to
(failsafe, *slow*) Berkowitz algorithm. I'd consider this a bug.

On my laptop:

(13:37) gp > A = matrix(100,100,i,j,random(Mod(1,1009)));
(13:36) gp > charpoly_newton(A);  \\ assumes p > n
time = 364 ms.
(13:36) gp > charpoly(A,, 2); \\ Hessenberg form, valid for any p
time = 121 ms.

Can be further optimized:

(13:37) gp > install(ZM_to_Flm,GL)
(13:37) gp > install(Flm_charpoly,GL);
(13:37) gp > B = ZM_to_Flm(lift(A),1009);
(13:37) gp > Flm_charpoly(B,1009);
time = 20 ms.

N.B. Over Z, the situation is even less favourable

(13:37) gp > A = lift(A);
(13:37) gp > charpoly_newton(A);
time = 11,689 ms.
(13:37) gp > charpoly(A);   \\ Berkowitz
time = 1,932 ms.
(13:37) gp > charpoly(A,, 4); \\ modular + Hessenberg
time = 1,129 ms.

> L'algo peut être généralisé/amélioré (évaluation des traces avec BSGS,
> extension aux caractéristiques p<dim) : dites-moi
> si ça peut être intéressant.
> 
> Dans tous les cas, je souhaiterais proposer pour Pari un certain
> nombre d'opérations sur les polynômes (reconstruction à partir
> des sommes de Newton, produits composés, puissances symétriques --
> c'est initialement ça que je voulais --...), modalités et
> syntaxe à voir par exemple durant l'atelier de janvier.

Sounds interesting. Reconstruction from Newton sums are already
implemented (in fact, possibly more than once), but not exported to GP,
which only knows about polsym().

In any cas, all implantations (with tests) are welcome :-)

Cheers,

    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-bordeaux1.fr/~belabas/
F-33405 Talence (France)       http://pari.math.u-bordeaux1.fr/  [PARI/GP]
`