I had an AI optimize my C program that finds the numbers that pass the first criterion of the test mentioned in the numbers game blog post of April 26 (in a specified range). I tested the various versions repeatedly and compared the output of passers (numbers passing the 1st criterion) with files of bona fide prime numbers produced by Kim Walisch’s primesieve library function. All the primes in the range 2 to 40 billion pass the first criterion, and all passers in the range are prime numbers.
There are a couple of caveats, in that the optimized C code excludes from consideration those composites with small prime factors in the range 2 to 199. As a result, some pseudoprimes for the official test from the numbers game Python code might be missed, due to being divisible by a prime from 2 to 199. The second caveat is that in the optimized C code of the AI, there is a heuristic to limit the number of tries at finding the parameter q , given the number to be tested, k, and the parameter p:
// Add q search limit
int q_search_count = 0;
while (!premier(q)) {
q += p;
q_search_count++;
if (q_search_count > Q_SEARCH_LIMIT) {
q = 0; // Indicate failure to find q
break;
}
}
// Skip this p if q search failed
if (q == 0) continue;Nonetheless, I would expect only a small number of pseudoprimes for the first criterion up to 40 billion and for the official test of April 26 (and maybe none). Below, I’ve copied documentation of the computation, including the comparison with primesieve’s output:
david@HPLaptop:~/mar31/apr29/VietaGame/trace/optimized/testing$ ls -l
total 18816820
-rwxr-xr-x 1 david david 20656 Apr 29 11:06 a.out
-rw-r--r-- 1 david david 19268378624 Apr 29 19:18 my_primes_40G
-rw-r--r-- 1 david david 9813 Apr 29 11:06 vieta_trace_optimized_enhanced.c
david@HPLaptop:~/mar31/apr29/VietaGame/trace/optimized/testing$
real 512m59.071s
user 511m41.596s
sys 0m55.694s
[1]+ Done time ./a.out > my_primes_40G
david@HPLaptop:~/mar31/apr29/VietaGame/trace/optimized/testing$ primesieve 40000000000 -p > ref_primes_40G
david@HPLaptop:~/mar31/apr29/VietaGame/trace/optimized/testing$ sha256sum my_primes_40G
3f92725f59b61d4585d49b9ff9199b51ec1007dc38ff8c052ddb0c0b3ce06059 my_primes_40G
david@HPLaptop:~/mar31/apr29/VietaGame/trace/optimized/testing$ sha256sum ref_primes_40G
3f92725f59b61d4585d49b9ff9199b51ec1007dc38ff8c052ddb0c0b3ce06059 ref_primes_40G
-----------------------------------------------------------------------------------------------------------------------------
david@HPLaptop:~$ ls -l ~/mar31/apr29/VietaGame/trace/optimized/testing
total 39123192
-rwxr-xr-x 1 david david 20656 Apr 29 11:06 a.out
-rw-r--r-- 1 david david 20031049601 Apr 29 19:40 my_primes_40G
-rw-r--r-- 1 david david 20031049601 Apr 29 19:43 ref_primes_40G
-rw-r--r-- 1 david david 9813 Apr 29 11:06 vieta_trace_optimized_enhanced.c
-----------------------------------------------------------------------------------------------------------------------------
david@HPLaptop:~$ sha256sum ~/mar31/apr29/VietaGame/trace/optimized/testing/vieta_trace_optimized_enhanced.c
321816c43103d96963caea3e8051cbb4af7f813535a4c3494a52072a25a96278 /home/david/mar31/apr29/VietaGame/trace/optimized/testing/vieta_trace_optimized_enhanced.c
(The computation took about 8.5 hours.) Finally, I copy below the C source code file vieta_trace_optimized_enhanced.c :
#include <stdio.h>
#include <stdint.h>
#include <math.h>
#include <stdlib.h>
// Computes the numbers that pass criterion 1, from START to STOP
#define START 2
#define STOP 40000000000
// Optimization parameters
#define Q_SEARCH_LIMIT 20000 // Maximum iterations for q search
#define SMALL_PRIMES_COUNT 46 // Number of primes up to 200
// Pre-computed list of small primes up to 200
const int64_t small_primes[SMALL_PRIMES_COUNT] = {
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
157, 163, 167, 173, 179, 181, 191, 193, 197, 199
};
// Optimized primality test
int premier(int64_t n)
{
// Handle small cases efficiently
if(n < 2) return 0;
if(n == 2 || n == 3) return 1;
if(n % 2 == 0 || n % 3 == 0) return 0;
// Check divisibility only by numbers of form 6k±1
int64_t limit = sqrt(n);
for(int64_t d = 5; d <= limit; d += 6)
{
if(n % d == 0) return 0;
if(n % (d + 2) == 0) return 0;
}
return 1;
}
// Optimized GCD using binary GCD algorithm
int64_t gcd(int64_t m, int64_t n)
{
// Ensure positive values
m = llabs(m);
n = llabs(n);
// Base cases
if(m == 0) return n;
if(n == 0) return m;
// Find common factors of 2
int shift;
for(shift = 0; ((m | n) & 1) == 0; ++shift) {
m >>= 1;
n >>= 1;
}
// Remove remaining factors of 2 from m
while((m & 1) == 0)
m >>= 1;
// From here on, m is always odd
do {
// Remove factors of 2 from n
while((n & 1) == 0)
n >>= 1;
// Swap if needed so m <= n
if(m > n) {
int64_t t = m;
m = n;
n = t;
}
// Subtract and reduce
n -= m;
} while(n != 0);
// Restore common factors of 2
return m << shift;
}
// 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)
{
int64_t temp[3][3] = {0}; // Use temporary array to avoid overwriting input
// 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 in temporary array
temp[i][j] = v;
}
}
// Copy result to output
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
P[i][j] = temp[i][j];
}
}
}
// Optimized matrix exponentiation
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];
// 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 and reduce mod n
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
base[i][j] = A[i][j] % n;
// 3) exponentiation by squaring
while (e > 0) {
if (e & 1) {
// result = result * base mod n
mul_mat_mod(result, result, base, n);
}
// base = base * base mod n
if (e > 1) { // Only square if we'll use it again
mul_mat_mod(base, base, base, n);
}
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];
}
// Optimized integer root calculation
int64_t integerroot(int64_t n, int exponent)
{
// Use a more direct calculation for common exponents
if (exponent == 2) return (int64_t)sqrt(n);
if (exponent == 3) return (int64_t)cbrt(n);
// For other exponents, use logarithm method
return (int64_t)round(exp(log((double)n) / exponent));
}
// Optimized perfect power check
int perfect_power(int64_t n)
{
// Quick check for small numbers
if (n <= 3) return 0;
// Check if n is a perfect square first (most common case)
int64_t sqrt_n = integerroot(n, 2);
if (sqrt_n * sqrt_n == n) return 1;
// Check if n is a perfect cube
int64_t cbrt_n = integerroot(n, 3);
if (cbrt_n * cbrt_n * cbrt_n == n) return 1;
// Check higher powers only for larger numbers
if (n > 1000000) {
// Check 4th power
int64_t root4 = integerroot(n, 4);
int64_t pow4 = root4 * root4;
pow4 = pow4 * pow4;
if (pow4 == n) return 1;
// Check 5th power
int64_t root5 = integerroot(n, 5);
int64_t pow5 = root5 * root5 * root5 * root5 * root5;
if (pow5 == n) return 1;
// Check 6th power
int64_t root6 = integerroot(n, 6);
int64_t pow6 = root6 * root6 * root6;
pow6 = pow6 * pow6;
if (pow6 == n) return 1;
}
// For very large numbers, check higher powers
if (n > 1000000000) {
for (int exponent = 7; exponent <= 60; exponent++) {
int64_t nroot = integerroot(n, exponent);
// Only verify if the root is at least 2
if (nroot >= 2) {
// Manually compute power to avoid floating-point issues
__int128 result = 1;
for (int i = 0; i < exponent; i++) {
result *= nroot;
if (result > n) break; // Early exit if we exceed n
}
if (result == n) return 1;
} else {
// If root < 2, higher powers won't be perfect powers
break;
}
}
}
return 0;
}
int test_Vieta(int64_t k)
{
// Quick check: perfect powers cannot be valid
if (perfect_power(k)) {
return 0;
}
// Trial division by small primes
for (int i = 0; i < SMALL_PRIMES_COUNT; i++) {
int64_t p = small_primes[i];
if (k == p) {
// k is a small prime itself
return 1;
}
if (k > p && k % p == 0) {
// k is divisible by a small prime
return 0;
}
}
int64_t m;
int64_t p;
int64_t q;
int64_t a;
int64_t b;
int irr;
int done = 0;
// Search for suitable prime p
for (m = 0; m < 1001 && !done; m++) {
p = m*m + m + 7;
// Only check if p is prime
if (premier(p)) {
if (gcd(p, k) == 1) {
// Find prime q >= k mod p
q = k % p;
// Add q search limit
int q_search_count = 0;
while (!premier(q)) {
q += p;
q_search_count++;
if (q_search_count > Q_SEARCH_LIMIT) {
q = 0; // Indicate failure to find q
break;
}
}
// Skip this p if q search failed
if (q == 0) continue;
// Check irreducibility
irr = 1;
for (a = 0; a < q && irr; a++) {
__int128 B = (__int128)a*a*a - (__int128)p*a - (__int128)p;
b = (int64_t)(B % q);
// Ensure b is positive
if (b < 0) b += q;
if (b == 0) {
irr = 0;
break;
}
}
if (irr == 1) {
done = 1;
}
}
}
}
if (done == 1) {
int64_t A[3][3] = {0};
int64_t Ak[3][3];
int64_t Ak2[3][3];
int64_t trace[3][3];
// Initialize matrix A
A[1][0] = 1;
A[2][1] = 1;
A[0][2] = p;
A[1][2] = p;
// Calculate matrix powers
pow_mat_mod(Ak, A, k, k);
pow_mat_mod(Ak2, Ak, k, k);
// Check if trace is zero
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
trace[i][j] = (A[i][j] + Ak[i][j] + Ak2[i][j]) % k;
if (trace[i][j] != 0) return 0; // Early exit if any element is non-zero
}
}
return 1;
}
return 0;
}
int main(void)
{
int64_t k;
// Main loop to check numbers in the specified range
for (k = START; k < STOP; k++) {
// Skip small primes check for k < 200 to avoid redundant checks
if (k < 200) {
if (test_Vieta(k)) {
printf("%ld\n", k);
}
continue;
}
// Trial division by small primes before full test
int is_divisible = 0;
for (int i = 0; i < SMALL_PRIMES_COUNT; i++) {
int64_t p = small_primes[i];
if (k == p) {
// k is a small prime itself
printf("%ld\n", k);
is_divisible = 1;
break;
}
if (k % p == 0) {
// k is divisible by a small prime
is_divisible = 1;
break;
}
}
if (is_divisible) continue;
// Only run the full test if k passes the trial division
if (test_Vieta(k)) {
printf("%ld\n", k);
}
}
return 0;
}
