Depth-first search for multiplicative sequences
From Polymath Wiki
This C program performs depth-first search for multiplicative sequences.
In the source code, define N to be the limit of your search and C to be the limit on discrepancy. You also need to set MAX_N_PRIMES and MAX_N_DIVISORS (but you'll be prompted to do so if they aren't high enough).
The initial mapping of primes (f) is set within the main program: edit to taste! f[i] = true means that f maps the ith prime to -1, where the 0th prime is 2.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#ifndef _WIN32
# include <stdbool.h>
#else
# define bool char
# define true 1
# define false 0
#endif
#define N 3000 // limit of search
#define MAX_N_PRIMES 600
#define MAX_N_DIVISORS 600
#define C 3 // limit on sums
static unsigned int primes[MAX_N_PRIMES];
static unsigned int n_primes;
static bool is_prime[N+1];
static bool a[N+1][MAX_N_PRIMES]; // table of factorizations
static unsigned int s[N+1][MAX_N_DIVISORS]; // s[n] = list of divisors of n
static unsigned int l[N+1]; // l[n] = number of divisors of n
static bool f[MAX_N_PRIMES], g[MAX_N_PRIMES]; // f = initital values, g = values under test, true = sends to -1
static char x[N+1]; // sequence under test
static tot[N+1]; // totals table
static void seq(unsigned int M) { // calculate x corresponding to g, up to and including M
unsigned int n, i;
bool *fact;
char v = 1;
for (n = 1; n <= M; n++) {
fact = a[n];
for (i = 0; i < n_primes; i++) if (fact[i] && g[i]) v = -v;
x[n] = v;
}
}
static bool update_x(unsigned int i) { // update x from primes[i] to primes[i+1]-1, if discrep OK, return false if not
bool ok = true;
bool *fact;
char v;
unsigned int p1 = primes[i+1];
unsigned int n = primes[i];
unsigned int j;
char tt;
while ((n < p1) && ok) {
fact = a[n];
v = 1;
for (j = 0; j <= i; j++) if (fact[j] && g[j]) v = -v;
x[n] = v;
tt = tot[n-1] + x[n];
if (abs(tt) <= C) tot[n] = tt;
else ok = false;
n++;
}
return ok;
}
static void print_seq(FILE *fp, char *seq, unsigned int n) {
unsigned int j;
fprintf(fp, "\n ");
for (j = 1; j<=n; j++) {
fprintf(fp, "%s", (seq[j] == 1) ? "+" : "-");
if (j%81==0) fprintf(fp, "\n ");
}
fprintf(fp, "\n");
}
static bool check(char *seq, unsigned int n) {
bool ok = true;
unsigned int d = 1;
while ((d <= n) && ok) {
unsigned int i = d;
char tot = 0;
while ((i <= n) && ok) {
tot += seq[i];
if (abs(tot) > C) ok = false;
i += d;
}
d++;
}
return ok;
}
int main(int argc, char *argv[]) {
unsigned int i, max_i, min_back, max_fwd, next_g_i, n1, i0;
bool hope, ok, more;
{ // use Eratosthenes to find primes up to N
unsigned int p = 2;
bool ok = true;
unsigned int i, index;
is_prime[0] = false;
is_prime[1] = false;
for (i = 2; i<=N; i++) is_prime[i] = true;
while (ok) {
unsigned int i, q;
bool found;
for (i = 2*p; i <= N; i += p) is_prime[i] = false;
q = p+1;
found = false;
while ((q <= N) && !found) {
if (is_prime[q]) {
found = true;
p = q;
n_primes++;
}
q++;
}
if (!found) ok = false;
}
if (n_primes > MAX_N_PRIMES) {
fprintf(stderr, "Aborting. Increase MAX_N_PRIMES to at least %d.\n", n_primes);
return -1;
}
index = 0;
for (i = 2; i <= N; i++) {
if (is_prime[i]) {
primes[index] = i;
index++;
}
}
}
{ // initialize divisors table
unsigned int n, d;
for (n = 1; n <= N; n++) {
unsigned int d_count = 0;
for (d = 1; d <= N; d++) {
if ((n%d) == 0) {
if (d_count > MAX_N_DIVISORS) {
fprintf(stderr, "Aborting. Increase MAX_N_DIVISORS to at least %d (n = %d).\n", d_count, n);
return -1;
}
s[n][d_count] = d;
d_count++;
}
}
l[n] = d_count;
}
}
{ // make table of factorizations
unsigned int n;
for (n = 1; n <= N; n++) {
unsigned int i, k = n;
for (i = 0; i < n_primes; i++) {
unsigned int p = primes[i];
bool v = false;
while (k % p == 0) {
v = !v;
k /= p;
}
a[n][i] = v;
}
}
}
{ // set initial values of f, g, x
unsigned int i, p;
for (i = 0; i < n_primes; i++) {
p = primes[i];
if (p==3) f[i] = true;
else f[i] = (p%3==2);
if ((p==709) || (p==881) || (p==1201) || (p==1697) || (p==1999) || (p==2837) || (p==3323) || (p==3457)) f[i] = !f[i];
}
for (i = 0; i < n_primes; i++) g[i] = f[i];
seq(N);
}
{ // initialize totals table, t[n] = x[1]+x[2]+...+x[n]
unsigned int n;
for (n = 0; n <= N; n++) tot[n] = 0;
tot[1] = x[1];
}
i = 0;
max_i = 0;
min_back = 0;
max_fwd = 0;
hope = true;
while ((i < n_primes) && hope) {
if (i < min_back) {
min_back = i;
max_fwd = min_back;
printf("%d ", primes[min_back]);
fflush(stdout);
}
ok = update_x(i);
if (!ok) {
if (g[i] == f[i]) {
g[i] = !f[i];
ok = update_x(i);
}
}
if (ok) {
//if (i > max_fwd) {
// max_fwd = i;
// printf(" %d", i);
//}
if (i > max_i) {
max_i = i;
min_back = max_i;
n1 = primes[i+1]-1;
printf("\n\nmax_i = %d\n", max_i);
printf("max_n = %d\n", n1);
printf("verify: %s\n", check(x, n1) ? "pass" : "fail");
print_seq(stdout, x, n1);
printf("primes sent to -1:\n");
for (i0 = 0; i0 <= i; i0++) if (g[i0]) printf("%d, ", primes[i0]);
printf("\n\nbacktracking...\n");
}
i++;
}
else /* backtrack */{
while (g[i] == !f[i]) {
g[i] = f[i];
i--;
}
g[i] = !f[i];
if ((i == 0) && (g[i] == f[i])) hope = false;
}
}
printf("\n\nEND\n");
return 0;
}
