Re: Counting real roots of integer polynomials

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

* Gordon Royle [2018-06-08 07:48]:
> However, I could not find a way to use polsturm effectively - the
> polsturm  in  Pari 2.9.5  needs a squarefree polynomial, and while the
> one in Pari  2.10.0 version will accept a non-square-free polynomial,
> it only returns distinct roots.
> 
> So to get the actual multiplicities, I could see no better way than
> using “factor” to first factorise the polynomial, and then deal with
> each factor separately.
> 
> This proved to be a bad idea, because for polynomials of the size that
> I am dealing with, “factor” takes about 5 times as long as just
> running “polrootsreal”. 

Yes, in pure GP you can't easily do better. Here's another way using
libpari's "squarefree factorization", which writes P in Z[X] as
\prod F[i]^E[i] with F[i] squarefree and pairwise coprime.

  install(ZX_squff,"G&");
  val(P) = 
  { my (x = variable(P), v = 0, v2);
    P /= x-1;
    v2 = valuation(P, x-2); if (v2, P /= (x-2)^v2);
    F = ZX_squff(P, E);
    for (i = 1, #F, v += E[i] * polsturm(F[i], [1,2]));
    [v, v2];
  }

  ? P = (x-1) * (2*x - 3)^3 * (x^2 - 2)^2 * (x-3)^3 * (x-2)^2;
  ? ZX_squff(P,E)
  %2 = [x - 1, x^3 - 2*x^2 - 2*x + 4, 2*x^2 - 9*x + 9]~
  ? E
  %3 = Vecsmall([1, 2, 3])
  ? val(P)
  %4 = [5, 2]

It should be faster than polrootsreal (and much faster than factor):

  ? for(i=1,10^4,polrootsreal(P))
  time = 2,903 ms.
  ? for(i=1,10^4,val(P))
  time = 1,071 ms.

Cheers,

    K.B.
--
Karim Belabas, IMB (UMR 5251)  Tel: (+33) (0)5 40 00 26 17
Universite de Bordeaux         Fax: (+33) (0)5 40 00 21 23
351, cours de la Liberation    http://www.math.u-bordeaux.fr/~kbelabas/
F-33405 Talence (France)       http://pari.math.u-bordeaux.fr/  [PARI/GP]
`