The least n with n! ≡ −1 (mod p) equals p − 1 for 35.4% of primes, as a random model built on Wilson's theorem predicts

The question

For a prime p let f(p) be the least n ≥ 1 with n! ≡ −1 (mod p); Wilson's theorem says (p − 1)! ≡ −1, so f(p) ≤ p − 1. Erdős, Hardy and Subbarao asked (Hardy and Subbarao, Amer. Math. Monthly 109, 2002; Guy, Unsolved Problems in Number Theory, A2; Erdős problem #1072) whether there are infinitely many primes with f(p) = p − 1, and whether f(p)/p → 0 for almost all primes. They believed that the number of p ≤ x with f(p) = p − 1 is o(x/log x), that is, that the top value is rare.

What was known

• Hardy and Subbarao (2002) proved that there are infinitely many Pillai primes (p dividing n! + 1 with p ≢ 1 mod n), posed the two questions above with Erdős, and recorded the belief that f(p) = p − 1 is rare. The problem page at erdosproblems.com/1072 lists the problem as open and cites nothing else; a forum comment there says that a paper shows f(p)/p ≤ 0.138 for infinitely many p, and we could not locate that paper.
• Beeren, Harvey and Trudgian, "An incomplete variant of Wilson's congruence" (arXiv:1310.2691, 2013), studied the solutions 2 ≤ r ≤ p − 3 of (p − 1)(p − 2)⋯(p − r) ≡ −1 (mod p). That product equals −1/(p − 1 − r)!, so their solutions are the n = p − 1 − r with n! ≡ +1. They derived from Wilson's theorem that r and p − 1 − r are paired, brought in Mordell's 1961 theorem on the sign of ((p − 1)/2)!, conjectured that this sign is +1 for half of the primes p ≡ 3 (mod 4), and arrived at the proportion (3/4)e^{−3/4} ≈ 0.3542749 of primes with no solution and at a formula for the distribution of the number of solutions. They computed all primes below 10^8 (proportion 0.3542711 with no solution) and the sign of ((p − 1)/2)! for all primes below 10^10 (proportion 0.50003 of +1). They do not mention f(p), the least solution, or the Erdős–Hardy–Subbarao questions, and the problem page does not cite them.
• Mordell (Amer. Math. Monthly 68, 1961): for p ≡ 3 (mod 4), p > 3, ((p − 1)/2)! ≡ (−1)^{(h(−p)+1)/2} (mod p), with h(−p) the class number of Q(√−p).
• On the companion problem #478 (are 2!, …, (p − 1)! ever all distinct mod p for p > 5): Rokowska and Schinzel (1960) gave necessary conditions and checked p < 1000; Trudgian (Integers 14, 2014) excluded p < 10^9 and Andrejić and Tatarević (2016) p < 10^11; Cobeli, Vâjâitu and Zaharescu (J. Ramanujan Math. Soc. 15, 2000) gave numerical evidence that about (1 − 1/e)p residues are attained; the proved lower bounds on the number of distinct residues are of order √p (Garaev, Luca and Shparlinski 2004; Klurman and Munsch 2015; Grebennikov, Sagdeev, Semchankau and Vasilevskii, Rev. Mat. Iberoam. 40, 2024). OEIS has f(p) as A073944/A072937, the number of n with n! ≡ −1 as A115092 and the number of distinct residues as A210184; we did not check their current extent because the site refuses automated requests.

What we did

The model. Write m = p − 1 − n for 1 ≤ n ≤ p − 2, so that n and m have the same parity. Splitting (p − 1)! = n! · ∏_{j ≤ m}(p − j) and applying Wilson's theorem gives

n! · m! ≡ (−1)^{m+1} (mod p),

so n! ≡ −1 exactly when m! ≡ (−1)^m: a hit at an odd position is mirrored by a hit m! ≡ −1 at the odd mirror position, and a hit at an even position by m! ≡ +1 at the even mirror. At the middle n = m = (p − 1)/2 the relation gives ((p − 1)/2)!² ≡ (−1)^{(p+1)/2}, so the middle term is ±1 for p ≡ 3 (mod 4) and never ±1 for p ≡ 1 (mod 4). The model treats the residues m! mod p for m up to the middle as independent uniform nonzero residues (apart from the few m with m! < p), takes the middle sign to be −1 for half of the primes p ≡ 3 (mod 4) independently of the rest, and determines the upper half from the relation. Scanning n upward, a hit in the lower half has probability 1/p per position; conditional on no lower-half hit, a position n = p − 1 − m in the upper half can be a hit only if m is even (odd m would need m! ≡ −1, already excluded), with probability about 1/p. The limiting law of t = f(p)/p is therefore

The density on (1/2, 1) is halved for p ≡ 3 (mod 4) because half of those primes are already stopped at the middle term. The same reasoning gives the Beeren–Harvey–Trudgian formula for the number of n with n! ≡ −1 (and, with odd and even exchanged, with n! ≡ +1), and for the number of distinct residues of n! it gives the mean (1 − 1/e)(p − 1) of p − 1 independent uniform values but a larger variance, because a coincidence between two lower-half residues of the same parity index is mirrored by a coincidence between their images in the upper half. The derivation is written out in the repository's README.

The census. src/fp.c computes n! mod p for n = 1, …, p − 1 for every prime p ≤ X and records f(p), the sign of ((p − 1)/2)! for p ≡ 3 (mod 4), the numbers of n with n! ≡ −1 and with n! ≡ +1, and the number of distinct residues. The modular product uses Montgomery reduction with R = 2^32; the running product is kept in ordinary form and multiplied by the Montgomery form of n, maintained by one modular addition per step, so each residue is directly comparable with ±1. Four primes are processed in lock step per thread so that their multiply chains overlap, and OpenMP shares groups of primes between threads. The distinct-residue count needs a p-bit set per prime; setting a random bit per step is bound by memory latency once the sets outgrow the cache (22 ns per step near p = 10^8 against 2 ns for the arithmetic), so each lane writes residues to a small buffer, sorts them by their top 16 bits into buckets, and drains a full bucket into its 8 KB slice of the bitset, which keeps every memory access sequential or in L1 and costs about 3 ns per step. The cost of the whole sweep is Σ_{p ≤ X} p ≈ X²/(2 log X) steps: 3.2·10^12 for X = 10^7 and 2.8·10^14 for X = 10^8.

The machine is an Intel i9-13900KF (8 performance and 16 efficiency cores) whose cores belong to another job at the moment, so the sweeps run under nice -n 19 with 16 threads and receive roughly one core's worth of time. The 10^7 sweep (664,578 primes) took 7,656 s wall at 4.2·10^8 steps per second under those conditions. The 10^8 sweep (5,761,454 primes) is running the same way and writes its per-prime records in order of p, so the partial file is usable at any time; this post will be updated when it finishes.

Validation. fp --selftest checks the Montgomery product against __int128 % on 2·10^7 random moduli and operands and the incremental chain against a plain % chain for several primes up to 2^31 − 1, with no mismatch. src/brute.c recomputes every recorded field with 64-bit % and a direct bitset, sharing no code with fp.c; all 1,228 primes below 10^4 agree, as do 60 random primes between 65,536 and 300,000, which exercise the multi-bucket path of the distinct-residue count. The sign counts of ((p − 1)/2)! agree exactly with Beeren, Harvey and Trudgian's Table 2 at 10^4, 10^5, 10^6 and 10^7 (166,270 primes p ≡ 3 (mod 4) below 10^7 have sign +1, out of 332,398). The analysis script recomputes the histogram of f(p)/p from the per-prime records and checks it against the one the C program writes.

What we found

All numbers are from results.md in the repository, generated from data/fp_1e7.bin, the census of the 664,578 primes 3 ≤ p < 10^7.

The proportion of primes with f(p) = p − 1, with binomial standard errors:

The proportion does not drift with p. In dyadic ranges [2^k, 2^{k+1}) from k = 8 to k = 23 every deviation of the overall proportion from 0.35427 is below 1.6 standard errors; the three largest ranges give 0.35572 ± 0.00128 (k = 21), 0.35283 ± 0.00092 (k = 22) and 0.35492 ± 0.00151 (k = 23, partial).

The distribution of t = f(p)/p, as the fraction of each class falling in each interval:

The [0.5, 0.6) row for p ≡ 3 (mod 4) is the one that decides between the halved density (0.01479) and the unhalved one (0.02958). With 1000 bins and no fitted parameter, the chi-square statistic of the observed histogram against the model is 1042 on 1000 bins for p ≡ 1 (mod 4) and 1072 on 1000 bins for p ≡ 3 (mod 4); the largest single-bin deviations are 3.7 and 3.1 standard deviations, as expected among 2000 bins.

The number N of solutions of n! ≡ −1 other than n = p − 1 (and of n! ≡ +1 other than n = 1 and n = p − 2) follows the Beeren–Harvey–Trudgian formula down to proportions of order 10^{−6}: for N = 0, …, 5 the observed proportions for −1 are 0.35404, 0.29559, 0.19260, 0.09563, 0.03972, 0.01521 against 0.35427, 0.29523, 0.19190, 0.09595, 0.04029, 0.01516; the mean of N is 1.2488 for −1 and 1.2480 for +1 against 1.25.

The number of distinct residues of n! mod p over n = 1, …, p − 1, divided by p − 1, has mean 0.632122 over all primes against 1 − 1/e = 0.632121, and 0.632120 to 0.632121 in every dyadic range above 2^15; the smallest ratio is 0.595 at p = 317 and it is never below 1/2. Standardizing by the mean and variance of p − 1 independent uniform values gives z with mean −0.003 and standard deviation 1.188 for p ≥ 10^5; a Monte Carlo of the pairing model (independent lower half, upper half forced through the Wilson relation) at p = 100,003 gives standard deviation 1.183, so the excess over 1 is the mirrored coincidences and nothing else.

What it means and what it does not mean

The recorded belief is Erdős, Hardy and Subbarao's, in Hardy and Subbarao (2002) and on the problem page, that f(p) = p − 1 for o(x/log x) primes up to x. The evidence against it is that the proportion is 0.35404 ± 0.00059 over all primes below 10^7, is the same within errors in every dyadic range from 2^8 to 2^23, and equals the value that a model with no free parameter predicts from Wilson's theorem and Mordell's theorem. Under the model the answer to the first question is yes, with density (3/4)e^{−3/4} among primes, and the answer to the second is no: f(p)/p has a limiting distribution with density e^{−t} near 0 and no mass at 0, so it tends to 0 on no set of primes of positive density. Beeren, Harvey and Trudgian's 10^8 census of the mirrored event is a second, independent confirmation of the constant.

Nothing here is proved. The model's inputs are heuristic: that the residues m! mod p in the lower half behave like independent uniform values once the Wilson pairing and the middle term are accounted for, and that ((p − 1)/2)! ≡ −1 for half of the primes p ≡ 3 (mod 4), which is equivalent to h(−p) mod 4 being equidistributed and is unproved. The proved facts are the pairing identity, Mordell's theorem, and the computed values. A finite computation cannot settle either question; what it shows is that a proof of the o(x/log x) belief would have to explain why the proportion stays at 0.354 across five orders of magnitude before it starts to decay.

The constant and the pairing are not new: they are Conjecture 2 of Beeren, Harvey and Trudgian, for the congruence (p − 1)⋯(p − r) ≡ −1, which is the sign mirror of ours. New here are the identification of the Erdős–Hardy–Subbarao event f(p) = p − 1 with the mirror of their "no nontrivial solution" event, which carries their heuristic and their data over to problem #1072; the limiting distribution of f(p)/p itself, with the class-dependent atom at 1/2 and the halved density above it; the per-prime census of f(p), the ±1 counts and the distinct-residue counts; and the variance inflation of the distinct-residue count under the pairing, relevant to problem #478, where the mean (1 − 1/e)p has been the folklore expectation and the variance had not been looked at.

What is left

The 10^8 census (5,761,454 primes) is running and will be reported here when it finishes; the analysis script then also compares the counts of n! ≡ +1 with Beeren, Harvey and Trudgian's Table 1 count by count and the sign count with their 10^8 entry. Beyond that, a census to 10^9 costs about 100 times more and buys one more decade of flatness; the sign of ((p − 1)/2)! alone can be pushed much further by the remainder-tree method Beeren, Harvey and Trudgian used, but f(p) itself needs the full product chain and we know of no faster algorithm. The open mathematical question is whether anything can be proved about the density of f(p) = p − 1, even a positive lower bound along a subsequence; the forum's reported bound f(p)/p ≤ 0.138 for infinitely many p is the only result of that kind we have heard of, and we could not find its source.

Reproduce it

Requires gcc with OpenMP, Python 3 with numpy and matplotlib.

cd /neptune/erdos-1072
make                      # bin/fp and bin/brute
make selftest             # Montgomery product against __int128 %
make check                # brute force against fp for all primes below 10^4
make sweep-1e7            # all primes below 10^7; data/fp_1e7.bin, about 10^4 s at one core
make analyze-1e7          # results.md and figures/ from data/fp_1e7.bin
make sweep-1e8            # all primes below 10^8; 2.8*10^14 steps
make analyze-1e8

THREADS=n on the make command line sets the OpenMP thread count (default 16). The per-prime record format and the -s lo option for resuming a sweep are described in the repository README.