Finding Wilson Primes 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 &lt;iostream&gt;
#include &lt;stack&gt;

using namespace std;

stack&lt;int&gt; findPrimesUnder(int limit){ //stack is used for O(k) complexity where k is constant
    stack&lt;int&gt; primes;

    for(int i = 2; i &lt;= limit; i++){
        int numFactors = 0;
        for(int j = i-1; j &gt; 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 &lt;&lt; "Enter limit: "; cin &gt;&gt; limit;

    stack&lt;int&gt; primes = findPrimesUnder(limit);
    stack&lt;int&gt; 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 &lt;&lt; "Prime: " &lt;&lt; primes.top() &lt;&lt; endl;
            cout &lt;&lt; "First Check: " &lt;&lt; firstWilsonCheck &lt;&lt; endl;
            cout &lt;&lt; "Second Check: " &lt;&lt; secondWilsonCheck &lt;&lt; endl;

            cout &lt;&lt; "------------------------" &lt;&lt; endl;
        }


        primes.pop();
    }


    cout &lt;&lt; "These are the Wilson Primes under " &lt;&lt; limit &lt;&lt; ":" &lt;&lt; endl;
    while(!wilsonPrimes.empty()){
        cout &lt;&lt; wilsonPrimes.top() &lt;&lt; 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:


How could this code be made to be more efficient?
How can I increase precision on larger numbers?
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 &gt; 1
  for(int i = 2; i*i &lt;= 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?