Finding Wilson Primes

Question

A Wilson prime is a prime number p such that p² divides (p−1)!+1. The first three Wilson Primes are 5, 13 and 563 and the fourth is larger than \$2×10^{13}\$. I was curious as to how much memory/processing power it would take to calculate the Wilson Primes using a brute force type method.

When I was coding this I kept efficiency in mind, and for this reason I used Stacks rather than Vectors.

#include <iostream>
#include <stack>

using namespace std;

stack<int> findPrimesUnder(int limit){ //stack is used for O(k) complexity where k is constant
    stack<int> primes;

    for(int i = 2; i <= limit; i++){
        int numFactors = 0;
        for(int j = i-1; j > 1; j--){//anything divided by itself is 1 so can be excluded for efficiency by subtracting 1 from i
                                     //also, anything divided by 1 is itself so is excluded for efficiency;
            if (float(i)/j == i/j){
                numFactors++;
            }
        }
        if(numFactors == 0)
            primes.push(i);

    }

    return primes;
}

int factorial(int n){
    if(n == 1){
        return 1;
    }
    else{
        return n*factorial(n-1);
    }
}

int main(){
    int limit;
    cout << "Enter limit: "; cin >> limit;

    stack<int> primes = findPrimesUnder(limit);
    stack<int> wilsonPrimes;

    bool descriptive = false;

    while(!primes.empty()){
        unsigned long long int firstWilsonCheck = (factorial(primes.top()-1)+1)/primes.top();//((p-1)! + 1)/p is always an int where p is prime
        double secondWilsonCheck = double(firstWilsonCheck)/primes.top();

        if(secondWilsonCheck == int(secondWilsonCheck))
            wilsonPrimes.push(primes.top());

        if(descriptive){
            cout << "Prime: " << primes.top() << endl;
            cout << "First Check: " << firstWilsonCheck << endl;
            cout << "Second Check: " << secondWilsonCheck << endl;

            cout << "------------------------" << endl;
        }


        primes.pop();
    }


    cout << "These are the Wilson Primes under " << limit << ":" << endl;
    while(!wilsonPrimes.empty()){
        cout << wilsonPrimes.top() << endl;
        wilsonPrimes.pop();

    }   

    return 0;
}

Output:

Enter limit: 100
These are the Wilson Primes under 100:
5
13
37
41
43
47
53
59
61
67
71
73
79
83
89
97

This code calculates the first two Wilson Primes fine (5 and 13), however it cannot calculate any after 13 because the numbers just get too big, I assume.

My questions are:

1. How could this code be made to be more efficient?
2. How can I increase precision on larger numbers?

Answer 1

When working with number theory problems like this, efficiency improvements are usually found more in math than in code.

Let us deal with the first problem: For large primes, the factorials overflow. We can deal with this by finding the remainder the factorial leaves on dividing by p^2 instead of the factorial itself.

int factorialremainder(int n, int mod) {
    if(n == 1) {
        return 1;
    } else {
        return (n*factorial(n-1))%mod;
    }
}

Now how do we check whether a number is Wilson prime? Easy. Just check whether (factorialremainder(p-1, p*p) + 1)%(p*p) == 0.

That brings us to the next point, you seem to not know how to use the remainder operation. You can find the remainder on dividing x by y using x%y. x is divisible by y iff x%y == 0. You don't need to do (float(i)/j == i/j) gymnastics anymore.

Finally, your primality testing can be improved a lot. Most numbers have small factors so it actually makes sense to run the loop over j in the reverse direction. Also, you don't need to go till i-1, square root of i is enough. So to find primes you can just do

bool isprime(int n) {
  // Assuming n > 1
  for(int i = 2; i*i <= n; ++i) {
    if(n%i == 0) return false;
  }
  return true;
}

Actually, since you are interested in finding all primes upto a certain range, you should try using some method like Sieve of Eratosthenes.

Also, stack is a very weird data structure to use for the problem. Why not just use a vector?