[Exercise] Triangle numbers divisors

Joe

This exercise is taken from a famous online programming exercise website. It's really interesting and worth adding to our collection.
Note that your code should execute reasonably fast, meaning you should get the solution in approx. 1 min on a modern computer.

The sequence of triangle numbers is generated by adding the natural numbers. So the 7th 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?

Ayman

This is a very nice exercise, I worked on it a couple of days ago here is the code in C#, it can be a bit more concise but it does the job.

        static void Main(string[] args)
        {
            int triangle = 1;
            int triangleValue = 1;
            int divisors = 0;
       
            while (divisors < 500)
            {
                Console.WriteLine("Triangle: {0} , Value: {1}, Divisors: {2}", triangle, triangleValue,divisors);

                divisors = 0; 

                triangleValue = getTriangleVal(triangle);

                divisors = getNumDivisors(triangleValue);
                triangle++;
            }
            
            Console.WriteLine("{0} divisors for {1}",divisors,triangleValue);
        }

        public static int getTriangleVal(int triangle)
        {
            int sum = 0;
            for (int i = 1; i <= triangle; i++)
                sum += i;
            return sum;
        }

        public static int getNumDivisors(int n)
        {
            List<int> rawPrime = getPrimeFactors(n);
            Dictionary<int, int> factorCount = new Dictionary<int, int>();
            int product = 1;

            foreach(var factor in rawPrime)
            {
                int count = rawPrime.Where(temp => temp.Equals(factor))
                    .Select(temp => temp)
                    .Count();

                if (!factorCount.ContainsKey(factor))
                    factorCount.Add(factor,count);
            }

            foreach (var factor in factorCount)
                product = product * (factor.Value + 1);

            return product;
        }

        public static List<int> getPrimeFactors(int n)
        {
            List<int> rawFactors = new List<int>();

            int prime = 2;

            while (prime <= n)
            {
                if (n % prime == 0)
                {
                    n = n / prime;
                    rawFactors.Add(prime);
                }
                else prime++;
            }
            return rawFactors;
        }

The program execution time is 3 seconds. Since we only need the number of divisors for the number no need to check each and everyone. We can do prime factorization on the triangle number, get the powers of the prime factors, add one to each, multiply them and we'll get the number of divisors.

m0ei

Nice one, same as Ayman I've done it coupe of days ago when solving Project Euler exercises.

Here's my solution in Python:

counter = 0
n = 1

while counter < 500:
    counter = 2
    n += 1

    triangle = n*(n+1)/2
    for i in xrange(2, int(triangle**0.5)):
        if triangle % i == 0:
            counter += 2

    print triangle

Ra8

I love Project Euler!
Add me if you want:

91920865310012_30465243fa220e4397d3dae974d3f4bf

My solution in C++ in 0.248 seconds

#include <iostream>
#include <cmath>

using namespace std;
long tri(int n)
{
    return n*(n-1)/2;
}
int divs(long n)
{
    int tot=2;
    int div=2;
    double sq=sqrt(n);
    while(div<sq)
    {
        if(n%div==0)tot+=2;
        div++;
    }
    if((long)(sq)==sq)tot++;
    return tot;
}
int main()
{
    long n;
    int i=0;
    do
    {
        i++;
        n=tri(i);
    }while(divs(n)<500);
    cout << n <<endl;
    return 0;
}

divs(n) counts the number of divisors:
i start with 2 because 1 and n are divisors
then if a number a is a divisor, n/a is also a divisor of n
i stop at the square root of n, if the square root is an integer add 1 to the total

Edit: Even Faster: 0.027 seconds

#include <iostream>
#include <cmath>

using namespace std;
long divs(long n)
{
    long long tot=2;
    long long div=2;
    double sq=sqrt(n);
    while(div<sq)
    {
        if(n%div==0)tot+=2;
        div++;
    }
    if((long long)(sq)==sq)tot++;
    return tot;
}
int main()
{
    long long n1,n2;
    int i=0;
    do
    {
        i++;
        if(i%2==0)
        {
            n1=(i+1);
            n2=i/2;
        }
        else
        {
            n1=(i+1)/2;
            n2=i;
        }
    }while(divs(n1)*divs(n2)<500);
    cout << (n1*n2) <<endl;
    return 0;
}

the triangle number is n(n+1)/2
i divide the even number between n and n+1 and then multiply the number of divisor of each to get the total.

for >1000 divisors: 842161320 0.16 seconds
for >2000 divisors: 49172323200 2.97 seconds

NuclearVision

Hey guys, i wrote a solution in python, it is very logic, however it does not seem to work. Any ideas?

def divlist(n):
    list = []
    for i in xrange(1,int((n/2)+1)):
        if n%i==0:
            list.append(i)
    list.append(n)
    return list
def T(n):
    return n*(n+1)/2
n=0
while True:
    n+=1
    if len(divlist(T(n)))>500:
        print T(n)
        break

Joe

@NuclearVision:

Your code is fine, it's just extremely slow. If you test it on low thresholds like 5 or 50 (instead of 500) you'll get correct results.

A big part of the exercise is to try to optimize the code. Take a look at the 3 pieces of code above you, as well as the explanations given by Ayman and Ra8. It should give you ideas on how to optimize.

Joe

Here's my solution in Python:

from functools import reduce
from collections import Counter
from operator import mul

def prime_factors(n):
    d = 2
    while (n > 1):
        while (n%d==0):
            yield d
            n /= d
        d = d + 1

def format_prime(n):
    return Counter(prime_factors(n))

def prime_count(c):
    return reduce(mul, (c[i] + 1 for i in c))

def triangle(i):
    a = format_prime(i)
    b = format_prime(i+1)
    c = a + b

    c[2] -= 1  # This is equivalent to a division by 2

    return c


def main():

    i = 2
    while prime_count(triangle(i)) <= 500:
        i += 1
    print " ", i, triangle(i), i * (i+1) / 2

if __name__ == "__main__":
    main()

I had originally commented the code, but it quickly got out of hand. And this led me to attempt my first literate programming essay.

In this article I explain the logic behind my code. I would love to get as much feedback and suggestions as possible.

For the record the code shares some ideas already exposed here, like:

- the number of divisors equals the product of exponents of prime factor decomposition plus 1.
- triangle(i) = i * (i+1) / 2

arithma

@Ra8: That trick you have is pretty neat. How'd you think about it? (Referring to the }while(divs(n1)*divs(n2)<500);)

//
//  main.cpp
//

#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <cmath>
using namespace std;

vector<int> primes;

class PrimeItr
{
    vector<int>::iterator itr;
public:
    PrimeItr(){
        itr = primes.begin();
    }
    
    int operator*(){
        return *itr;
    }
    
    PrimeItr & operator++(){
        if(itr+1 == primes.end()){
            int n = *itr + 1;
            for(;;){
                bool divisible = false;
                for(vector<int>::iterator itr = primes.begin(); itr != primes.end(); itr++){
                    if(n % (*itr) == 0) {
                        divisible = true;
                        break;
                    }
                }
                
                if(!divisible){
                    primes.push_back(n);
                    itr = primes.end()-1;
                    return *this;
                }
                
                n++;
            }
        }
        else {
            itr++;
            return *this;
        }
    }
};

int primePower(int64_t& n, int prime){
    int count = 0;
    
    while(n % prime == 0){
        
        int64_t multi = prime;
        int c = 1;
        while(n % (multi*multi) == 0){
            multi *= multi;
            c *= 2;
        }
        
        count += c;
        n /= multi;
    }
    
    return count;
}

int divs(int64_t n){
    int result = 1;
    
    // get prime/count
    PrimeItr itr;
    while(n > 1){
        int prime = *itr;
        if(n % prime == 0){
            result *= primePower(n, prime) + 1;
        }
        
        ++itr;
    }
    
    return result;
}

int main(){
    primes.push_back(2);
    
    int64_t sum = 0;
    int i = 1;
    int64_t d;
    do {
        sum += i;
        d = divs(sum);
        if(i%1000==0)
            cout << i << " triangle numbers tried so far" << endl;
        i++;
    } while (d < 1000);
    cout << "triangle# = " << sum << " & divisors = " << d << endl;
    cout << "largest prime " << *(primes.rbegin()) << endl;
    cout << "primes explored so far = " << primes.size() << endl;
    return 0;
}

Ra8

arithma wrote@Ra8: That trick you have is pretty neat. How'd you think about it? (Referring to the }while(divs(n1)*divs(n2)<500);)

Well this works only in this case since n and n+1 will never have a common divisor.
if you had for example 4 ={1,2,4} and 6={1,2,3,6} 24={1,2,3,4,6,8,12,24} in this case it will not work because 2 is a common divisor

NuclearVision

@rahmu can you please explain every byte in this?

def prime_factors(n):
d = 2
while (n > 1):
while (n%d==0):
yield d # Especially those
n /= d # two lines
d = d + 1

Thanks in advance rahmu, i just love learning programming.

arithma

Ra8 wrote
arithma wrote@Ra8: That trick you have is pretty neat. How'd you think about it? (Referring to the }while(divs(n1)*divs(n2)<500);)
Well this works only in this case since n and n+1 will never have a common divisor.
if you had for example 4 ={1,2,4} and 6={1,2,3,6} 24={1,2,3,4,6,8,12,24} in this case it will not work because 2 is a common divisor

The smart thing about it is notices how i & i+1 are coprime, and hence all their divisors are coprime as well.
Very observant is the least I can say.

I tried to beat its performance by keeping a set of already visited primes (while iterating through triangle numbers) but it' not helping at all.

#include <iostream>
#include <string>
#include <vector>
#include <set>
#include <algorithm>
#include <cmath>
using namespace std;

int count_reused = 0;
int count_explord = 0;
set<int> primes;

int primePower(int64_t& n, int prime){
    int count = 0;
    
    while(n % prime == 0){
        
        int64_t multi = prime;
        int c = 1;
        while(n % (multi*multi) == 0){
            multi *= multi;
            c *= 2;
        }
        
        count += c;
        n /= multi;
    }
    
    return count;
}

int divs(int64_t n){
    int result = 1;
    int sqr = sqrt(n);
    
    set<int>::iterator itr = lower_bound(primes.begin(), primes.end(), sqr);
    if(itr == primes.end())
        itr--;
    
    while(n > 1){
        int prime = *itr;
        if(n % prime == 0){
            count_reused++;
            result *= primePower(n, prime) + 1;
        }
        
        if(itr == primes.begin())
            break;
        
        itr--;
    }
    
    int div = 2;
    while(n > 1 && div < n){
        if(n % div == 0){
            count_explord++;
            primes.insert(div);
            result *= primePower(n, div) + 1;
        }
        
        div++;
    }
    
    return result;
}

int main(){
    primes.insert(2);
    
    int64_t sum = 0;
    int i = 1;
    int64_t d;
    do {
        sum += i;
//        d = divs(sum);
        d = i % 2 == 0 ? divs(i/2)*divs(i+1) : divs((i+1)/2)*divs(i);
        if(i%1000==0){
            cout << i << " triangle numbers tried so far" << endl;
            cout << sum << " << triangle number max" << endl;
        }
        i++;
    } while (d < 1000);
    cout << "reused " << count_reused << endl;
    cout << "explored " << count_explord << endl;
    cout << "triangle# = " << sum << " & divisors = " << d << endl;
    cout << "largest prime " << *(primes.rbegin()) << endl;
    cout << "primes explored so far = " << primes.size() << endl;
    return 0;
}

Analyzing through, this is what I came up with:
The first advantage for Ra8's approach is the fact that for triangle number n, it only need to find divisors for its square root. For odd i for example, n=i*(i+1)/2. i in this case is comparable to square root of n. The algorithm for finding the number of divisors for i only has to go up to sqrt(i].
So the complexity of finding divisors for a triangle number in this case is only: O(sqrt(sqrt(n))

I would have sworn before trying it that keeping aside prime numbers for later would even beat the above asymptotically, but it doesn't. It returns almost instantly for finding a triangle number with more than 500 divisors. But checking how it starts to behave near the larger numbers, it's no where near as performant. I have no idea why.

Joe

NuclearVision wrote@rahmu can you please explain every byte in this?

Sorry for the late answer. It behaves a little bit like this:

def pf(n):
    d = 2
    primes = []
    while (n>1):
        while (n%d == 0):
            primes.append(d)
            n /= d
        d = d + 1

    return primes

The difference between the two pieces of code are subtle and a bit difficult to grasp, but take a look anyway at this very good explanation of Python yield keyword.

NuclearVision

@rahmu i understood that yield creates a generator of the form [ func(j) for j in sequence ], that can be used once.
But i still don't get the n/=d., which probably means n=d/n, but we can't change the value of n, and let's say we changed it, is this some type of recursion?

Thanks for the explanation anyway.

Joe

NuclearVision wrotei still don't get the n/=d., which probably means n=d/n,

No.

In Python the two following expressions are equivalent:

• a (op)= b
• a = a (op) b

Thus n /= d is equivalent to n = n/d.

NuclearVision wrotebut we can't change the value of n

The function is not changing the value of n, it's changing the value of a copy of n. What do you think the output of the following code will be?

def foo(n):
    n = 5
    return n

n = 10
foo(n)
print(n)

Answer: 10.
This is because the variable n is never passed to function foo. Instead a copy of n, containing the same value, is sent.

Note that this behavior is specific to Python, and you cannot make the same assumptions about other languages. You can learn more by reading the Wikipedia article on Evaluation Strategy.

Also note that Python has a somewhat complex evaluation strategy (similar to the ones used by Java or - surprise - Ruby). You can read more about it here.

Understanding the code better
My favorite way of understanding a piece of code I cannot immediately understand, is to add print() calls. So I modified the code to be like this:

def pf(n):
    d = 2
    primes = []
    while (n > 1):
        while (n % d == 0):
            print d, n, primes
            primes.append(d)
            n /= d
        d = d + 1

    return primes

(Notice the print call on line 6)

You can now experiment with different input inside the python interpreter. I'm assuming you work on Windows, and I don't know how to launch the interpreter there, so maybe somebody else can help you with that. In any case, here are some sample outputs:

>>> pf(24)
2 24 []
2 12 [2]
2 6 [2, 2]
3 3 [2, 2, 2]
[2, 2, 2, 3]

>>> pf(38)
2 38 []
19 19 [2]
[2, 19]

>>> pf(236)
2 236 []
2 118 [2]
59 59 [2, 2]
[2, 2, 59]

longbit

so i did mine in java

public class TriangleNumber {
    public static boolean divNumber(int n){
        int counter=0;
        int sum=0;
        int end =(int) Math.sqrt(n);
        end++;//because looping to i<end is faster than looping to i<=end
        for(int i = 1 ; i < end; i++){
            if(n%i==0){
                counter+=2;
                sum+=i;
            }
        }
       
        if((counter>500)){
            System.out.println("the result is: "+n);
            return false;
        }

        else
            return true;
    }
    public static void main(String[] args){
       int val=0;
        int inc =0;
        boolean iterrate=true;
        while(iterrate){
            val+=inc;
            iterrate=divNumber(val);
            inc++;

        }

}
}

NuclearVision

ok here's what i got to, finally. Inspired from Ra8's

from math import sqrt
def t(n):
    return n*(n+1)/2
def div(n):
    j=0
    for i in xrange(1,int(sqrt(n))):
        if n%i==0:
            j+=1
    if sqrt(n)-int(sqrt(n))==0:
        return j*2 +1
    else:
        return j*2
i=1
while True:
    i+=1
    if div(t(i))>500:
        print t(i)
        break

CodingFreak

My solution in C++, takes 1/3 of a second to solve the exercise. Ra8's solution looks interesting, I'll check it out later.

#include <iostream>
#include <cmath>
using namespace std;

int number_of_divisors(int n) {

    int count = 2;
    int square_root = sqrt(n);

    for (int i = 2; i <= square_root; i++)
        if (n%i == 0)
            count++; 

    return (n%square_root == 0) ? count*2-1 : count*2;

}

int main() {

    long n = 1; 
    for (int i = 2; number_of_divisors(n) <= 500; i++)
        n += i;

    cout << n << endl;
    return 0;

}

NuclearVision

#include <iostream>
#include <math.h>
using namespace std;
long t(int n)
{  return n*(n+1)/2;
}
int divcount( long n)
{ int j=2,i=2,lim=sqrt(n); // the number it self and 1.
    for(i;i<lim;i++)
      { if (n%i==0)
	  {j+=2;
	  }
      }
  return j;
}
int main()
{ int n=1;
  while (true)
    { n++;
      if (divcount(t(n))>500)
	{ cout<<t(n);
          break;
	}
    }
  return 0;
}

m-sabra

public class triangle {
	public static void main(String [] args){
		int count=0;
		int triangle=0;
		int factors=1;
		int test;
		while (factors <= 500){
                        count=count+1;
                        triangle=triangle+count;
                        factors=1;
                        for(int a=1;a<(Math.sqrt(triangle));a++){
                                test=triangle%a;
				if ( test == 0){
                                   factors=factors+2;
				}
			}
		}
	System.out.print(+triangle);
	}
}

PatrickGhazal

Hey guys I'm new here so I don't really know how to post code the way you do, but anyway I've been trying to solve this for some time now and I can't figure out what's wrong, maybe one of you can help me :) here's my Java code :

public class Problem12
{
  public static void main(String[] args)
  {
    long i = 1L;
    long triangle = 0L;
    int nbrOfDiv = 0;
    int highestNbrOfDiv = 0;
    while (highestNbrOfDiv <= 500)
    {
      triangle = newTriangle(triangle, i);
      nbrOfDiv = divisors(triangle);
      if (highestNbrOfDiv < nbrOfDiv) highestNbrOfDiv = nbrOfDiv;      
      i++;
    }
    System.out.print(triangle);
  }
  
  public static long newTriangle(long triangle, long i)
  {
    return triangle+i;
  }
  
  public static int divisors(long number)
  {
    int nbrOfDiv = 0;
    for (long i = 1; i <= (long)(Math.sqrt(number))+1; i++)
    {
      if (number%i == 0) nbrOfDiv++;
    }
    return nbrOfDiv+1;
  }
}

I get 842161320 but this is apparently incorrect.

Would anyone happen to know what's wrong?

Thanks !