Solution to Project Euler Problem #37 in C++: Truncatable Primes

Question

I was recently looking at some of my old Project Euler solutions, and saw that my implementation for Problem 37 could be improved drastically.

The problem reads as follows:

The number 3797 has an interesting property. Being prime itself, it is possible to continuously remove digits from left to right, and remain prime at each stage: 3797, 797, 97, and 7. Similarly we can work from right to left: 3797, 379, 37, and 3.

Find the sum of the only eleven primes that are both truncatable from left to right and right to left.

NOTE: 2, 3, 5, and 7 are not considered to be truncatable primes.

Here is my C++ solution to the problem:

#include <vector>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <cstdlib>
#include <string>
std::vector<long long int> primesUpto(long long int limit) // Function that implements the Sieve of Eratosthenes
{
 std::vector<bool> primesBoolArray(limit, true);
 std::vector <long long int> result;
 
 primesBoolArray[0] = primesBoolArray[1] = false;
 long long int sqrtLimit = std::sqrt(limit) + 1;
 
 for (size_t i = 0; i < sqrtLimit; ++i)
 {
 if (primesBoolArray[i])
 {
 for (size_t j = (2 * i); j < limit; j += i)
 {
 primesBoolArray[j] = false;
 }
 }
 }
 for (size_t i = 0; i < primesBoolArray.size(); ++i)
 {
 if (primesBoolArray[i])
 {
 result.push_back(i);
 }
 }
 return result;
}
bool isTruncPrime(long long int number, const std::vector<long long int>& primeList)
{
 std::string numberString = std::to_string(number);
 for (int i = 1; i < numberString.size(); ++i)
 {
 std::string truncLeft = numberString.substr(0, i); // The truncated prime from the left
 std::string truncRight = numberString.substr(i, numberString.size() - 1); // The truncated prime from the right
 if (!
 (
 std::binary_search(primeList.begin(), primeList.end(), std::atol(truncLeft.c_str())) &&
 std::binary_search(primeList.begin(), primeList.end(), std::atol(truncRight.c_str()))
 ) // If either of the two truncated sides are not prime
 )
 {
 return false;
 }
 }
 return true; // All truncated parts are prime, so the number is a truncatable prime
}
int main()
{
 const std::vector<long long int> primesUptoMillion = primesUpto(1'000'000LL); // Represents all the primes up to 1 million
 
 int numberTruncatablePrimes = 0;
 long long int currentNumber = 9; // 2, 3, 5, and 7 are not included in the search for truncatable primes
 long long int truncatablePrimeSum = 0;
 
 while (numberTruncatablePrimes != 11) 
 {
 if (
 std::binary_search(primesUptoMillion.begin(), primesUptoMillion.end(), currentNumber) && // If the number itself is prime
 isTruncPrime(currentNumber, primesUptoMillion) // If the number is also a truncatable prime
 )
 {
 ++numberTruncatablePrimes; // Increase amount of truncatable primes
 truncatablePrimeSum += currentNumber; // Add the number's sum
 }
 currentNumber += 2; // Only odd numbers can be prime other than 2, so no need to look at every number
 }
 std::cout << truncatablePrimeSum << "\n";
}

Here is how I run the code:

g++ Problem037.cpp -std=c++14 -O2

Here is code for execution, timing, and output:

$ time ./a.out
748317
real 0m0.042s
user 0m0.040s
sys 0m0.004s

Here are the things I want from this Code Review:

• The function isTruncPrime() does a lot of string manipulation to verify its parameter. Is there a way to improve the algorithm?
• While the code does run pretty quickly, is there any code that slows down the rest significantly that I could improve on?
• Are there any better or more idiomatic ways to format/structure my program so it follows C++ style?

• polfosol: Maybe add the project-euler tag? I tried to do so, but my edit was too short and wasn't applied.
• Sigur: Do you have a ref for the only 11 part?
• Arnav Borborah: @polfosol, Apparently, whenever I add the project-euler, it gets replaced with programming-challenge. I suspect the tag must be a synonym.
• Arnav Borborah: @Sigur, If you are referring to the number of truncatable primes, the problem says to find the sum of the only 11 primes.
• Sigur: If I well understood, there are only 11 such numbers, right? So I asked for a ref to prove this fact. It is off topic question.

Answer 1

Simple things

You've misspelt std::size_t in a couple of places. I find I often do that, particularly when I've been writing C as well.

Conversion of primes table

I don't think there's a good reason to collapse primesBoolArray into a vector when it's perfectly usable as it is (and faster to use than a binary search). In fact, I'd go further, and make it live up to its name (since we always call the function with a constant):

// Get an array of bool - containing true at the prime indexes
template<std::size_t N>
std::array<bool,N> primesUpto()
{
 // Use the Sieve of Eratosthenes
 std::array<bool,N> primes;
 primes[0] = primes[1] = false;
 std::fill(primes.begin()+2, primes.end(), true);
 static const long long int sqrtLimit = std::sqrt(N) + 1;
 for (std::size_t i = 2; i < sqrtLimit; ++i)
 if (primes[i])
 for (std::size_t j = i+i; j < N; j += i)
 primes[j] = false;
 return primes;
}

String manipulation

We can use division to truncate numbers. Right to left is most obvious:

for (; n; n/=10)
 test(n);

Left to right can be done like this:

for (std::size_t x = 10; x < n; x*= 10)
 test(n%x);

This becomes

template<std::size_t N>
bool isTruncPrime(std::size_t n, const std::array<bool,N>& primes)
{
 for (std::size_t x = 10; x < n; x*= 10)
 if (!primes[n%x]) return false;
 for (; n; n/=10)
 if (!primes[n]) return false;
 return true;
}

These two changes reduce runtime on my system from 0.022 seconds to 0.007 seconds.

We can afford to search the whole problem space (assuming we already know there are no 7-digit truncatable primes):

int main()
{
 static const auto primes = primesUpto<1'000'000>();
 auto sum = 0ull;
 // Start at 11, because single-digit numbers aren't truncatable
 // Advance by 2, because even numbers are non-prime (except 2, but 2<11)
 for (std::size_t n = 11; n < primes.size(); n += 2) {
 if (isTruncPrime(n, primes))
 sum += n;
 }
 std::cout << sum << "\n";
}

This form allows us to parallelize the computation to (possibly) eke out a tiny bit more speed:

#pragma omp parallel for reduction(+:sum)

On my system, this makes it much slower - you'd need a larger problem space to benefit, I think.

Reduce the search space

Instead of testing all odd numbers, we can reduce the number of calls to isTruncPrime. We know that all such numbers end in a prime, but can't end in 2 or 5 (since they are composite), so the last digit can be only 3 or 7. Similarly, the first digit must be prime, and intermediate digits must be odd (but not 5):

static const auto b0 = { 0, 20000, 30000, 50000, 70000 };
static const auto b1 = { 10000, 30000, 70000, 90000 };
static const auto c0 = { 0, 2000, 3000, 5000, 7000 };
static const auto c1 = { 1000, 3000, 7000, 9000 };
static const auto d0 = { 0, 200, 300, 500, 700 };
static const auto d1 = { 100, 300, 700, 900 };
static const auto e0 = { 20, 30, 50, 70 };
static const auto e1 = { 10, 30, 70, 90 };
for (std::size_t a: { 0, 200000, 300000, 500000, 700000 }) {
 for (std::size_t b: a ? b1 : b0) {
 for (std::size_t c: a+b ? c1 : c0) {
 for (std::size_t d: a+b+c ? d1 : d0) {
 for (std::size_t e: a+b+c+d ? e1 : e0) {
 for (std::size_t f: { 3, 7 }) {
 const auto n = a + b + c + d + e + f;
 if (isTruncPrime(n, primes)) {
 sum += n;
 }
 }
 }
 }
 }
 }
}

This knocks another 70% off the runtime on my system.

---

Modified code

#include <array>
#include <algorithm>
#include <cmath>
#include <cstddef>
// Get an array of bool - containing true at the prime indexes
template<std::size_t N>
std::array<bool,N> primesUpto()
{
 // Use the Sieve of Eratosthenes
 std::array<bool,N> primes{};
 primes[0] = primes[1] = false;
 std::fill(primes.begin()+2, primes.end(), true);
 constexpr std::size_t sqrtLimit = std::sqrt(N) + 1;
 for (std::size_t i = 2; i < sqrtLimit; ++i)
 if (primes[i])
 for (std::size_t j = i+i; j < N; j += i)
 primes[j] = false;
 return primes;
}
template<std::size_t N>
constexpr bool isTruncPrime(std::size_t n, const std::array<bool,N>& primes)
{
 for (std::size_t x = 10; x < n; x*= 10)
 if (!primes[n%x]) return false;
 for (; n; n/=10)
 if (!primes[n]) return false;
 return true;
}

#include <iostream>
int main()
{
 static const auto primes = primesUpto<1'000'000>();
 auto sum = 0ull;
 static const auto b0 = { 0, 20000, 30000, 50000, 70000 };
 static const auto b1 = { 10000, 30000, 70000, 90000 };
 static const auto c0 = { 0, 2000, 3000, 5000, 7000 };
 static const auto c1 = { 1000, 3000, 7000, 9000 };
 static const auto d0 = { 0, 200, 300, 500, 700 };
 static const auto d1 = { 100, 300, 700, 900 };
 static const auto e0 = { 20, 30, 50, 70 };
 static const auto e1 = { 10, 30, 70, 90 };
 for (std::size_t a: { 0, 200000, 300000, 500000, 700000 }) {
 for (std::size_t b: a ? b1 : b0) {
 for (std::size_t c: a+b ? c1 : c0) {
 for (std::size_t d: a+b+c ? d1 : d0) {
 for (std::size_t e: a+b+c+d ? e1 : e0) {
 for (std::size_t f: { 3, 7 }) {
 const auto n = a + b + c + d + e + f;
 if (isTruncPrime(n, primes)) {
 sum += n;
 }
 }
 }
 }
 }
 }
 }
 std::cout << sum << "\n";
}

• Arnav Borborah: Wow, the idea that I don't have to separate the primes from non primes is actually brilliant! I am going to apply that in some of my other prime number solutions as well. What's really great about it is that it lets me get rid of all my find algorithms. I have a question, however. Why is the sqrt limit made static?
• Miles Budnek: std::fill isn't constexpr yet, so you shouldn't mark primesUpto as constexpr. Or you could replace the std::fill call with a simple loop to maintain constexpr-ness.
• Toby Speight: Thanks @Miles - I shouldn't assume that just because it compiles for me, it's right. I'll fix that now.

Answer 2

Warnings

It is usually a good idea to compile your code with all warnings enabled to see if anything is wrong. In your case, nothing really worrying but we may as well fix it.

Press ENTER or type command to continue
euler37.cpp: In function ‘std::vector<long long int> primesUpto(long long int)’:
euler37.cpp:16:26: warning: comparison between signed and unsigned integer expressions [-Wsign-compare]
 for (size_t i = 0; i < sqrtLimit; ++i)
 ^
euler37.cpp:20:40: warning: comparison between signed and unsigned integer expressions [-Wsign-compare]
 for (size_t j = (2 * i); j < limit; j += i)
 ^
euler37.cpp: In function ‘bool isTruncPrime(long long int, const std::vector<long long int>&)’:
euler37.cpp:39:23: warning: comparison between signed and unsigned integer expressions [-Wsign-compare]
 for (int i = 1; i < numberString.size(); ++i)

You just need to use unsigned integers.

Algorithm

You could use a different approach to solve this problem. Instead of trying to find an efficient way to check if a number is a truncatable prime, you could try to find a way to generate them. You could use the fact that the all right-truncable-primes of size n can be built easily from the list of right-truncable-primes of size n-1. The same logic applies to left-truncable-primes.

Here is a code example based on Toby Speight's code:

template<std::size_t N>
constexpr bool isLeftTruncPrime(std::size_t n, const std::array<bool,N>& primes)
{
 for (std::size_t x = 10; x < n; x*= 10)
 if (!primes[n%x]) return false;
 return true;
}
template<std::size_t N>
constexpr bool isRightTruncPrime(std::size_t n, const std::array<bool,N>& primes)
{
 for (; n; n/=10)
 if (!primes[n]) return false;
 return true;
}
template<std::size_t N>
constexpr bool isTruncPrime(std::size_t n, const std::array<bool,N>& primes)
{
 return isRightTruncPrime(n, primes) && isLeftTruncPrime(n, primes);
}
#include <iostream>
int main()
{
 const auto limit = 1'000'000;
 static const auto primes = primesUpto<limit>();
 auto expected_sum = 748317ull;
 if (1)
 {
 auto sum = 0ull;
 // Start at 11, because single-digit numbers aren't truncatable
 // Advance by 2, because even numbers are non-prime (except 2, but 2<11)
 for (std::size_t n = 11; n < primes.size(); n += 2) {
 if (isTruncPrime(n, primes))
 sum += n;
 }
 std::cout << sum << "\n";
 std::cout << (sum == expected_sum) << "\n";
 }
 if (0)
 {
 auto sum = 0ull;
 std::vector<int> truncprimes;
 truncprimes.push_back(2);
 truncprimes.push_back(3);
 truncprimes.push_back(5);
 truncprimes.push_back(7);
 while (!truncprimes.empty())
 {
 std::vector<int> truncprimes2;
 for (auto tp: truncprimes)
 {
 auto tens = tp*10;
 for (size_t d = 1; d < 10; d+=2)
 {
 auto n = tens + d;
 if (n < limit && primes[n])
 {
 truncprimes2.push_back(n);
 if (isLeftTruncPrime(n, primes) /*&& n > 10 */)
 sum += n;
 }
 }
 }
 truncprimes = truncprimes2;
 }
 std::cout << sum << "\n";
 std::cout << (sum == expected_sum) << "\n";
 }
}

Then, because of the small amount of prime check we perform, pre-computing a sieve is not worth the pain. Thus, you don't even need to assume a maximum value. You don't even need to assume there will be 11 numbers found.

The new code is way faster than the original code:

bool isPrime(std::size_t num)
{
 if(num <= 1) return false; // <- 1 is not prime 
 for (std::size_t i = 2; i * i <= num; i++)
 if (num % i == 0) return false;
 return true;
}
bool isLeftTruncPrime(std::size_t n)
{
 for (std::size_t x = 10; x < n; x*= 10)
 if (!isPrime(n%x)) return false;
 return true;
}
#include <iostream>
int main()
{
 auto expected_sum = 748317ull;
 if (1)
 {
 auto sum = 0ull;
 std::vector<int> truncprimes;
 truncprimes.push_back(2);
 truncprimes.push_back(3);
 truncprimes.push_back(5);
 truncprimes.push_back(7);
 while (!truncprimes.empty())
 {
 std::cout << truncprimes.size() << "\n";
 std::vector<int> truncprimes2;
 for (auto tp: truncprimes)
 {
 auto tens = tp*10;
 for (size_t d = 1; d < 10; d+=2)
 {
 auto n = tens + d;
 if (isPrime(n))
 {
 truncprimes2.push_back(n);
 if (isLeftTruncPrime(n))
 sum += n;
 }
 }
 }
 truncprimes = truncprimes2;
 }
 std::cout << sum << "\n";
 std::cout << (sum == expected_sum) << "\n";
 }
}

Conclusion

I really like this problem and this solution because it shows that going straight into the Eratosthene sieve is not always the perfect solution.

• Arnav Borborah: Wow, this is also a great solution! I like the suggested alternative. However, when I compare your code to Toby's, the performance is about the same. (Yours is 0.07ms, while his comes to 0.08ms). The command I used is g++ Problem037.cpp -std=c++14 -O2. In the end, I will probably use Toby's solution to implement this, since that is easier to understand for me in general. Thanks for pointing out the signed/unsigned mismatch!

Answer 3

Cutting the number of divisions in the prime check in half cuts the running time roughly in half too (measured with time using clang 5 w/ -O3 anyway):

bool isPrime(std::size_t num) {
 if (num <= 1) return false; // <- 1 is not prime
 if (2 == num) return true;
 if (0 == (num & 1)) return false;
 for (std::size_t i = 3; i * i <= num; i += 2)
 if (num % i == 0) return false;
 return true;
}

Perhaps not as important in this case, but in general it is usually a good idea to avoid going to the heap too often, so instead of copying truncprimes2, move it outside the loop and then:

std::swap(truncprimes, truncprimes2);
truncprimes2.clear();

That way the heap should only get involved when the vectors need to grow.

For measuring this sort of thing, a dedicated microbenchmarking library is going to give more accurate results than relying on time. For example, this one works well: https://github.com/google/benchmark

With a recent enough C++ standard, the truncprimes initialization can be shorter:

std::vector<size_t> truncprimes{ 2, 3, 5, 7 };

EDIT: There are a few other things that can help performance.

Adding "-march=native" to the compiler flags doubled performance (with clang 5 on a 1st gen Xeon E3-1240).

Using 32-bit integers instead of 64-bit integers tripled performance (VS2017 on a 5 year old laptop).

Using a more complicated prime check can help too. With VS2017 on the same laptop as above, this is about 3x faster than the "int" version of the above isPrime() and even 20% faster than looking up primes in a pre-sieved std::vector<bool>; I didn't try std::array<bool, N>. Presumably this has more to do with memory latency on this laptop than from the sieve lookup polluting the data cache, but that is just speculation on my part.

static bool check_by_30(int num) {
 if (num <= 1) return false;
 if (2 == num) return true;
 if (0 == (num & 1)) return false;
 if (3 == num) return true;
 if (0 == (num % 3)) return false;
 if (5 == num) return true;
 if (0 == (num % 5)) return false;
 if (7 == num) return true;
 if (0 == (num % 7)) return false;
 if (11 == num) return true;
 if (0 == (num % 11)) return false;
 if (13 == num) return true;
 if (0 == (num % 13)) return false;
 if (17 == num) return true;
 if (0 == (num % 17)) return false;
 if (19 == num) return true;
 if (0 == (num % 19)) return false;
 if (23 == num) return true;
 if (0 == (num % 23)) return false;
 if (29 == num) return true;
 if (0 == (num % 29)) return false;
 // 30 = 2 * 3 * 5 
 for (int i = 30; i * i <= num; i += 30)
 {
 if (num % (i + 1) == 0) return false;
 if (num % (i + 7) == 0) return false;
 if (num % (i + 11) == 0) return false;
 if (num % (i + 13) == 0) return false;
 if (num % (i + 17) == 0) return false;
 if (num % (i + 19) == 0) return false;
 if (num % (i + 23) == 0) return false;
 if (num % (i + 29) == 0) return false;
 }
 return true;
}

A somewhat simpler variant is about as fast as a sieve lookup:

static bool check_by_6(int num)
{
 if (num <= 1) return false;
 if (2 == num) return true;
 if (0 == (num & 1)) return false;
 if (3 == num) return true;
 if (0 == (num % 3)) return false;
 if (5 == num) return true;
 if (0 == (num % 5)) return false;
 for (int i = 6; i * i <= num; i += 6)
 {
 if (num % (i + 1) == 0) return false;
 if (num % (i + 5) == 0) return false;
 }
 return true;
}

Note that computing all the primes in a given range as fast as possible is not the same as trying to determine if a random small number is prime. The distribution of numbers passed to the prime check will make a difference too (e.g., uniformly random vs. incrementing).

• SylvainD: Hi, these are good comments! Just to be sure, they apply to my answer, not too the original code?
• Arnav Borborah: @Josay, you are right, these changes make your code much faster. Now when I ran it I got (using time), I got 0.04ms. Of course, using a microbenchmarking library will obviously be more accurate, but it really is unneeded in this case. (I just wanted to see relative performance)
• Henric Jungheim: Note that even when the numbers involved are too big for 32-bit integers, the performance difference may be enough to warrant having the 64-bit version calling a 32-bit version at runtime whenever the argument is small enough.
• Henric Jungheim: @Josay, yes, to your changes.