Archive for June 19th, 2009

Problem 11 – Products in a Grid

June 19, 2009

Problem 11 seems somewhat straightforward, once you’ve got the grid read in from a file or so. I should probably feel pretty guilty about not doing these more generally, to accept products of a chosen number (instead of 4) of adjacent terms. Right now, though, I just feel tired.

Here’s what I came up with, as a first go:

def solve(grid):
    ret = 0
    for i in xrange(0,len(grid)):
        for j in xrange(0,len(grid[0]) - 4 + 1):
            p = grid[i][j] * grid[i][j+1] * grid[i][j+2] * grid[i][j+3]
            ret = max(ret, p)
            p = grid[j][i] * grid[j+1][i] * grid[j+2][i] * grid[j+3][i]
            ret = max(ret, p)
            p = grid[j][j] * grid[j+1][j+1] * grid[j+2][j+2] * grid[j+3][j+3]
            ret = max(ret, p)
    return ret

In the inner loop, I compute the products in each direction, starting at the index [i][j]. I realized you can improve this inner loop slightly, by not multiplying each product by grid[i][j] initially. Since all of them are being multiplied by that factor, you might as well find the biggest of the products without that factor, then multiply by it, and then take the max with the previous best result. And if grid[i][j] is 0, there’s no point in doing any products starting there. So you can replace the inner loop above with the following:

            if(grid[i][j]):
                p = max(grid[i][j+1] * grid[i][j+2] * grid[i][j+3],
                        grid[j+1][i] * grid[j+2][i] * grid[j+3][i],
                        grid[j+1][j+1] * grid[j+2][j+2] * grid[j+3][j+3])
                ret = max(ret, p*grid[i][j])

I then realized that this problem could be reduced to problem 8. In that problem, you were just looking for the biggest product among adjacent things in a list. The code I used, today, to solve that, is basically the code from the end of my post last time I talked about problem 8:

def bigprod(list, num):
    """ Largest product of num consecutive elements of list

    Assumptions: len(list) >= num, list contains only positive integers
    """
    p = reduce(lambda x,y:x*y, list[0:num-1], 1)
    f = 1 # a fake (for now) first digit of previous product
    ret = 0

    for idx in xrange(0, len(list) - num + 1):
        p = p // f
        p *= list[idx+num-1]
        f = list[idx]
        ret = max(ret, p)
    return ret

So it remains to see how to convert the grid to a list. We can get sort of “horizontally adjacent” terms by just working across rows. If we stick a 0 between rows, then products of adjacent terms in the list are only non-zero when they came from terms that were adjacent, horizontally, in the grid. We can do essentially the same thing for vertically adjacent, and diagonally adjacent, though the indexing takes slightly more thinking about for the diagonal bits. Here’s my code implementing this idea:

def solve(grid):
    """ The biggest product of 4 adjacent numbers in the grid

    Method: re-write grid as a bunch of lists, have bigprod to the work
    Assumption: the grid is square
    """
    adjH, adjV, adjMD, adjDU, adjDL = [], [], [], [], []
    for i in xrange(0, len(grid)):
        for j in xrange(0, len(grid)):
            adjH.append(grid[i][j])
            adjV.append(grid[j][i])
            if(i and i < len(grid) - 4 + 1 and j < len(grid) - i):
                adjDU.append(grid&#91;j&#93;&#91;j+i&#93;)
                adjDL.append(grid&#91;j+1&#93;&#91;j&#93;)
        adjH.append(0)
        adjV.append(0)
        adjDU.append(0)
        adjDL.append(0)
        adjMD.append(grid&#91;i&#93;&#91;i&#93;)

    # stick all of the adjacents together, they already end in 0 (besides MD)
    list = adjH + adjV + adjDU + adjDL + adjMD

    ret = 0

    while(len(list) >= 4):
        endidx = len(list)
        try:
            endidx = list.index(0)
        except:
            pass

        if(endidx >= 4):
            ret = max(ret, bigprod(list[0:endidx], 4))

        list = list[endidx+1:]

    return ret

I’ve not done a particularly formal performance comparison of these, but an informal look indicates that the first solution is much faster. It completes the problem on the assigned input (on my machine) in about .002 seconds, while the second solution takes .003 seconds. Big difference, no? I also fed them random grids of larger sizes. Looking at grids with size 200, the first solution takes along the order of .2 seconds, while the second is more like 1.5 seconds.

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.