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:

enter image description here

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?

• ratchet freak: you only need to factor n and n+1 (first divide the even one by 2); they are coprime and thus have mutually exclusive factors so nd(n*(n+1))=nd(n)*nd(n+1)
• Ali Caglayan: Here is a tip. Project Euler problems are designed so that they can run under a second (usually, there is always an exception). Analysis of the problem can lead to a better algorithm. Some people (who are extremely talented and/or have too much time) can do project Euler with pen and paper.

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.

• vincentChen: is there any reason you use autoinstead of int ? IIRC, variable live from the start of their definition until the end of their block. Putting auto in front of them is redundant since that is the default anyway. Please do correct me.
• TemplateRex: this is a C++11 feature: auto means that the type will be deduced from the initalizer. E.g. auto i = 1; will give i the type of int. The old auto keyword from C has been deprecated. If you cannot use C++11, just replace auto with int everywhere in the above code.
• vincentChen: Thanks for clarifying. Another question: how is separate functions is faster over nested loop, in this case mine. I see you called in between functions. If my question has been asked before, please cite the source here. Thanks!
• TemplateRex: @vincentChen my advice is to first code for clarity, then for speed. A separate function that returns by value instead of modifying a reference, makes the code easier to reason about. It is also easier to test a small function that does only thing. You can use such small functions as building blocks for larger algorithms. You shouldn't worry too much about function call overhead. The compiler is likely to inline it anyway when given a suitable optimization flag.

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.

• TemplateRex: Upvoted, such algorithmic optimizations are what Project Euler is all about.
• vincentChen: very complicated, I wrapped my head around for a day now, still can't understand how you averaging the result. And how the TIN sieve works? what is the mathematical concept behind it? Thanks for answering my question.
• user33306: I added a simple explanation for you.