Re: Fast test for integer-rooted polynomial

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

Hello,

  here's my take on this:

1) about the characteristic polynomial: I don't see how to improve that step.
For such small matrices as the one presented, it doesn't help to compute it
modulo a prime (we use a multimodular algorithm which in practice stops after
a single prime...)

2) instead of factor(p) + couting irreducible factors, you could use
nfroots(,p) [ notice the comma "," ! ] then sum valuations. But since
your roots are heavily constrained, you can just directly sum valuations :

test2(a) = my(p = charpoly(a)); sum(i = 0, #a, valuation(p, x-i)) == #a;

On your given input I get

(11:52) gp > for(i=1,10^5,test(a))
time = 4680 ms.

(11:52) gp > for(i=1,10^5,test2(a))
time = 1604 ms.

So you can expect to test ~ 62k matrices per second. A direct C
implementation would probably optimize the valuations test to essentially 0
(i.e. a global speedup of about a factor 2, see below).

Another possibility would be to replace p = charpoly(a) by p /= gcd(p,p')
and test that if you expect many multiple roots. But it doesn't gain much
in your example, since timings decompose as follows:

(11:52) gp > for(i=1,10^5, ) \\ loop overhead
time = 4 ms.
(11:52) gp > for(i=1,10^5, p = charpoly(a))
time = 656 ms.
(12:53) gp > for(i=1,10^5, sum(i = 0, 9, valuation(p, x-i)))
time = 877 ms.

(12:53) gp > for(i=1,10^5, q = gcd(p,p'))
time = 373 ms.
(12:53) gp > for(i=1,10^5, sum(i = 0, 9, valuation(q, x-i)))
time = 464 ms.

Cheers,

    K.B.


* Denis Simon [2019-04-10 11:30]:
> Hi Gordon, 
> 
> writing your test as a function, this looks like 
> 
> test(a) = my(p,f); p=charpoly(a);f=factor(p);vecmax(apply(poldegree,f[,1]))==1 
> 
> If your matrices have a high probability to fail the test, then you could 
> do a first check modulo a prime : 
> 
> testp(a) = my(p,f); p=charpoly(a);f1=factormod(p,31,1);if(vecsum(f1[,2])<#a,0,f=factor(p);vecmax(apply(poldegree,f[,1]))==1) 
> 
> Denis SIMON. 
> 
> > De: "Watson Ladd" <watsonbladd@gmail.com>
> > À: "Gordon Royle" <gordon.royle@uwa.edu.au>
> > Cc: "pari-users" <pari-users@pari.math.u-bordeaux.fr>
> > Envoyé: Mercredi 10 Avril 2019 03:15:48
> > Objet: Re: Fast test for integer-rooted polynomial
> 
> > On Tue, Apr 9, 2019, 6:10 PM Gordon Royle < [ mailto:gordon.royle@uwa.edu.au |
> > gordon.royle@uwa.edu.au ] > wrote:
> 
> >> Dear Pari/GP users
> 
> >> This is a very simple task, but I want to do it billions of times, and so I want
> >> it to be as fast as possible.
> 
> >> I have written a suitable Pari/GP routine to do it, but just naively, and I
> >> wonder if there is a clever way to do it.
> 
> >> INPUT: Billions of symmetric matrices of smallish order (say at most 16 x 16)
> >> each of which has small integer entries.
> 
> >> OUTPUT: The input matrices that have only integer eigenvalues.
> 
> >> From the nature of the problem, I know that any integer eigenvalues lie in
> >> {0,1,2, …, n} where n is the order of the matrix.
> 
> >> What I do currently:
> 
> >> - get Pari to compute the characteristic polynomial (charpoly())
> >> - factor the characteristic polynomial
> >> - check that the factors are all linear
> 
> >> Here is an example:
> 
> >> a=
> >> [7,0,-1,-1,-1,-1,-1,-1,-1; 0,7,-1,-1,-1,-1,-1,-1,-1; -1,-1,7,0,-1,-1,-1,-1,-1;
> >> -1,-1,0,7,-1,-1,-1,-1,-1; -1,-1,-1,-1,7,0,-1,-1,-1; -1,-1,-1,-1,0,4,0,0,0;
> >> -1,-1,-1,-1,-1,0,7,-1,-1; -1,-1,-1,-1,-1,0,-1,7,-1; -1,-1,-1,-1,-1,0,-1,-1,7]
> >> p=charpoly(a)
> >> f=factor(p)
> >> print(f)
> 
> >> produces
> 
> >> [x - 9, 2; x - 8, 3; x - 7, 2; x - 4, 1; x, 1]
> 
> >> showing that the eigenvalues are just 9, 9, 8, 8, 8, 7, 7, 4, 0 which are
> >> integers as required.
> 
> >> This works fine, but of course, faster is better, which leads me to my
> >> questions:
> 
> >> Q1: Is there perhaps a better way to compute the char poly, such as the other
> >> methods supplied by Pari, perhaps using the flag=4 option to take advantage of
> >> the integer entries?
> 
> >> Q2: Is there perhaps a better way to check that the polynomial is a product of
> >> linear factors without actually factorising it? For example, for each i in
> >> {0,1,2, …, n}, check if p(i) = 0 and then divide by (x-i) if so. Repeat.
> 
> >> Q3: Or is there a way to check the eigenvalues without even computing the
> >> characteristic polynomial? (I looked at qfjacobi but it returns floating point
> >> numbers for the eigenvalues)
> 
> >> Thanks in advance for any suggestions
> 
> >> Gordon
> 
> > Why not compute the characteristic polynomial and then evaluate at 16 points and
> > compute multiplicaties via the derivative?

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