What's the most efficient way to find Wilson Prime number using python?

Question

I have written a code of Wilson prime numbers and my code is working for most of the numbers but it's giving OverflowError: int too large to convert to float for very large numbers. Is there any way to write Wilson prime number code for very large numbers. The main problem is for checking Wilson prime Wilson primes it should satisfy the following condition. Where P represents a prime number.

Then ((P-1)! + 1) / (P * P) should give a whole number. And as you can see factorials are involved in this procedure, so for very large numbers it's pretty difficult.

My Code :

def am_i_wilson(n):
    import math
    n1 = math.sqrt(n)
    n1 = math.ceil(n1)
    c = 0

    def fact(n):
        num = 1
        for i in range(2,n+1):
            num = num*i
        return num

    if n <= 1:
        return False

    for i in range(2, n1 + 1):
        if n%i == 0:
            c+ = 1

    if c != 0:
        return False

    x = (fact(n-1)+1)/((n**2)*1.0)

    return x.is_integer()

In my code, I am returning True if the number is Wilson Prime else False. Here n is the number to check if it's Wilson prime or not.

Answer 1

I think this is the most efficient way

import math
def am_i_wilson(num):
    if num < 2 or not all(n % i for i in range(2, num)):
        return False
    return (math.factorial(num - 1) + 1) % (num ** 2) == 0

or you can try this too

import math
def am_i_wilson(n):
    if n <= 2:
        return False
    fact=math.factorial(n-1)
    #this conditional checks that the number is prime or not
    #this condition is called wilson theorem in number theory
    if (fact+1)%n==0:
        x = (fact+1)%(n**2)
        if x==0:
            return True
        else:
            return False
    else:
        return False

Answer 2

if anyone has any better method then please answer it.

Your program primarily relies on testing for primes and computing factorials. You separate out the factorial logic but embed an inefficient prime test -- it keeps testing remainders after it knows the number isn't prime! I would separate both out so that they can be tested and optimized independently of the Wilson prime test itself:

def factorial(n):
    number = 1

    for i in range(2, n + 1):
        number *= i

    return number

def am_i_prime(n):
    if n < 2:
        return False

    if n % 2 == 0:
        return n == 2

    for divisor in range(3, int(n ** 0.5) + 1, 2):
        if n % divisor == 0:
            return False

    return True

def am_i_wilson(n):
    return am_i_prime(n) and (factorial(n - 1) + 1) % n ** 2 == 0

A more efficient approach, given a fixed target to test up to, would be to implement a prime sieve and for each prime encountered, while you're computing the sieve, test if it's a Wilson prime.

Answer 3

I've been experimenting with prime sieves recently. I did a quick modification (i.e. hack) to one of them written by Robert William Hanks and came up with this. Output first:

$ ./wilson_primes.py 10000 [5, 13, 563]

...so I suspect the Wikipedia article about Wilson primes is correct ;-)

import sys
def fact(n):
    num = 1
    for i in range(2, n+1):
        num *= i
    return num

def is_wilson(n):
    return (fact(n-1)+1) % n**2 == 0

def rwh_primes1(n):
    """ Returns  a list of primes < n """
    sieve = [True] * (n/2)
    for i in range(3,int(n**0.5)+1,2):
        if sieve[i/2]:
            sieve[i*i/2::i] = [False] * ((n-i*i-1)/(2*i)+1)
    # return [2] + [2*i+1 for i in xrange(1,n/2) if sieve[i]]
    for i in range(1,n/2):
        if sieve[i]:
            p = 2*i + 1         # convert index to normal prime
            if is_wilson(p):    #
                yield p         #

if len(sys.argv) > 1:
    N = int(float(sys.argv[1]))
else:
    N = 10000           # default: 1e4 10,000

print [p for p in rwh_primes1(N)]

First I tried just the fact() function and was pleasantly surprised to see it can produce huge results. But it is very slow compared to the original prime sieve. Perhaps it could be made to run faster by remembering the last factorial computed and re-use that to skip part of next factorial computation.

EDIT

I changed fact() to remember its last result, as follows:

last_fact = 1
last_n = 1

def fact2(n):
    global last_fact, last_n
    num = last_fact
    for i in range(last_n+1, n+1):
        num *= i
    last_n = n
    last_fact = num
    return num

def is_wilson(n):
    return (fact2(n-1)+1) % n**2 == 0

That did speed it up quite a bit. cProfile shows that is_wilson() is now the bottleneck. I can't think of an easy way to make it faster.