# Find primes using Wilson's Theorem

The following code is an instantiation of wilson's theorem. It uses a simple factorial and division theorem to check for primality. My issue is with the performance of the code, and there are two most prominent issues. Firstly, the cumtime for which the remainder is taken far exceeds the rest of the code; Secondly, there is a noticeable time crunch when adding 1 to the numerator due to what I believe is a problem when updating the reference to another object for very large integers.

Is there any possible way I could fix these issues, except, of course, making this program faster by coding in cpp?

n = 100000

def absolute(x: int):
if (x == 2):
a = Diction[x - 1] * x * (x + 1)
Diction.pop(x - 1)
Diction.update({(x + 1): a})
return 2
try:
a = Diction[x-1] * x
return a
finally:
a *= (x+1)
Diction.pop(x-1)
Diction.update({(x+1):a})

def main():
for i in range(3, n+1,2):
if (((absolute(i - 1) + 1 ) % i) == 0):
prime.append(i)

if __name__ == '__main__':
import time
import cProfile
start = time.perf_counter()

Diction = {1: 1}
prime = [2]
cProfile.run('main()')
end = time.perf_counter()
print(len(prime), end-start)


I expect that some other possible data holder, which uses a more optimal format, would work better in this instance - possibly numpy or an implementation of numba.

• See the addition at the bottom of my accepted answer below. Apparently by using gmpy2, performance increases by a huge factor (between 4 and 10). Nov 27, 2022 at 8:33

Apart from performance, there are some general review points:

• Don't put code after the if __name__ == '__main__': other than just calling main().
• Don't rely on globals to share data between functions (Diction, prime). Use function arguments and return values to exchange data.
• Use meaningful, descriptive variable names. prime contains primes (plural), so it's better to call it primes. Diction is referring to a datatype, which is not a descriptive variable name.
• Use capitals for module-wide constants (or, in this case: use function parameters instead of a constant n).
• cProfile already gives you the run time. No need to compute that again using time.

## Main topic: performance

If I understand your code correctly, you use Diction in combination with absolute to keep track of $$\(i - 1)!\$$. This seems overly complicated to me. You can keep track of the factorial inside the loop. This saves max_prime / 2 - 1 calls to dict operations pop and update.

Implementing all tips above, this gives the following (type hinted and docstring documented) code:

import cProfile

def generate_primes(max_prime: int) -> list:
"""Returns a list of primes until max_prime (included)

This function uses Wilkinson's theorem.

Args:
max_prime (int): upper limit for the primes to be calculated

Returns:
list: list of primes from 2 up until max_primes (included)
"""
primes = [2]
factorial = 2  # factorial of i-1

for i in range(3, max_prime + 1, 2):
if (factorial + 1) % i == 0:
primes.append(i)
factorial *= i * (i + 1)

return primes

def main():
max_prime = 100000

profiler = cProfile.Profile()
primes = profiler.runcall(generate_primes, max_prime)
profiler.print_stats()

print(len(primes))

if __name__ == '__main__':
main()


Running this gives roughly the same performance as your original code. :-(

As you already suspected: there are two calculations that probably account for most of the compute time: the modulo and the product of the factorial (which becomes rapidly very, very large). To measure this, I created two functions that perform these calculations, so cProfile can profile them. This gives the following output:

109591 function calls in 10.326 seconds

Ordered by: standard name

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
1    0.872    0.872   10.326   10.326 slow_prime_test.py:12(generate_primes)
49999    5.592    0.000    5.592    0.000 slow_prime_test.py:4(modulo_is_0)
49999    3.862    0.000    3.862    0.000 slow_prime_test.py:8(multiply)
9591    0.000    0.000    0.000    0.000 {method 'append' of 'list' objects}
1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}


So we're spending 54% of the time calculating the modulo and 37% of the time calculating (large) products.

Ideally, there is an algorithm that can calculate $$\((n - 1)! + 1) % n\$$ without actually calculating the factorial. The + 1 makes this very tricky and I don't know any algorithm that can do that. Without the + 1 we could use Legendre's formula to calculate $$\n! % p\$$.

You can probably speed your code up a bit using C++ or any other compiled high-performance language, but that fact remains that Wilson's theorem is not a very good performing theorem to generate primes. This is also stated on the Wikipedia page of the theorem:

Wilson's theorem has been used to construct formulas for primes, but they are too slow to have practical value.

## Use gmpy2 to speed up big integer calculations

I found the Python module gmpy2, which increases performance of big integer calculations. Just adding an from gmpy2 import mpz and replacing factorial = 2 with factorial = mpz('2') improves performance on my computer by a factor of 4.5.

Using gmpy2: True
9592
109591 function calls in 4.035 seconds

Ordered by: standard name

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
1    0.533    0.533    4.035    4.035 slow_prime_test.py:13(generate_primes)
49999    3.077    0.000    3.077    0.000 slow_prime_test.py:5(modulo_is_0)
49999    0.424    0.000    0.424    0.000 slow_prime_test.py:9(multiply)
9591    0.001    0.000    0.001    0.000 {method 'append' of 'list' objects}
1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

Using gmpy2: False
109591 function calls in 18.458 seconds

Ordered by: standard name

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
1    1.438    1.438   18.458   18.458 slow_prime_test.py:13(generate_primes)
49999   14.895    0.000   14.895    0.000 slow_prime_test.py:5(modulo_is_0)
49999    2.123    0.000    2.123    0.000 slow_prime_test.py:9(multiply)
9591    0.001    0.000    0.001    0.000 {method 'append' of 'list' objects}
1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}


On tio.run, I even get a performance improvement of around factor 10. I guess this depends on the computer platform. This TIO link includes the final code I used:

Try it online!

## Edit (2): using Legendre's formula

Since you're calculating all primes, we can use the fact that when testing if $$\n\$$ is prime, we already have all primes smaller than $$\n\$$. We can exploit this using Legendre's formula. Instead of calculating $$\((n - 1)! + 1) % n\$$, we can use the prime factorisation of $$\n - 1\$$ and the powmod function of gmpy2. This way, we avoid using very large integers, and (possibly) improve the execution speed. I implemented this, but sadly, it performs much slower than the gmpy2 optimised code mentioned above.

def wilsons_prime_test(n: Union[int, gmpy2.mpz],
primes: list) -> bool:
"""Checks if ((n - 1)! + 1 ) % n == 0 using the given set of primes < n

Args:
n (int or mpz): number of which to check if it is prime
primes (list): prime list which must give all primes < n

Returns:
bool: _description_
"""

result = gmpy2.mpz('1')

# Loop over all primes
for p in primes:
# Find largest power of p that divides n - 1
tmp_n = n - 1
largest_power = 0
while tmp_n:
tmp_n //= p
largest_power += tmp_n

# Multiply the result with prime ^ largest_power modulo n
result *= gmpy2.powmod(p, largest_power, n)

# Return if (result + 1) % n == 0
return (result + 1) % n == 0

• Why does the + 1 make it tricky? (a + 1) % b is just ((a % b) + 1) % b, or (a % b) + 1 unless that is equal to b, else zero. Nov 27, 2022 at 11:29
• You are right! That doesn't make it more difficult. I implemented Wilson's prime test using Legendre's formula, but it performs (much) worse than the gmpy2 optimised code. :-( Nov 27, 2022 at 14:16