Posts Tagged ‘prob12’

Problem 12 – Triangular Numbers

June 19, 2009

I went about counting divisors of a number t by looking at the prime factorization of t. If t=\prod p_i^{e_i}, then t has \prod 1+e_i divisors, because any divisor has a sort of “sub”-factorization, which is to say the same primes just to smaller (or possibly the same) powers (or 0, that’s why the “1+” is in each term in the product).

So, to determine the number of factors, we first find the prime factorization, and then do a quick product. Suppose we have a method “factor” which takes an integer n, and returns a list of tuples of the form (prime, power). So, for example, “factor(12)” would return “[ (2,2), (3,1) ]”. We can then pass such an object to the following function to calculate the number of divisors:

def numberDivisors(f):
    """ take a factorization and return the number of divisors """
    return reduce(lambda x,y:x*y, map(lambda t:t[1]+1,f), 1)

Now, I want to use this to find the number of divisors for triangular numbers. We could just loop over the triangular numbers, factoring each in turn. However, that’s pretty redundant. The n-th triangular number (let’s denote that by T_n) is \frac{1}{2}n(n+1). Exactly one of n or n+1 is even, so T_n=\frac{e}{2}\cdot o, where e is even and o is odd. If we factor \frac{e}{2} and o separately, we can multiply the factorizations together by simply adding powers. While we’re at it, if we store the factorization of \frac{e}{2} and o, and T_n, we’ll be able to use those factorizations to factor later triangular numbers. Probably if we thought about it for a while, we could get away with only storing the factorization of o and T_n, because the factorization of \frac{e}{2} won’t be used again (I think).

Here’s my implementation of these ideas:

def solve(goal):
    """ Find the first triangle number with more than 'goal' divisors

    Return the number and the time taken to find it
    Assumption: goal is a positive integer
    """
    t = time.time()

    # lets store some factorizations to kick things off
    factorization = { 2:[[2,1]], 3:[[3,1]], 4:[[2,2]], \
                      5:[[5,1]], 6:[[2,1],[3,1]] }

    # since goal is a positive integer, by assumption, we are looking
    # for a triangle number with >=2 divisors

    divs, n, tri, epsilon = 2, 2, 3, 0

    # loop invariant: T_n=tri has divs divisors, epsilon = 0 or 1
    # and n+epsilon is even
    while(divs <= goal):
        n, epsilon = n+1, 1-epsilon
        tri = n * (n+1) // 2
        lesser, greater = (n+epsilon)//2, n+1-epsilon
        if not (lesser in factorization):
            factorization&#91;lesser&#93; = factor(lesser)
        if not (greater in factorization):
            factorization&#91;greater&#93; = factor(greater)

        # concatenate the lists of factors
        combfact = factorization&#91;lesser&#93; + factorization&#91;greater&#93;

        # we know how to factor the product
        factproddict = {}
        for p,e in combfact:
            if not (p in factproddict):
                factproddict&#91;p&#93; = e
            else:
                factproddict&#91;p&#93; += e

        factorization&#91;tri&#93; = factproddict.items() # convert dictionary => list
        divs = divisors(factorization[tri])

    t = time.time() - t
    return (tri,t)

I guess things could possibly be simplified a little if the “factor” function that I use there returned a dictionary instead of a list, but currently I’ve got it set up to return a list. By the way, the factorization loop I’m using is the “factorAndSieve” code from Numerical Recipes, slightly modified to return a list of tuples, instead of just a (possibly repeating) list of primes. If I’d been paying slightly more attention, he’s got a lengthy post about computing the number of divisors for an integer as well.

I decided to compare the two methods of solving this problem, as I am one to do (the graphs are just so much fun). The first method just factors every triangle number in turn until you have enough divisors (called “Basic” in the following graph). The second is the one with the code above, where we store factorizations (called “Store”). In the following graph, the x-axis indicates the number of divisors that you are hoping for (tickmark labelled n corresponding to 10n divisors), and the y-axis is some scaled “ln(time)”.

Just for some sort of definite idea about timing, the last tickmark, corresponding to 500 divisors (the stated problem) takes my computer around 10 seconds for the Basic script, and less than 1 second for the Store script.

Of course, the Store script uses more storage space. I think we could remove factors as we went to save some space. At the n-th triangular number, I think we never need factorizations of anything less than n/2, so we could remove those as we went. Just as a side note, when you look for the first triangular number with more than 500 divisors, the Store script above has 21524 stored factorizations, with numbers ranging from 2 to the answer (which I got to be up in the 76 million range). If anybody can point me in the right direction to learn about how to test the space usage of my scripts, I’d love to hear about it. Something similar to the built in time module, but for space, would be awesome.

The function taking n to the number of divisors of T_n is probably a vaguely interesting function. I plotted it with inputs up to 100:

and with inputs up to 300:

That first big spike in the second graph is at input 224, indicating T_{224}=25200 has 90 divisors. Apparently.

Advertisements

Divisors

June 18, 2009

I haven’t posted in a while either.  So to compliment Chris’ post, I thought I’d put up the code I used to count the number of divisors in problem 12.  It isn’t very pretty, but it gets the job done

>>> import math
>>> def div(x):
          v=2
          y=2
          s=math.sqrt(x)
          if s==int(s):
                    v=3
          while y

Triangle numbers

June 17, 2009

I feel bad for not having worked on any problems since Problem 6. So here is my Mathematica solution to Problem 12.

Triangle[n_] := n (n + 1)/2;
n = 1;
While[Length[Divisors[Triangle[n]]] < 500, n++]; Print[Triangle[n]] [/sourcecode] As usual, I cheated by not programming in Python or Sage. But I figure this is better than nothing. I suspect I would get a lot more out of this programming project if I thought about how to write my own Divisors algorithm. Oh well.