# 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 <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?
• FWIW, as a general rule, the cache-friendly std::vector will be faster than the cache-busting std::list (which is what std::stack uses by default). In this case it doesn't matter, because you're not doing anything significant with the data structure; you're just iterating over it once and then you're done. In cases where it doesn't matter, I strongly recommend std::vector, so that you get into the habit of using vector in general. Jan 14 '17 at 20:07
• Okay! Point taken! I used a stack because I learnt in computing class that stacks are more efficient than vectors. However, like many things I learn in computing class, this isn't useful, practical, or even true in reality. Jan 14 '17 at 23:10
• Lists might be theoretically more efficient, but in practice, on modern architectures, they are less efficient than arrays (like vector) because of "locality of reference" (Googleable). These are the kinds of things that comp sci classes don't teach you, but matter a lot in the real world. Accessing memory is expensive, so you need to use the cache to your advantage. A contiguous structure like an array does that, whereas a non-contiguous list does precisely the opposite, setting yourself up to fight against the cache. Jan 15 '17 at 18:07

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?

• Thanks for the great feedback! I'll be implementing this, however I still feel like I'm gonna run into the same problem of numbers getting too big. How can I accurately deal with BIG numbers in C++? Jan 14 '17 at 23:18
• How would you deal with BIG numbers on paper? Just take that approach and translate it into C++. Raziman has already shown you one very good idea, which is to keep the numbers from getting super big in the first place. But since you know you're going to have to deal with numbers on the order of 10^14, which is 2^42, you might start thinking about ways to keep even those small numbers smaller. Maybe the Chinese Remainder Theorem could help; I haven't thought about it much. Jan 14 '17 at 23:25