I had an AI write a Python program to test random numbers in the numbers game from my last post. On the command line, one gives two numerical arguments: the number of iterations and the upper limit on the numbers. So python3 vietatest.py 20 1000000000000 will run 20 tests on random numbers in the range 1 to 1,000,000,000,000. For each number tested, the program tells us whether the number is prime or composite and whether it passes at least one of the three criteria in the Vieta test or numbers game from my last post. Here is some sample output:
$ python3 vietatest.py 20 1000000000000
Running 20 tests with random numbers up to 1000000000000
found composite 466924670922 test failed
found composite 633750414149 test failed
found composite 577008152114 test failed
found composite 724421856951 test failed
found composite 784494692235 test failed
found composite 330084995397 test failed
found composite 372438535220 test failed
found composite 999664171310 test failed
found prime 255867532901 test passed
found composite 304354078213 test failed
found composite 330127649409 test failed
found composite 947476385992 test failed
found prime 699576688403 test passed
found composite 241267665178 test failed
found prime 156245948081 test passed
found composite 86559569357 test failed
found composite 144520444347 test failed
found composite 506757159939 test failed
found composite 82031664421 test failed
found composite 789164509604 test failed
It’s possible to pipe the ouput through the “grep” command, for example to show only passed tests:
$ python3 vietatest.py 1000 1000000000000 | grep pass
found prime 592829144477 test passed
found prime 779670098717 test passed
found prime 800337025871 test passed
found prime 161685138929 test passed
found prime 905429021597 test passed
found prime 626344428433 test passed
found prime 598064193049 test passed
found prime 914714148151 test passed
found prime 45784555741 test passed
found prime 399256172219 test passed
found prime 17909282701 test passed
found prime 10401940531 test passed
found prime 775921544899 test passed
found prime 936841424767 test passed
found prime 891480479789 test passed
found prime 201978541213 test passed
found prime 794911375771 test passed
found prime 198629484367 test passed
found prime 861181539391 test passed
found prime 725218381319 test passed
found prime 475368217741 test passed
found prime 210428758213 test passed
found prime 118762312003 test passed
found prime 207032828089 test passed
found prime 424524614243 test passed
found prime 352771073807 test passed
found prime 479098207057 test passed
found prime 297058960439 test passed
found prime 963102649409 test passed
found prime 919950731213 test passed
found prime 430325011513 test passed
found prime 605072504321 test passed
found prime 51898246489 test passed
As of this writing, the range of integers up to 1000 billion or 1,000,000,000,000 has not been exhaustively tested. Below I’ve copied the program listing. If you should find a composite that passes or a prime that fails the test, I’d be very interested. Should you have ideas for improvements, please don’t hesitate to leave a comment. Update I’ve added a Pari/gp script that also tests random numbers, in the range 1 to 1,000,000,000,000. It is 8 times faster than the Python script. Numbers can be tested individually with the test function, and as an example, randomtries(40,1) will show just the numbers where there is agreement between Pari/gp’s isprime function and the experimental test function (at present, it only checks one of the three criteria: the trace criterion). In the same vein, randomtries(40,0) will show just the numbers where isprime and the experimental test function disagree. The discussion of the Pari/gp script continues after the Python program listing.
#!/usr/bin/env python3
import sys
import random
from sympy import isprime, integer_nthroot
from math import gcd
def poly_mul(a, b, p, k):
"""
Multiply two polynomials a, b mod f(x)=x^3 - p*x - p, with coeffs mod k.
a, b are tuples (c0, c1, c2) representing c0 + c1 x + c2 x^2
"""
r = [0, 0, 0, 0, 0] # up to x^4
for i, ai in enumerate(a):
for j, bj in enumerate(b):
r[i+j] = (r[i+j] + ai * bj) % k
# reduce x^4 -> p*x + p and x^3 -> p*x + p
# x^4 term: r[4] x^4 => add r[4]*p to r[2], add r[4]*p to r[1]
r[2] = (r[2] + r[4]*p) % k
r[1] = (r[1] + r[4]*p) % k
# x^3 term: r[3] x^3 => add r[3]*p to r[1], add r[3]*p to r[0]
r[1] = (r[1] + r[3]*p) % k
r[0] = (r[0] + r[3]*p) % k
return (r[0], r[1], r[2])
def poly_pow(e, p, k):
"""
Fast exponentiation of alpha (0,1,0) to exponent e in the polynomial ring
"""
result = (1 % k, 0, 0)
base = (0, 1 % k, 0)
while e > 0:
if e & 1:
result = poly_mul(result, base, p, k)
base = poly_mul(base, base, p, k)
e >>= 1
return result
def is_perfect_power(n):
"""
Check if n is a perfect power
"""
if n < 4:
return False
# try exponents 2..log2(n)
max_e = n.bit_length()
for e in range(2, max_e+1):
root, exact = integer_nthroot(n, e)
if exact:
return True
return False
def vieta_check(k):
"""
Perform the Vieta criterion check on number k
Returns True if k passes the test (likely prime), False otherwise
"""
if k < 2:
return False
if k == 2:
return True
if is_perfect_power(k):
return False
# Try m=0..1000
for m in range(0, 1001):
p = m*m + m + 7
if not isprime(p):
continue
if gcd(k, p) != 1:
continue
# find first prime q ≡ k (mod p)
rem = k % p
t = 0
q_search_limit = 20000 # Limit search for q based on benchmarking
while True:
q = rem + t*p
if isprime(q):
break
t += 1
if t > q_search_limit:
# If we can't find a suitable q within the limit, try next p
break
if t > q_search_limit:
continue
# If f is reducible mod q, skip to next p
reducible = False
for a in range(q):
if (a*a*a - p*a - p) % q == 0:
reducible = True
break
if reducible:
continue # try the next p, do not fail immediately
# f is irreducible mod q → now do the Vieta checks
alpha = (0, 1 % k, 0)
beta = poly_pow(k, p, k)
gamma = poly_pow(k*k, p, k)
# s1 criterion: alpha + beta + gamma = 0
s1 = tuple((alpha[i] + beta[i] + gamma[i]) % k for i in range(3))
if s1 == (0, 0, 0):
return True
# s2 criterion: alpha*beta + alpha*gamma + beta*gamma = -p
ab = poly_mul(alpha, beta, p, k)
ag = poly_mul(alpha, gamma, p, k)
bg = poly_mul(beta, gamma, p, k)
s2 = tuple((ab[i] + ag[i] + bg[i]) % k for i in range(3))
if s2 == ((-p) % k, 0, 0):
return True
# s3 criterion: alpha*beta*gamma = p
tmp = poly_mul(alpha, beta, p, k)
s3 = poly_mul(tmp, gamma, p, k)
if s3 == ((p % k), 0, 0):
return True
# f was irreducible but none of the Vieta relations held
return False
# Exhausted all m up to 1000 without finding an irreducible f
return False
def main():
"""
Main function to handle command line arguments and run tests
"""
# Check command line arguments
if len(sys.argv) != 3:
print("Usage: python vietatest.py <num_tests> <max_value>")
print("Example: python vietatest.py 1000 1000000000000000")
sys.exit(1)
try:
num_tests = int(sys.argv[1])
max_value = int(sys.argv[2])
except ValueError:
print("Error: Both arguments must be integers")
sys.exit(1)
if num_tests <= 0 or max_value <= 0:
print("Error: Both arguments must be positive integers")
sys.exit(1)
print(f"Running {num_tests} tests with random numbers up to {max_value}")
# Perform the tests
for i in range(num_tests):
# Generate a random number between 1 and max_value
rnd_num = random.randint(1, max_value)
# Check if it's prime using sympy's isprime
is_prime = isprime(rnd_num)
# Run the Vieta check for all numbers
test_result = vieta_check(rnd_num)
# Report the result
if is_prime:
if test_result:
print(f"found prime {rnd_num} test passed")
else:
print(f"found prime {rnd_num} test failed")
else:
if test_result:
print(f"found composite {rnd_num} test passed")
else:
print(f"found composite {rnd_num} test failed")
if __name__ == "__main__":
main()
Here is some sample output with the Pari/gp randomtries function with the flag set to 1:
? randomtries(40,1)
200290928302 0 0
456881108018 0 0
95169935946 0 0
651528487876 0 0
414058790821 0 0
105133864048 0 0
432715469186 0 0
510704827866 0 0
683170736239 0 0
780679789321 0 0
746367834623 0 0
809667053694 0 0
909471654786 0 0
761866423075 0 0
604459132405 0 0
830282576796 0 0
717620605875 0 0
268902180467 0 0
820713653027 0 0
9200143880 0 0
31048419449 0 0
590216686594 0 0
971007338045 0 0
353777357846 0 0
112332712073 1 1
128073138713 0 0
633452463673 0 0
131760037038 0 0
611093610179 0 0
641610787791 0 0
464248737010 0 0
860345709215 0 0
668014508677 0 0
291447470440 0 0
82876941693 0 0
301963294162 0 0
456609754989 0 0
192469659642 0 0
534234173995 0 0
221310402134 0 0
And here are the definitions of the 5 functions in the script. test is the main experimental testing function. It calls findp , findq and score. findp finds the p-parameter for a given test number, while findq finds the q-parameter for a given number. These last two functions are usually very fast, although sometimes the search will take hundreds of tries. The other function called by test is the scoring function “score”. This is the workhorse of the entire setup, and does computations with 3×3 companion matrices of cubic polynomials. My plan is to give a quite detailed conceptual description of the experimental test in the coming weeks. I probably will not give proofs, but instead state the ideas needed to run the experimental test. Please feel free to leave a comment.
test(n)={if(n<2,return(0));if(ispower(n),return(0));p1=findp(n);q1=findq(n);ans=score(n,p1,q1);return(ans)}
findp(k)={for(m=0,1000,p=m^2+m+7;if(isprime(p)&&((k%p)>0),r=k%p;while(!isprime(r),r=r+p);if(polisirreducible(Mod(x^3-p*x-p,r)),q=r;break)));return(p)}
findq(k)={for(m=0,1000,p=m^2+m+7;if(isprime(p)&&((k%p)>0),r=k%p;while(!isprime(r),r=r+p);if(polisirreducible(Mod(x^3-p*x-p,r)),q=r;break)));return(q)}
score(n,p,q)={A=matrix(3);A[2,1]=Mod(1,n);A[3,2]=Mod(1,n);A[1,3]=Mod(p,n);A[2,3]=Mod(p,n);An=A^n;An2=An^n;tr=A+An+An2;s=0;pass=1;for(i=1,3,for(j=1,3,pass=pass&&tr[i,j]==Mod(0,n)));if(pass,s=1);return(s)}
randomtries(k,fl)={for(iter=1,k,num=random(10^12);a=isprime(num);b=test(num);if(fl==(a==b),print(num," ",a," ",b)))}
I’m adding a PDF file that provides some mathematical details relevant to the Pari/gp script functions test, findp, findq and score.
