Problem 12 of Project Euler has a wording which is somewhat different than previous problems. However, as we shall see deriving efficient solutions for the problem, we can use theory which is very similar to some of the previous problems. The problem reads

The sequence of triangle numbers is generated by adding the natural numbers. So the 7

^{th}triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first ten terms would be:1, 3, 6, 10, 15, 21, 28, 36, 45, 55, …

Let us list the factors of the first seven triangle numbers:

1: 1

3: 1,3

6: 1,2,3,6

10: 1,2,5,10

15: 1,3,5,15

21: 1,3,7,21

28: 1,2,4,7,14,28

We can see that 28 is the first triangle number to have over five divisors.What is the value of the first triangle number to have over five hundred divisors?

When I first read the problem description, I rather quickly saw a direct approach to find the answer to the problem, where I use trial division to find the number of divisors to a number. This is the first approach I will describe in the this blog post. After that I will describe two modifications to solution strategy to speed up the algorithm.

## Brute Force with Trial Division

The trial division is pretty straight forward. The main piece of the code looks like

int number = 0; int i = 1; while(NumberOfDivisors(number) < 500){ number += i; i++; }

where the “intelligence” lies in the function it calls called NumberOfDivisors. The function uses trial division of all numbers up to and including the square root of the original number.

The code for that looks as

private int NumberOfDivisors(int number) { int nod = 0; int sqrt = (int) Math.Sqrt(number); for(int i = 1; i<= sqrt; i++){ if(number % i == 0){ nod += 2; } } //Correction if the number is a perfect square if (sqrt * sqrt == number) { nod--; } return nod; }

I don’t see much reason to go into the details of describing the code as it speaks for it self. The result is

The first triangle number with over 500 divisors is: 76576500 Solution took 214 ms

So while the code is easy to read, it isn’t terribly fast to execute.

## Prime Factorisation

Now that we have established a method to find the answer. The question is how to speed it up. It is pretty obvious that the time consuming part of the code is the NumberOfDivisors functions. So lets see if we can improve that.

The property we need in order to find an efficient method to derive the number of divisors of a number is based on the prime factorisation of the number which we also dealt with in Problem 3.

Any number N is uniquely described by its prime factorisation

Where p_{n} is the a distinct prime number, a_{n} is the exponent of the prime and *K* is the set of all prime numbers less than or equal to the square root of *N*. Taking any unique combination of the prime factors and multiplying them, will yield a unique divisor of *N*. That means we can use combinatorics to determine the number of divisors based on the prime factorisation. The prime *p _{n}* can be chosen 0,1,…,a

_{n}times. That means in all we can choose

*p*in

_{n}*a*+1 different ways.

_{n}The total number of divisors of *N* , *D(N)* can be expressed as

Now we just need a function to carry out the factorisation.

We used a prime factorisation in the solution to Problem 3, I have modified it to only use primes, instead of odd numbers, and I have modified it to actually count the divisors instead of returning a divisor.

The code looks like

private int PrimeFactorisationNoD(int number, int[] primelist) { int nod = 1; int exponent; int remain = number; for (int i = 0; i < primelist.Length; 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; }

It requires a list of primes, it starts from the lowest prime, and finds the exponent for that prime. Just as in Problem 3, if there the remainder reaches 1, there are no more divisors, and if there is a left over, that left over will be a prime divisor as well.

In the main loop we will need to generate a list of primes as well. This is done using the Sieve of Eratosthenes developed in Problem 10. With the line of code

int[] primelist = ESieve(75000);

I am uncertain of the upper limit of the prime divisors, so I tried to set it to 75000, to have a reasonably conservative bound.

Using prime factorisation yields the same result, but significantly faster than using trial division.

The first triangle number with over 500 divisors is: 76576500 Solution took 16 ms

## Coprime Property

The last improvement of the algorithm, that we will deal with today, is a decomposition of the problem, into two smaller problems, which are significantly easier to solve.

In problem 6, we established that the sum of natural numbers can be written analytically, such that

where *N _{k}* is the k’th triangle number.

The property that we will use is that k and k+1 are coprimes (which we also dealt with in Problem 9). You can convince your self of that in the general case of k, using Euclid’s alogrithm.

Since *k* and *k+1* are coprimes it means that their set of primes factors are distinct, and

D(N_{k}) = D(k/2)D(k+1) if k is even and

D(N_{k}) = D(k)D((k+1)/2) if k is odd.

In the first case, if *k* and *k+1* are coprimes, so is *k/2* and *k+1*, since 2 will be in the set of prime factors for *k*, and therefore not for *k+1*. Since otherwise they would not be coprimes. The problem of factorisation has now been split into prime factorisation for two smaller numbers, which is a significantly easier task.

Furthermore if we code it a bit smart, we can reuse the prime factorisation for k+1 in the subsequent iteration, and thus we only need to factorise one number in each iteration.

The prime generator and the prime factorisation we have already implemented can be reused, but we need to rework the main algorithm. I have chosen an implementations that looks like

int number = 1; int i = 2; int cnt = 0; int Dn1 = 2; int Dn = 2; int[] primelist = ESieve(1000); while (cnt < 500) { if (i % 2 == 0) { Dn = PrimeFactorisationNoD(i + 1, primelist); cnt = Dn * Dn1; } else { Dn1 = PrimeFactorisationNoD((i + 1) / 2, primelist); cnt = Dn*Dn1; } i++; } number = i * (i - 1) / 2;

Based on the explanation the code should be fairly self explanatory. One of the more obvious observations is that I have tightened the upper bound of the prime number generation signifanctly.

The execution of the algorithm yields

The first triangle number with over 500 divisors is: 76576500 Solution took 1 ms

For a number with 500 divisors, there wont be much difference, but I lost the patience of the second algorithm at 1200 divisors. Where this version still ran smoothly at more than 5000 divisors. So for scalability, the last version is a significant improvement.

## Wrapping up

I covered three iterations of the algorithm, each improving the execution of the algorithm significantly by using different areas of number theory and combinatorics. As usual the source code is available for download. Did it make sense?

Hi,

the code is awesome… hav a doubt though, for triangle numbers with over 6000 digits, we get negative results… is that correct??

Nazia

Hi Nazia

Thank you for the comment, it is much appreciated.

For all methods you will at some point end up with a negative result since the integer that stores the result can only hold a number less than 2147483647.

Depending on which of the methods you use you will meet this restriction sooner or later. I have checked the co-prime version of the code and it fails for 6000 divisors. Changing most of the integers to longs pushes the limits for when the method will fail. The limit for longs is 9223372036854775807. If that is too small a limit and if you want to have a piece of more scaling code (apart from the prime generation) you could convert the necessary variables to BigInteger as I use in the solution to problem 13

I hope this answers your current question, otherwise just ask again 🙂

/Kristian

Thanks for the quick reply Kristian…. i tried converting it to BigInteger, but came up with errors in the “ESieve()” method. Would highly appreciate if you could share code with the BigInteger version?

In the meanwhile, I converted to “ulong” and was able to generate results for triangle numbers with upto 20000 divisors, i hav set the upper bound of the prime number generation to “45” and it still works… do you think this would be safe for numbers with more than 1 lac divisors?

also can we have some calculation to determine the upper bound of the prime number generation based on the number of divisors that we are searching for ?

Hi Nazia

Let me give you a short answer now and then I will promise to get back to you with a longer answer including code.

I haven’t found an upper bound on the prime generation based on the number you want to test. However, using the last code you cannot check numbers larger than twice the fourth root of the bound of the prime number generation, so if you set the bound for prime generation so if you set it for 45, you can’t ensure that you have found all factors for numbers larger then 2 million. Setting the limit to 2.000.000 allows you to test numbers up to 8.0 × 10^24. How many divisors that bound can handle I don’t know.

Saying that, I don’t think the sieve method is the right approach if you want to generate large numbers, since it will have it’s limitation at primes around 2.7 million due to arrays being indexed by an integer. So I think we will need to chance the approach for prime generation.

If someone has an upper bound on the primes needed for a certain number of divisors, I would be very happy to hear about it.

I will get back with some code which is based on another type of prime generation, most likely the method I used in Problem 5.

/Kristian

thanks Kristian.. will be waiting…

Hi Nazia

I needed a bit more space and tools to solve it, so I couldn’t answer you in the comment. I made a blog post on the subject which you can find here Problem 12 – revisited.

There is source code for download in the blog post as well.

Hope that inspires you in the right direction

/Kristian

My program took about two hours to give me the answer

Can anybody explain why is this so

I have following code in c++ compiled with vc++ 10

run on 2.1 GHz dual processor

#include

#include

#include

using namespace std;

int factors(long int);

int remainder(long int j,long int num);

void main(void)

{

long int num=1;

int fact,n=2;

while(1)

{

fact=factors(num);

if(fact>500)

break;

num+=n;

n++;

}

cout<<"ans= "<<num<<endl;

_getch();

}

int factors(long int num)

{

long int j;

int n=0;

for(j=1;j<=(num/2);j++)

{

if(remainder(j,num))

n++;

}

return n;

}

int remainder(long int j,long int num)

{

int qout,rem;

qout=num/j;

rem=num-(qout*j);

if(rem==0)

return 1;

else

return 0;

}

Can I have your brain please?

Your algorithms are fantastic. My best one was comparable to your “Brute Force” one at best.

Thanks for the great blog!

Thanks for the compliment, and no you can’t have my brain since that would turn me into a zombie and then we would start the end of the world… Oh wait that is coming in a week anyway according to everyone.

I am just a beginner (19 yrs) at programming and am having a tough time figuring out algorithms , if you could help me with this one that would be awesome. What I fail to understand is why is the loop only till the sqrt of the number:

for(int i = 1; i<= sqrt; i++){

if(number % i == 0){

nod += 2;

}

}

and why is the no.of divisors being doubled here:

for (int i = 0; i number) {

return nod * 2;

}

To the first question:

For the number n, we have that every divisor to below the square root p, will have a corresponding divisor above the square root q, such p*q = n. Otherwise p would not be a divisor. Since we know that we don’t really need to check any numbers above the square root. We simply add two every time we find a divisor so we count both p and q.

To the second question:

In that algorithm I only count the divisors below the square root to begin with. Therefore I multiply by 2 in the end as of the first question. I hope this helps you on your way.

Hey, when you say like ‘ms 251’ you mean 251 milliseconds right?

Yes. The ms is a suffix though.

Sorry but could you please help me Mr.Kristian.

I’m 19 years old student in Thailand and I need to solve this problem in C++.

Is it possible to explain how to store and count the divisors of the value in c++.

Sorry for bad language though.

Nope, I can’t help you on that. It sounds like an excellent exercise for you to figure out on your own.

hello sir ,please tell me

how to make a triangle for only prime no

i am give the 27 ,all the prime no in

1

3 5

7 11 13

17 19 23 27

How long does the initial, naive implementation take?

214ms

/*The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be 1 + 2 + 3 + 4 + 500 + 6 + 7 = 28. The first ten terms would be: 1, 3, 6, 10, 1500, 21, 28, 36, 4500, 500500, …

Let us list the factors of the first seven triangle numbers:

1: 1

3: 1,3

6: 1,2,3,6

10: 1,2,500,10

1500: 1,3,500,1500

21: 1,3,7,21

28: 1,2,4,7,14,28

We can see that 28 is the first triangle number to have over five divisors.

What is the value of the first triangle number to have over five hundred divisors?*/

#include

#include

using namespace std;

int main()

{

vectorivec;

long long i=1,j=1,count=0,c=0,temp,n;

long long sum=0;

cout<>n;

while(true)

{

sum=sum+i;

//cout<<sum<<":";

for(j=1;j<=sum;j++)

{

c=sum%j;

if(c==0 && count<n)

{

// cout<<j<<",";

count++;

}

temp=count;

}

//cout<<"\n";

count=0;

if(temp==n)

{

cout<<sum;

break;

}

i++;

}

}

i can't find the sum when no of divisor is greater than 300

please reply

I don’t understand why nod+=2. I did something very similar but:

int nod = 0;

for(…)

{

for(…)

{

if(number%factor == 0) { num++; }

if(num == 500) { break; }

}

num = 0

}

The first for is just moving through the Triangular Numbers I already put in a generic list (triangleNumber[i]).

Thanks!

Table 1

Highly composite triangular numbers generated by a range of divisors (100-500)

Re: Project Euler 12: Triangle Number with 500 Divisors

http://img15.hostingpics.net/pics/283714pe12table1trinbr.jpg

Table 1 shows that the highly composite triangular numbers generated by a range of divisors (100-500) are all even numbers ending with at least one zero.

Common factors: 2x2x3x5x7.

Bigger highly composite triangular numbers:

OEIS list (5).

___

Sources

1) Kristian’s algorithm

http://www.mathblog.dk/files/euler/Problem12.cs

2) Prime Factorization Calculator

http://www.mathblog.dk/tools/prime-factorization/

3) Microsoft Visual C# 2010 Express

4) Microsoft Office Excel (2007)

5) Highly composite triangular numbers

http://oeis.org/A076711

And. ehm.. this is my code.. this works… but my hardware limits its array.. please check this.. if this code is true, i may go to problem 13…

struct listofprimes

{

int prime[100];

int top;

listofprimes()

{

prime[0]=2;

top=0;

}

};

static listofprimes pr;

static int historyof[100];

int where(int x)

{

int i=0;

while(pr.prime[i]*pr.prime[i] 2)

pr.prime[++pr.top]=x;

return 0;

}

int main(int argc, char* argv[])

{

historyof[1]=1;

int i=1,x=0;

while(x<=4)

{

int y=where(i+1);

if(y==0)

historyof[i+1]=2;

else

{

int a=i+1,b=0;

while(a%y==0)

{

a=a/y;

b++;

}

historyof[i+1]=historyof[(i+1)/y]/b*(b+1);

}

x=historyof[i]*historyof[i+1];

int n,m=0;

if(i%2==0)

n=i;

else

n=i+1;

while(n%2==0)

{

n=n/2;

m++;

}

x=x*m/(m+1);

cout<<i<<endl;

i++;

}

cout<<i*(i-1)/2;

getch();

return 0;

}

Table 1

Triangular numbers, square triangular numbers and digital roots.

http://img4.hostingpics.net/pics/256564pe12table1.jpg

Excerpt from ref.(1):

“It is interesting to note that digital root of all EVEN Square Triangular Numbers i.e. 36, 41616, 48024900, 55420693056 .. etc is always 9 and digital root of all ODD Square Triangular Numbers i.e. 1, 1225, 1413721, 1631432881, … etc is always 1. Also Square Triangular Numbers can never end in 2, 3, 4, 7, 8 or 9.”

http://www.shyamsundergupta.com/triangle.htm

___

Results in Table 1 shows that the predictions [1] do not

apply in many cases. We can also see that square triangular numbers can end in 4 or 9.

______

Sources:

1) Fascinating Triangular Numbers

http://www.shyamsundergupta.com/triangle.htm

2) Mathematische Basteleien

Triangular Numbers

http://www.mathematische-basteleien.de/triangularnumber.htm

This is brilliant. I made the mistake of trying to do this problem using the principles from the sieve or erastothenes and a humongous array. Thanks to this blogpost I’ve realized just what a horrible idea re-initializing a massive array to all 1’s at the start of every loop is.

The idea of only looking below square root and just adding two… and that’s just your “brute force” solution. Thank you.

here is my code this is giving correct output if i change cnt < 500 to ccnt < 100 but if i make it same it gives TLE can you tell what is wrong

#include

#include

#include

#include

#include

#include

using namespace std;

vector primelist = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997};

long getCount(long number){

long nod = 1;

long exponent;

long remain = number;

for (long i = 0; i number) {

return nod * 2;

}

exponent = 1;

while (remain % primelist[i] == 0) {

exponent++;

remain = remain / primelist[i];

}

nod *= exponent;

if (remain == 1) {

return nod;

}

}

return nod;

}

int a[1000003] = {0};

int main() {

long long num,cnt = 0;

int Dn1 = 2;

int Dn = 2;

int i = 2;

while (cnt > t;

while(t–){

long n;

cin >> n;

for(long i = 0 ; i n){

cout << i << "\n";

break;

}

}

}

return 0;

}

How come you do cnt = Dn * Dn1? Why isn’t it cnt = Dn or cnt = Dn1?

I just want to ask about the Brute force solution, why do you iterate the “nod” by +=2 ? Thanks.