Project Euler 12
March 29, 2019
When a question about Project Euler 12 came up today at Reddit, I came here to link my solution to the problem, only to find out that we have never done that problem. So today we will.
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,28We 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?
Your task is to write a program to solve Project Euler 12. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
I did many of the Euler problems some years ago.
I posted the code on ideone.
A naive divisor counting approach works better if we make use of the fact that triangular numbers are of the form n*(n+1)/2: so just need to count the divisors of n and 2n+1 or 2n-1:
def ndivisors(n): return sum(1 for i in range(1,n+1) if n%i == 0) n,last = 1,0 while True: count = ndivisors(n//2 if n%2 == 0 else n) if count*last > 500: break n,last = n+1,count print(n,n*(n-1)//2,count*last)Takes just a few seconds to get to 500 divisors.
A Haskell version.
A Julia version…
function NumberOfDivisors(n::Int64)
m = n
c = 2
d = n / 2
d_ = round(Int64, d)
end
function main(n::Int64 = 1000000)
x = 0
nd = 0
end
The script can be remedied so as to cover larger integers (BigInt type), if needed. Although somewhat slower in that case, it is still fast enough to be practical.
Another Algorithm:
We need to calculate the number of divisors of n(n+1)/2. This number is the product of two
numbers, one of which is even. Divide the even one by 2. Now we find the number of
divisors of the product n1 * n2. This can be done by finding the factors of n1 and n2, merging
them, then calculating the number of divisors from the number of each factor. But it should
be noted that as we scan through the values of n, we may have already calculated the factors
of the smaller of n1 and n2 and sometimes even the factors of the larger. Therefore we maintain
a list of the factors of each n we have seen.
This means we shall factor roughly the same number of n’s that we calculated for n(n+1)/2 but
each one we calculate is smaller (roughly = sqrt(n(n+1)/2) ) than the product n(n+1)/2. Thus it should be faster than
calculating the divisors of n(n+1)/2. On my computer it is 7X faster (21 msec).
in C#
static void Main(string[] args)
{
Stopwatch s = new Stopwatch();
int N = 15000;
List[] faclist = new List[15000];
s.Start();
for(int i = 0; i < N; i++)
faclist[i] = new List();
for(int i = 2; i < N; i++)
{
int n1 = i;
int n2 = i+1;
if((n1 & 1) == 0) //even?
n1 /= 2;
else
n2 /= 2;
if(faclist[n1].Count == 0)
faclist[n1] = factor(n1);
if(faclist[n2].Count == 0)
faclist[n2] = factor(n2);
if(numdivisors(merge(faclist[n1], faclist[n2])) > 500)
{
Console.WriteLine(“{0} {1}”, i, n1 * n2);
break;
}
}
s.Stop();
Console.WriteLine(s.Elapsed);
} //end Main
@Ernie: note that the two factors a,b of a triangular number t = ab = n(n+1)/2 are either n/2, n+1 or n, (n+1)/2 and since n and n+1 must be coprime, so must a and b. My solution above uses that fact as does this Rust program, which uses a sieve to efficiently find the number of divisors of 2..N (I’m just starting with Rust so may not be stylistically optimal, improvements welcome):
fn main() { const N : usize = 10000000; let mut a = vec![1;N]; // build divisor count here let mut b = vec![0;N]; // stores power count for each prime for p in 2..N { if a[p] == 1 { // p is prime let mut pp = p; // powers of p while pp < N { // Increment power count for multiples of // powers of p let mut i = pp; while i < N { b[i] += 1; i += pp; } pp *= p; } // Number of factors is (n0+1)..(nk+1) // where ni = multiplicity of ith prime factor let mut i = p; while i < N { a[i] *= b[i]+1; b[i] = 0; i += p; } } } // a[i] == ndivisors(i) let mut tmax = -1; // n and n+1 are coprime, therefore so are n/2 and n+1 (or n, (n+1)/2), // and for coprime a and b, ndivisors(a*b) = ndivisors(a)*ndivisors(b) // so can use that to find ndivisors of n*(n+1)/2, a triangular number for n in 1..N-1 { let (n0,n1) = if n%2 == 0 { (n/2,n+1) } else { (n,(n+1)/2) }; let t = a[n0]*a[n1]; let s = n0*n1; if t > tmax { println!("{} {} {}",n,s,t); tmax = t; } } }Looking at the first 1000000 triangular numbers tells us eg. that the 7289919th triangular number, 26571463158240, has 7680 factors, which as it’s a record breaker, appears in https://oeis.org/A076711, “Highly Composite Triangular Numbers” (at position 43: https://oeis.org/A076711/b076711.txt).
Thank you Matthew for that observation. I only need to change one line of code (eliminate the merge) and the programme becomes 40% faster.
Here’s a brute-force solution in C. It takes 0.270 seconds to run on my desktop (with or without compiler optimization flags).
#include <stdio.h> #include <stdint.h> #include <stdlib.h> typedef uint64_t u64; int main(void) { u64 x = 0; for (u64 i = 0; ; ++i) { x += i; u64 n = 0; for (u64 j = 1; j * j <= x; ++j) { if (x % j == 0) { ++n; if (j * j != x) ++n; } if (n > 500) { printf("%lu\n", x); return EXIT_SUCCESS; } } } }Mumps version: Converted Matthew’s python to Mumps. Thanks, @Matthew.
@ProgrammingPraxis: My apologies – I meant to blank out the middle answer, and forgot to do so prior to posting my answer.