Project Euler #12 - highly divisible triangular number

Question

The following is a solution to Project Euler problem number 12: which is to find the first triangular number with more than 500 divisors.

The code takes far too long to complete (over 10 minutes). Is there any way to optimize the run-time and to make it more elegant?

#include <iostream>
#include <chrono>
#include <ctime>

inline void triangle(int &num, int &interval)
{
    num = ((interval * interval) + interval ) / 2;
}

int main()
{
    long int input;
    int num = 1, interval = 1, divisor = 1, no_divisor = 0;

    std::cout << "Enter factor limit for Triangle number: ";
    std::cin >> input;

    auto start = std::chrono::steady_clock::now();

    for (;;)
    {
        while (divisor <= num)
        {
            if (0 == num % divisor) //check factor;
            {
                std::cout << "factor: " << divisor << "\n";
                no_divisor++;
                if (no_divisor > input) //divisor count over input number, jump to exit;
                {
                    goto exit;
                }
            }
            divisor++;
        }
        std::cout << "\n";
        interval++;
        triangle(num, interval);
        divisor = 1; //reset divisor back to 1;
        no_divisor = 0; //reset divisor count back to 0;
    }
    exit: std::cout << "the first Triangle number to have over " << input << " divisors : " << num << "\n";

    auto end = std::chrono::steady_clock::now();
    auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);

    std::cout << "\nIt took me " << elapsed.count() << " milliseconds." << "\n";

    return 0;
}

Optimization I've done:

1. using std::cout << "\n" instead of std::cout << endl to prevent flushing output buffers. Source: this
2. using inline void triangle(int &num, int &interval) speed up 10% time
3. using n(n+1)/2 according to this

Speculation:

I figured if the question asked for over 500 factors, which for example: Triangular number 28 has over 5 factors:

So I only need to iterate up to half of the limit. But how can I implement that into my code to only iterate up to half as many factor?

Answer 1

Project Euler problems typically don't benefit from micro-optimizations such as inlining functions or flushing standard output. The biggest speed gains often come from algorithmic improvements. Here there are two key mathematical facts that can make this computation tractable:

1. A number with prime factorisation n = Product_i p_i ^ e_i (where p_i are primes and e_i their exponents) has a number of divisors equal to d = Product_i (e_i + 1).
2. A triangle number is equal to n * (n + 1) / 2. These factors have no prime factors in common and only one of them has a factor of two. This means that the number of divisors of a triangle number can be factored in the number of the divisors of n and n+1.

In order to code this algorithm, a few more guidelines to writing better code:

1. write functions that take an int and return an int. Don't use Pascal like procedures by modifying a reference argument. This makes your code harder to reason about.
2. separate computation from I/O: an algorithm should ideally have no side effects, and be silent. If you want to query the result, just print that.

Below some working code:

#include <chrono>
#include <ctime>
#include <iostream>

int triangle(int n)
{
    return ((n + 1) * n) / 2;    
}

// Math fact: a prime decomposition n = Prod_i p_i^e_i yields
// Prod_i (e_i + 1) different divisors
int num_divisors(int n)
{
    auto divisors = 1;  // 1 has one divisor
    {
        auto f = 2;
        auto count = 0;
        while (n % f == 0) {
            ++count;
            n /= f;     // divide by factor of 2
        }
        divisors *= (count + 1);
    }

    // <= n guarantees that a prime has 2 divisors
    for (auto f = 3; f <= n; f += 2) {
        auto count = 0;
        while (n % f == 0) {
            ++count;
            n /= f;     // divide by odd factor
        }
        divisors *= (count + 1);
    }
    return divisors;
}

int unique_divisors(int n)
{
    if (n % 2 == 0)
        n /= 2;
    return num_divisors(n);
}

// triangle(n) = (n + 1) * n / 2
// only one of them has a factor of two, and they have no prime factors in common
// we can cache one of the divisor counts while iterating over triangle numbers
int first_triangle_with_more_than_n_divisors(int n)
{
    auto i = 1;
    for (auto f1 = unique_divisors(i), f2 = unique_divisors(i + 1); f1 * f2 <= n; ++i)   {
        f1 = f2;                        // re-use result from previous iteration
        f2 = unique_divisors(i + 2);    // only one new computation
    }
    return i;
}

int main()
{
    auto start = std::chrono::steady_clock::now();

    auto const n = 500;
    auto res = first_triangle_with_more_than_n_divisors(n);

    std::cout << "the first Triangle number to have over " << n << " divisors : " << triangle(res) << "\n";

    auto stop = std::chrono::steady_clock::now();
    auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(stop - start);

    std::cout << "\nIt took me " << elapsed.count() << " milliseconds." << "\n";
}

Live Example that runs in less than 0.2 seconds.

Answer 2

If you'd like it even faster consider this:

#include <math.h>
#include <iostream>
#include <chrono>

using namespace std;

int TR_triangle(int n)
{
    return ((n + 1) * n) / 2;
}

// Math fact: a prime decomposition n = Prod_i p_i^e_i yields
// Prod_i (e_i + 1) different divisors
int TR_num_divisors(int n)
{
    auto divisors = 1;  // 1 has one divisor
    {
        auto f = 2;
        auto count = 0;
        while (n % f == 0) {
            ++count;
            n /= f;     // divide by factor of 2
        }
        divisors *= (count + 1);
    }

    // <= n guarantees that a prime has 2 divisors
    for (auto f = 3; f <= n; f += 2) {
        auto count = 0;
        while (n % f == 0) {
            ++count;
            n /= f;     // divide by odd factor
        }
        divisors *= (count + 1);
    }
    return divisors;
}

int TR_unique_divisors(int n)
{
    if (n % 2 == 0)
        n /= 2;
    return TR_num_divisors(n);
}

// triangle(n) = (n + 1) * n / 2
// only one of them has a factor of two, and they have no prime factors in common
// we can cache one of the divisor counts while iterating over triangle numbers
int TR_first_triangle_with_more_than_n_divisors(int n)
{
    auto i = 1;
    for (auto f1 = TR_unique_divisors(i), f2 = TR_unique_divisors(i + 1); f1 * f2 <= n; ++i)   {
        f1 = f2;                        // re-use result from previous iteration
        f2 = TR_unique_divisors(i + 2);    // only one new computation
    }
    return i;
}
int* TIN_ESieve(int upperLimit)
{
    int sieveBound = (int)(upperLimit - 1) / 2;
    int upperSqrt = (int)((sqrt((double)upperLimit) - 1) / 2);
    bool PrimeBits[1000] = {false};
    int* numbers = new int[1000];
    int index = 0;
    for(int i = 2; i <= upperSqrt; i++)
    {
        if(!PrimeBits[i])
        {
            for(int j = i * i; j <= sieveBound; j += i)
            {
                PrimeBits[j] = true;
            }
            numbers[index++] = i;
        }
    }
    for(int i = upperSqrt + 1; i <= sieveBound; i++)
    {
        if(!PrimeBits[i])
        {
            numbers[index++] = i;
        }
    }
    numbers[index] = 0;
    return numbers;
}
int TIN_PrimeFactorisationNoD(int number, int primelist[])
{
    int nod = 1;
    int exponent;
    int remain = number;

    for(int i = 0; primelist[i] != 0; i++)
    {

        // In case there is a remainder this is a prime factor as well
        // The exponent of that factor is 1
        if(primelist[i] * primelist[i] > number)
        {
            return nod * 2;
        }

        exponent = 1;
        while(remain % primelist[i] == 0)
        {
            exponent++;
            remain = remain / primelist[i];
        }
        nod *= exponent;

        //If there is no remainder, return the count
        if(remain == 1)
        {
            return nod;
        }
    }
    return nod;
}
int main()
{
    auto start = chrono::steady_clock::now();
    for(int c1 = 1; c1 <=10; c1++)
    {
        auto sieve = TIN_ESieve(1000);
        int i1 = 2;
        int cnt = 0;
        int Dn1 = 2;
        int Dn = 2;
        while(cnt < 500)
        {
            if(i1 % 2 == 0)
            {
                Dn = TIN_PrimeFactorisationNoD(i1 + 1, sieve);
            }
            else
            {
                Dn1 = TIN_PrimeFactorisationNoD((i1 + 1) / 2, sieve);
            }
            cnt = Dn * Dn1;
            i1++;
        }
        if(c1 == 1)
            cout << "the first Triangle number to have over " << 500 << " divisors : " << i1 * (i1 - 1) / 2 << "\n";
    }
    auto stop = chrono::steady_clock::now();
    auto elapsed = chrono::duration_cast<std::chrono::milliseconds>((stop - start)/10);

    cout << "\nUsing TIN methods it averaged " << elapsed.count() << " milliseconds." << "\n\n";
    start = chrono::steady_clock::now();
    for(int c2 = 1; c2 <= 10; c2++)
    {
        auto const n = 500;
        auto res = TR_triangle(TR_first_triangle_with_more_than_n_divisors(n));
        if(c2 == 1)
            std::cout << "the first Triangle number to have over " << n << " divisors : " << res << "\n";
    }
    stop = chrono::steady_clock::now();
    elapsed = chrono::duration_cast<std::chrono::milliseconds>((stop - start) /10);

    cout << "\nUsing TR methods it averaged " << elapsed.count() << " milliseconds." << "\n";
    return 0;
}

This Live Example compares the 2 with a significant difference in execution times.

This was ported from a C# method I found. I can't remember where though. My apologies to the originator.

The sieve, based on the Sieve of Eratosthenes, basically uses an array of booleans initialized to false. You start with the first prime,2, and change every successive boolean that's a multiple to true. Repeat with every element that is false when the loop gets there. When the loop is done, the index of every element that is false is a prime. What I've done is added an optimization and added the primes to the list in the main loop then used another loop that starts where that one left off. The trick is to use booleans so that you're only checking true/false instead of decoding an integer.