I’m releasing beta code for a numbers game written in Python. The code was written by the Chatbot ChatGPT. Update: I’ve added a C language program after the Python program. The program interactively solicits the user to enter non-negative integers. Numbers with a hundred or more digits are acceptable. I don’t know how to run Python programs on Windows machines directly. The Python code can be run from a Linux installation. If you don’t have Linux available but would be interested in trying out Linux on Windows, I would recommend trying the Windows Subsystem for Linux method. This can be found as official Microsoft-distributed software. There are WSL1 and WSL2 with WSL2 the more mature version, so I’d recommend Windows Subsystem for Linux version 2. Along with the Microsoft-distributed WSL software, there are various Linux distributions to choose from. I like the Ubuntu distribution, so I’d recommend that one if you’re new to Linux. Managing software packages in Ubuntu or other Linux distributions is covered in articles on web sites such as Geeks for Geeks. I’ll assume that you have a Python interpreter installed on your Linux system, and that you have basic familiarity with Linux and the shell (bash shell).
I’m copying below an interactive session of the numbers game in Python:
david@HPLaptop:~/mar31/apr26/goldcode$ python3 numbers_game05a.py
Enter k (non-negative integer): 55
q= 13
q= 3
0
Enter k (non-negative integer): 137
q= 11
1
Enter k (non-negative integer): 343
0
Enter k (non-negative integer): 1001
q= 13
0
Enter k (non-negative integer): 1999
q= 11
1
Enter k (non-negative integer): 1729
q= 101
q= 149
0
Enter k (non-negative integer): 1000000007
q= 13
q= 19
1
Enter k (non-negative integer): 2047
q= 3
0
Enter k (non-negative integer): 8191
q= 29
q= 53
q= 2
1
Enter k (non-negative integer): 6
q= 13
q= 19
0
Enter k (non-negative integer): 26599
q= 13
q= 53
q= 37
q= 107
0
Enter k (non-negative integer): 170141183460469231731687303715884105727
q= 29
q= 23
1
Enter k (non-negative integer): 1
0
Enter k (non-negative integer): 0
0
Enter k (non-negative integer): 2
1
Enter k (non-negative integer): 5
q= 5
1
Enter k (non-negative integer): 257
q= 5
1
Enter k (non-negative integer): 65537
q= 3
1
Enter k (non-negative integer): 667
q= 2
0
Enter k (non-negative integer): 485329
q= 5
0
Enter k (non-negative integer): 917701
q= 29
q= 5
q= 191
q= 101
q= 37
0
Enter k (non-negative integer): 119029
q= 29
q= 53
q= 13
0
Enter k (non-negative integer): 1746289
q= 13
q= 103
q= 37
q= 73
0
Enter k (non-negative integer): 1879861
q= 11
0
Enter k (non-negative integer): 5818561
q= 47
q= 191
q= 89
0
Enter k (non-negative integer):
The program reports whether a number k passes the test with a single 1 or 0 on a line by itself. For example, 667 doesn’t pass the test, and we see that also 667 is composite because 23×29 = 667. If you find a composite that passes the test, I’d be very interested to know. Also, if for some reason a prime fails the test, I’d be interested although this would probably be a bug in the code. Also, if you have features that you would like to see in the numbers game, please leave a comment. Below, I’m reproducing the Python code from the file numbers_game05a.py in /home/david/mar31/apr26/goldcode.
#!/usr/bin/env python3
import sys
from sympy import isprime, integer_nthroot
from math import gcd
# Multiply two polynomials a, b mod f(x)=x^3 - p*x - p, with coeffs mod k.
def poly_mul(a, b, p, 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])
# fast exponentiation of alpha (0,1,0) to exponent e in the above ring
def poly_pow(e, p, k):
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):
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_test(k):
if k < 2:
return 0
if k == 2:
return 1
if is_perfect_power(k):
return 0
# if isprime(k):
# return 1
# 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 ≡ p (mod k)
rem = k % p
t = 0
while True:
q = rem + t*p
if isprime(q):
break
t += 1
# print q
print("q=",q)
# 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 # <-- important: 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
s1 = tuple((alpha[i] + beta[i] + gamma[i]) % k for i in range(3))
if s1 == (0,0,0):
return 1
# s2
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 1
# s3
tmp = poly_mul(alpha, beta, p, k)
s3 = poly_mul(tmp, gamma, p, k)
if s3 == ((p % k), 0, 0):
return 1
# f was irreducible but none of the Vieta relations held
return 0
# Exhausted all m up to 1000 without finding an irreducible f
return 0
if __name__ == "__main__":
try:
while True:
line = input("Enter k (non-negative integer): ").strip()
if not line:
continue
try:
k = int(line)
if k < 0:
print("Please enter a non-negative integer.", file=sys.stderr)
continue
except ValueError:
print("Invalid input; enter an integer.", file=sys.stderr)
continue
print(vieta_test(k))
except (EOFError, KeyboardInterrupt):
print("\nGoodbye!")
sys.exit(0)
I now have a C language program that finds the numbers between 2 and 10,000 that satisfy the first of the three criteria, which is known as the trace or s_1 criterion.
#include <stdio.h>
#include <stdint.h>
#include <math.h>
// Computes the numbers that pass criterion 1, from START to STOP
#define START 2
#define STOP 10000
int premier(int64_t n)
{
int64_t d=2;
int pass=1;
if(n<2)
{
return 0;
}
if(n==2)
{
return 1;
}
while(d*d<=n)
{
if((n%d)==0)
{
pass = 0;
break;
}
d++;
}
return pass;
}
int64_t gcd(int64_t m, int64_t n)
{
int64_t r,s,nexts;
if(m<n)
{
r = n;
s = m;
}
else
{
r = m;
s = n;
}
while(s>0)
{
nexts = r%s;
r = s;
s = nexts;
}
return r;
}
// Multiply two 3×3 matrices X and Y modulo n, storing the result in P.
// - P, X, Y: int64_t[3][3] arrays
// - n: the modulus (may be up to 2^63–1)
void mul_mat_mod(int64_t P[3][3], int64_t X[3][3], int64_t Y[3][3], int64_t n)
{
// For each row i of X
for (int i = 0; i < 3; i++) {
// For each column j of Y
for (int j = 0; j < 3; j++) {
// Use a 128‐bit accumulator to avoid overflow
__int128 sum = 0;
// Compute the dot‐product of row i of X with column j of Y
for (int k = 0; k < 3; k++) {
sum += ((__int128)X[i][k]) * ((__int128)Y[k][j]);
}
// Reduce modulo n
int64_t v = (int64_t) (sum%n);
// If negative, bring back into [0, n–1]
if (v < 0) v = v + n;
// Store the low 64 bits of the (now reduced) sum
P[i][j] = v;
}
}
}
// compute R = A^e mod n by repeated squaring
void pow_mat_mod(int64_t R[3][3], int64_t A[3][3], int64_t e, int64_t n)
{
int64_t base[3][3];
int64_t result[3][3];
int64_t newresult[3][3];
int64_t newbase[3][3];
// 1) init result = identity matrix I₃
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
result[i][j] = (i == j) ? 1 : 0;
// 2) copy A -> base
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
base[i][j] = A[i][j] % n; // reduce if you like
// 3) exponentiation by squaring
while (e > 0) {
if (e & 1) {
// result = result * base mod n
mul_mat_mod(newresult, result, base, n);
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
result[i][j] = newresult[i][j];
}
}
}
// base = base * base mod n
mul_mat_mod(newbase, base, base, n);
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
base[i][j] = newbase[i][j];
}
}
e >>= 1;
}
// 4) copy result -> R
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
R[i][j] = result[i][j];
}
int64_t integerroot(int64_t n,int exp)
{
int64_t start=(int64_t) roundl(expl(logl((long double)n)/((long double)exp)));
return start;
}
int perfect_power(int64_t n)
{
int64_t nroot;
int exp;
int val;
int i;
int pow;
val = 0;
for(exp=2;exp<60;exp++)
{
nroot = integerroot(n,exp);
pow = 1;
for(i=0;i<exp;i++)
{
pow=pow*nroot;
}
if(pow==n)
{
val = 1;
break;
}
}
return val;
}
int test_Vieta(int64_t k)
{
int64_t m;
int64_t p;
int64_t q;
int64_t a;
int64_t b;
int irr;
int done = 0;
if(perfect_power(k))
{
return 0;
}
for(m=0;m<1001;m++)
{
p = m*m + m + 7;
if(premier(p))
{
if(gcd(p,k)==1)
{
q = k%p;
while(!premier(q))
{
q = q + p;
}
irr = 1;
for(a = 0; a < q; a++)
{
__int128 B = (__int128)a*a*a - (__int128)p*a - (__int128)p;
b = (int64_t)(B%q);
while(b<0)
{
b = b + q;
}
b = b%q;
if(b==0)
{
irr = 0;
break;
}
}
if(irr==1)
{
done = 1;
}
}
}
if(done==1)
{
break;
}
}
if(done == 1)
{
int64_t A[3][3];
int64_t Ak[3][3];
int64_t Ak2[3][3];
int64_t trace[3][3];
int pass = 1;
int i,j;
for(i=0;i<3;i++)
{
for(j=0;j<3;j++)
{
A[i][j] = 0;
}
}
A[1][0] = 1;
A[2][1] = 1;
A[0][2] = p;
A[1][2] = p;
pow_mat_mod(Ak,A,k,k);
pow_mat_mod(Ak2,Ak,k,k);
for(i=0;i<3;i++)
{
for(j=0;j<3;j++)
{
trace[i][j] = (A[i][j]+Ak[i][j]+Ak2[i][j])%k;
}
}
for(i=0;i<3;i++)
{
for(j=0;j<3;j++)
{
pass=pass&&(trace[i][j] == 0);
}
}
return pass;
}
else
{
return 0;
}
}
int main(void)
{
int64_t k;
int j;
for(k=START;k<STOP;k++)
{
if(test_Vieta(k))
{
printf("%ld\n", k);
}
}
return 0;
}
