A search for Wilson primes

A Wilson prime is a prime p such that (p-1)! = -1 mod p^2. We report on a search for Wilson primes up to 2 * 10^13, and describe several new algorithms that were used in the search. In particular we give the first known algorithm that computes (p-1)! mod p^2 in average polynomial time per prime.

💡 Research Summary

This paper presents a large-scale computational search for Wilson primes up to 2×10¹³, accompanied by significant algorithmic innovations. A Wilson prime is a prime number p satisfying the congruence (p-1)! ≡ -1 (mod p²), a stronger version of Wilson’s Theorem. Only three such primes (5, 13, and 563) were previously known.

The primary challenge in searching for Wilson primes is the computational cost of calculating the Wilson quotient w_p = ((p-1)! + 1)/p mod p for each prime p. Prior algorithms exhibited exponential time complexity in the input size log p, typically requiring time linear in p or proportional to √p.

The authors’ main theoretical contribution is the first algorithm that computes w_p in average polynomial time per prime. Specifically, Theorem 1 states that the Wilson quotients for all primes p ≤ N can be computed in total time N log³⁺ᵋ N. This implies an average time of essentially log⁴ p per prime. The algorithm cleverly exploits redundancies in computing (p-1)! for many different primes p. Instead of calculating (p-1)! modulo p² separately for each p, it computes factorials like N! modulo the product of p² for many primes simultaneously, then efficiently recovers the individual residues using a remainder tree and Fast Fourier Transform (FFT)-based integer arithmetic.

Since the algorithm of Theorem 1 has substantial memory requirements, the authors develop a more practical time-space trade-off described in Theorem 2. This algorithm computes w_p for primes in an interval (M, N] with tunable space complexity O(N-M). It consists of two stages: Stage 1 computes M! modulo the product of p² for primes in (M, N], and Stage 2 is a modified version of the Theorem 1 algorithm. This allows the search to be partitioned into manageable chunks based on available RAM.

Using the implementation of Theorem 2, the authors conducted a exhaustive search for all Wilson primes below 2×10¹³. The computation consumed over 1.1 million CPU core hours but found no new Wilson primes, confirming that 5, 13, and 563 remain the only known examples. To manage the enormous data output, they only saved the residues w_p for which |w_p| ≤ p/50000 (about 0.004% of primes), comprising approximately 27 million primes. These residues are made publicly available. The paper includes tables of primes with the smallest absolute |w_p| and the smallest relative |w_p/p|, showcasing “near misses.”

The paper also provides a thorough historical review of Wilson prime searches, from early manual computations in the 19th and early 20th centuries to computer-assisted searches in the late 20th century. It discusses heuristic expectations based on the assumed uniform distribution of w_p modulo p, which suggests Wilson primes are exceedingly rare but infinite in number. Furthermore, the authors describe and extend known number-theoretic identities that can reduce the computation of (p-1)! mod p² to computing a smaller factorial ((p-1)/e)! mod p² for certain primes p where e is a small divisor of p-1, building on earlier work by Crandall, Dilcher, and Pomerance.

Underpinning the algorithms are detailed complexity analyses for fundamental operations—integer multiplication/division and prime enumeration via a space-efficient sieve—in the Turing machine model.

In summary, this work represents a major advance in computational number theory by demonstrating how to structure large-scale calculations to achieve average-case polynomial time, enabling a search of unprecedented scale for a problem with inherent exponential difficulty when processing primes individually. It sets a new benchmark for the Wilson prime search and provides a framework for similar exhaustive computations in number theory.