Project Euler 12

March 29, 2019

Let’s start with a function to find the factors of a number n:

(define (factors n) ; 2,3,5-wheel
  (let ((wheel (vector 1 2 2 4 2 4 2 4 6 2 6)))
    (let loop ((n n) (f 2) (w 0) (fs (list)))
      (if (< n (* f f)) (reverse (cons n fs))
        (if (zero? (modulo n f))
            (loop (/ n f) f w (cons f fs))
            (loop n (+ f (vector-ref wheel w))
                  (if (= w 10) 3 (+ w 1)) fs))))))

We use a method called wheel factorization; at O(sqrt n), it’s not the fastest factoring algorithm available, but is sufficient for Project Euler. Next we use the factorization of n to determine the number of divisors of n by counting multiplicities of factors:

(define (numdiv n)
  (let ((fs (factors n)))
    (let loop ((prev (car fs)) (fs (cdr fs)) (f 2) (d 1))
      (cond ((null? fs) (* d f))
            ((= (car fs) prev) (loop prev (cdr fs) (+ f 1) d))
            (else (loop (car fs) (cdr fs) 2 (* d f)))))))

Now we can generate the triangle numbers, count the divisors of each, and stop when we’ve seen enough divisors:

(define (euler12 limit)
  (let loop ((n 1) (tri 1))
    (if (< limit (numdiv tri))
        (list n tri (numdiv tri))
        (loop (+ n 1) (+ tri n 1)))))

On my system, using Chez Scheme, that finds the answer in about a fifth of a second:

> (time (euler12 500))
(time (euler12 500))
    13 collections
    203 ms elapsed cpu time, including 0 ms collecting
    208 ms elapsed real time, including 0 ms collecting
    55714592 bytes allocated, including 54737696 bytes reclaimed
(12375 XXXXXXXX 576)

In the spirit of Project Euler, I elided the answer. You can run the program at https://ideone.com/LDfhXa.

Advertisements

Pages: 1 2

10 Responses to “Project Euler 12”

  1. Paul said

    I did many of the Euler problems some years ago.
    I posted the code on ideone.

  2. matthew said

    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.

  3. Globules said

    A Haskell version.

    import Data.List (find)
    import Math.NumberTheory.ArithmeticFunctions (divisorsList)
    
    tris :: [Integer]
    tris = 1 : zipWith (+) [2..] tris
    
    main :: IO ()
    main = print $ find ((>500) . length . divisorsList) tris
    
    $ ./euler12
    Just 76576500
    
  4. Zack said

    A Julia version…

    function NumberOfDivisors(n::Int64)
    m = n
    c = 2
    d = n / 2
    d_ = round(Int64, d)

    if d == d_
        c += 2
    end
    
    for x = 3:d_
        d = n / x
        d_ = round(Int64, d)
    
        if d == d_
            if x != d
                c += 2
            else
                c += 1
            end
    
            m = d_ - 1
        end
    
        if x >= m; break; end
    end
    
    return c
    

    end

    function main(n::Int64 = 1000000)
    x = 0
    nd = 0

    for i = 1:n
        x += i
        if x > 10000000; nd = NumberOfDivisors(x); end  # No point searching for a solution in such small numbers!
        if i%100 == 0; println(x, " ", nd); end                    # This is just to make sure the program doesn't stall, due to the computational burden
        if nd > 500; return x, nd; end
    end
    
    return x, nd
    

    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.

  5. Ernie said

    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

  6. matthew said

    @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).

  7. Ernie said

    Thank you Matthew for that observation. I only need to change one line of code (eliminate the merge) and the programme becomes 40% faster.

  8. Daniel said

    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;
          }
        }
      }
    }
    
  9. Steve said

    Mumps version: Converted Matthew’s python to Mumps. Thanks, @Matthew.

    Source code
    ndivisors ;
     ;
     n count,flg,last,n
     s (flg,last)=0,n=1
     f  d  q:flg
     . s count=$$ndivisors2($s(n#2=0:n\2,1:n))
     . i count*last>500 s flg=1 q
     . s n=n+1,last=count
     w !,n," ",n*(n-1)\2," ",count*last
     q
     ;
    ndivisors2(n) ;
     n i,sum
     s sum=0
     f i=1:1:n d
     . i n#i=0 s sum=sum+1
     q sum
    
    
    Run
    YDB>w !,$h d ^ndivisors w !,$h
    
    65144,30115
    12376 76576500 576
    65144,30126
    
  10. Steve said

    @ProgrammingPraxis: My apologies – I meant to blank out the middle answer, and forgot to do so prior to posting my answer.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: