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?

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";
}

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.

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).