Re: Fast test for integer-rooted polynomial

[Gordon Royle <gordon.royle@uwa.edu.au>]

Dear All

Thank you all for your valuable suggestions.

I have done some more experimenting and timing based on your ideas and have come up with some that, experimentally, seems to be the best-so-far.

First: my matrices will almost never have integral eigenvalues and so the number of hits is very small and the number of misses very large; the largest case for which I have complete information is 11 x 11 matrices with 8400 hits out of a total of just over 1 billion matrices, which is just less than 0.001 %.

So I created a test set of 10^6 matrices of order 12 x 12, seeded it with a couple of winners (to check that the functions recognise hits).

Then I tried the following four tests 

test0(a) = my(p,f); p=charpoly(a);f=factor(p);vecmax(apply(poldegree,f[,1]))==1;

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

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);

test2(a) = my(p=charpoly(a)); for (i=0,#a,while(subst(p,x,i)==0,p\=(x-i))); p==1;

Test0 is Denis Simon’s version of my original program

Test1 is Karim Belabas’s version using valuation 

Testp is Denis Simon’s “check first modulo prime”

Test2 is my own version that just repeatedly factors out linear factors of the right type until there is nothing left

Outcome:

- the naive method using charpoly and factor took about 3m 30s

- the valuation method took MUCH LESS time, only 1m 20s

- the “check first modulo prime” took about 1m 30s 

- my “repeatedly divide” method also took 1m 20s 

The time taken in computing JUST the characteristic polynomials is 1m 12, and so the overhead after the characteristic polynomials is pretty 

So somehow we need to avoid computing the char poly if possible.

I have one idea on this… suppose A is a winning matrix. Then

- it has integer roots in [0..#A]

- it has 0 as a simple eigenvalue (all of them do), and so

- A+I has integer roots in [1..#A+1]

Therefore  the DETERMINANT of A+I must factorise into integers in the range [1..#A+1], so I can replace a characteristic polynomial calculation with a determinant calculation.

heuristictest3(a) = my(f=factor(matdet(a+matid(#a))));vecmax(f[,1])<=#a+1

This now takes only 55seconds to process my test set, but of course it produces some “false positives” which are matrices that pass the test based on the determinant, but actually don’t have integer eigenvalues at all. For my sample of 1000000 matrices, it found 21 false positives, which take negligible time to remove at the end.

So good progress so far… 

Thanks to all

Gordon

On 10 Apr 2019, at 9:09 am, Gordon Royle <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