Depth-first search for multiplicative sequences

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