Problem 19 – First Sundays

July 25, 2010 by

I don’t think I can state the question of problem 19 any more succinctly than the problem author:

How many Sundays fell on the first of the month during the twentieth century (1 Jan 1901 to 31 Dec 2000)?

I spent a few minutes perusing the datetime docs, and about an hour playing with a handful of variations on how to do this problem.

As usual, I like to start out with a stupid algorithm, just to get the right answer, so I know when other tries work. There’s not much simpler than starting on 1 Jan 1901 and adding days until you get to 31 Dec 2000, seeing when you’re on a day that’s a Sunday and the first of the month:

   count = 0
   curdate = date(1901, 1, 1)
   daydelta = timedelta(1) # 1 day
   while curdate.year < 2001:
      if curdate.weekday() == 6 and curdate.day == 1:
         count += 1
      curdate += daydelta

   print count

Of course, we’d do better to increment by weeks, so that we’re always on Sunday (as long as we start on a Sunday). Alternatively, we could increment by as many days as are in the month we’re in, but then we have to keep track of leap years. That never seems as clean, at least when I do it:

   count = 0
   daysinmonth = [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
   curdate = date(1901, 1, 1)
   while curdate.year < 2001:
      if curdate.weekday() == 6:
         count += 1
      curdate += timedelta(daysinmonth[curdate.month]-1)
      if curdate.year % 4 == 0 and (curdate.year % 100 != 0 or curdate.year % 400 == 0):
         curdate += timedelta(1)

   print count

Of course, I’m generally most amused if I can solve the problem in a single line. Allowing longish lines (perhaps ‘single statement’ would be more correct)…

print len(filter(lambda d:d.weekday()==6, [date(1901,1,1).replace(month=n,year=y) for n in range(1,13) for y in range(1901,2001)]))

Problem 40 – Digits of an Irrational

May 11, 2010 by

I scored myself a programming internship for the summer, so I thought it might be a good idea to actually do some programming. This will have to count…

In problem 40 we consider the irrational number defined by concatenating the positive integers, in order, and putting them to the right of the decimal place. We are asked for the product of some particular subset of digits.

As a solution that’s fairly easy to check for correctness, we can build up a string containing the irrational until we have enough digits, and then just pick out the digits we want. The stated problem only picks out powers of 10 as the determined digits, and so I threw together:

def solve(boundpow):
    irr = ""
    idx = 1
    bound = 10**boundpow
    while len(irr) < bound:
        irr += str(idx)
        idx += 1
    ret = reduce(lambda x,y:x*y,[int(irr[10**n - 1]) for n in xrange(0,boundpow+1)])
    return ret

I do like the use of reduce and list comprehension here. And I came up with the following that seems to be a quicker way to generate the irr string (plus it’s only one line):

irr = "".join(map(str, xrange(1,bound//(boundpow-1))))

This generates longer strings, but takes less time. It bothers me, a little, that I have to do “”.join instead of just sum. But what do I know. Anyway, it seems pretty pointless to keep track of the entire string. If larger powers of 10 were required, we’d start getting in trouble. Also, why require powers of 10?

So I have a slightly more general solution, which seems to also work (on the little testing I’ve done). Basically the idea is that I don’t need to remember all the digits that came before, when I go to concatenate the next integer. I just need to know the length of the string of digits before. And to be more general, my solve function takes in a list of integers (and assumes it is sorted in ascending order) indicating the position of the digits requested (position 1 being the first digit after the decimal).

def solve(digseq):
    lensofar = 0
    idx = 1
    ret = 1
    while len(digseq):
        nextstr = str(idx)
        nextlen = lensofar + len(nextstr)
        while len(digseq) and nextlen >= digseq[0]:
            ret *= int(nextstr[digseq[0] - lensofar - 1])
            digseq.pop(0)
        lensofar = nextlen
        idx += 1
    return ret

Seems to work out fine. Actually slower than the first solution, but whatever. It catches up if you go out one more power of 10 (though not if you use the faster generation of irr), and it’s more flexible.

Clearly there’s a more elegant solution, sorting out what integer the n-th digit would be part of basically algebraically. But I only care so much.

Problem 39 – Perimeters of Right Triangles

February 17, 2010 by

In problem 39 we’re supposed to look at right triangles with integral sides and perimeter no bigger than 100. Specifically, among such perimeters, which has the most representing right triangles?

The given bound seemed small enough that I figure brute force was a good place to start. Loop through lengths for the legs, calculate what the hypotenuse would be, keep track of solutions you find, organized by perimeter, and see who wins. The first (correct) program I wrote ran in 5 seconds. With 5 minutes tweaking the code without much thought, I had the runtime down to half a second. 1000 really isn’t a difficult bound (fwiw though, 10000 takes nearly a minute on my computer). Here’s my code:

def solve(bound):
    solns = [0 for n in range(0,bound+1)] # storage
    ret = 0 # the winning perimeter
    for smleg in xrange(1,bound//2 + 1):
        # the larger leg will still be shorter than the hypotenuse
        for lgleg in xrange(smleg+1, (bound-smleg)//2 + 1):
            hyp = int((smleg**2 + lgleg**2)**.5)
            peri = smleg + lgleg + hyp
            if peri <= bound and smleg**2 + lgleg**2 == hyp**2:
                solns[peri] += 1
                if solns[peri] > solns[ret]:
                    ret = peri
    return ret

I have a hard time letting such ugly code go, but I’m pretty sure I’m supposed to be spending my time doing other things anyway, so perhaps I should… move on to problem 40 🙂 Nearing a level-up on Project Euler…

Problem 38 – Pandigital Products

February 14, 2010 by

In problem 38 we’re supposed to find the biggest 9 digit pandigital number that is obtained by concatenating consecutive multiples of a fixed number, starting with the number itself. Sorta a mouthful, but the examples in the problem text make it pretty clear, I think.

10 minutes thought, or so, might convince you that (if the answer isn’t the example given in the problem text) the number you’ll be taking multiples of must be a four digit number, bigger than 9182. So you really don’t have too much to loop through, which is nice. Here’s what I came up with:

def solve():
    biglhs = 9182
    digits = map(str, xrange(1,10))
    for lhs in xrange(9183,10000):
        need = digits[:] # make a copy
        # get rid of digits in the lhs
        for d in str(lhs):
            if need.count(d): need.remove(d)
            else: break # quit early, MASSIVE speed improvement
        # if lhs is unique digits, and rhs is the remaining digits
        if len(need) == 5 and need == sorted(map(str, str(2*lhs))):
            biglhs = lhs
    biggest = "%s%s"%(biglhs,2*biglhs)
    return biggest

Quick problems are fun.

Problem 37 – Truncatable Primes

February 14, 2010 by

In problem 37 we are supposed to find the eleven primes whose decimal representations can be continuously truncated from the left, or continuously truncated from the right, producing primes the whole time. And the single digit primes don’t count. I’ll call these numbers truncatable primes. The truncation means the the right-most digit must always be a single digit prime, and the same for the left. Furthermore, 2s and 5s can only show up as the left-most digit. So, in fact, the right-most digit must be either 3 or 7.

So, with an oracle telling you about primes, you could just loop until you’ve found eleven truncatable primes. This ran in under 10 seconds on my computer, which is probably fine. But then I spent a while (much longer than I had hoped) writing a faster version. It’s better than real work.

Here’s what I’ll do. Keep track of a list of left-trunctable (can remove digits from the left) primes (I’ll call them LTs), beginning with the list [3, 7]. Loop through this list, looking for ways to add digits to the left and still have a LTs. When I find a new LT, add it to the list, and if it is also right-truncatable, keep it in a list of truncatable (both sides) primes.

Organizing the search for new LTs seems to be the interesting part of this setup. You can get more LTs, from a LT, by adding a single digit prime to the left, separated by a (possibly empty) string of 1s and 9s. However, you don’t want to go looking for all of the extensions to a given seed all at once (essentially a depth-first approach), because there might be infinitely many.

I decided to do my search for more LTs based on the length of their decimal representation. Begin by looking for 2 digit LTs, then 3 digits, and so on, keeping track of all the LTs you’ve found so far. When you want to look for LTs with the next number of digits, go through the list you’ve made so far and then only look for LTs that extend a given LT and also only have the correct number of digits. Since you’ve already got an LT, and will be adding a single digit prime to the very right, you’re essentially looking at how many 1s and 9s to throw in the middle.

In my code, I’ve got two functions “isltrunc” and “isrtrunc” that return True or False, if the input is a left-truncatable, or right-truncatable, prime (well, more correctly: a string containing the decimal representation of the prime). I’ve got another method that looks for LT extensions of a given LT, based on how may 1s and 9s to put in the middle. It looks like:

def ltruncs(p, midlen, oracle=PrimeOracle(100)):
    # first build all of the possible strings of 1s and 9s of length midlen
    intstrs = [""]
    while len(intstrs[0]) < midlen:
        intstrs = map(lambda s:"1"+s,intstrs) + map(lambda s:"9"+s,intstrs)

    growth = ["%s%s%s"%(d,i,p) for d in [2,3,5,7] for i in intstrs]
    return filter(lambda s:isltrunc(s,oracle), growth)

(By the way, my “isltrunc” and “isrtrunc” functions also take an oracle parameter)

Finally, my main code goes like this:

def solve():
    isprime = PrimeOracle(1000000)
    ret = set([]) # the truncatables
    seeds = ["3","7"] # the LTs
    strlen = 2
    while len(ret) < 11:
        # look for LTs of length strlen
        fringe = [] # new LTs we find
        for pstr in seeds:
            midlen = strlen - len(pstr) - 1
            if midlen >= 0:
                growth = ltruncs(pstr, midlen, isprime)
                fringe += growth

                # if we find new truncatables, add 'em to our store
                ret = ret.union([int(g) for g in growth if isrtrunc(g,isprime)])
                if len(ret) >= 11:
                    break
        strlen += 1
        seeds += fringe
    return sum(ret)

After building the PrimeOracle, this takes very little time to run.

I think I should be getting in the habit of using ’ to delineate strings, instead of “. And perhaps throwing appropriate unicode handling where it belongs. I have to get straight where that is though. Of course, none of this affects any of this work.

Problem 36 – Simultaneous Palindromes

February 10, 2010 by

In problem 36 we are asked to find all the numbers below a bound that are simultaneously palindromic in both base 10 and base 2. Not the most interesting problem, but for something like completeness sake (and probably to delay doing other things) I thought I’d post about it. I probably learned things, so it was probably worthwhile.

The problem text points out that the leading digit can’t be 0, in either base, so we know that all of the numbers to consider are odd. That’ll likely cut down list iteration by a factor of two. No matter how we set things up, we’ll likely want to reverse a string. Here’s what I rigged up:

def reverse(str):
    return "".join([str[len(str)-1-n] for n in xrange(0, len(str))])

Then I got curious if there was something built-in. There is a built-in way to reverse a list, but not a string. There’s also a “reversed” method that’ll reverse an iterator. And strings have iterators, so you could also do

def reverse(str):
    return "".join(reversed(str))

which seems to have the same affect as

def reverse(str):
    return "".join([i for i in reversed(str)])

It’d probably be fun to compare all of these for speed, but the strings we’re doing are only so big, so it seems to not matter a whole lot for this problem.

We’ll also likely need to test if a string is a palindrome. Here’s what I did:

def ispali(str):
    idx = 0
    while 2*idx <= len(str):
        if not str[idx] == str[len(str)-1-idx]:
            return False
        idx += 1
    return True

The idea being that you only need to go half-way in to see if a string is palindromic. I realized you could also just return “str == reverse(str)”. Mine looks to be slightly faster 🙂

Ok, so, the final bit is converting an integer to binary. It’s built-in (who would have guessed). “bin(n)” converts an integer to a string containing the binary representation of n. It also tacks a ’0b’ at the beginning, so we’ll drop that with “[2:]”.

I think all the bits are in place now. It’s easy enough to write the stupid loop that checks every single thing number up to the bound. In fact, you can do it in a line, if you’ve got an upper bound variable called “bound”:

sum([n for n in xrange(1,bound) if ispali(str(n)) and ispali(bin(n)[2:])])

Surely there’s more efficient code.

I decided to go about it by building up strings I already knew were palindromic decimal numbers. That is, make the right half, and flip it to get the left half. Of course, all of the strings will have even length, and so you’ll never get, for example, a 3 digit number. So I also toss in a pivot digit, which can either be empty or any single digit.

After a while I realized this had a slight glitch. Namely, if your right half is a valid decimal, it won’t have 0s as leading digits, so your palindromes will never have 0s in the middle (or more than one 1, if you set up 0 as a pivot digit), when clearly they should. So in addition to flipping the right side to make the left side, I also allow for middle pivots that are larger strings of 0s.

Here’s what I have:

def solve(digits):
    """ Sum of decimal and binary palindromes up to a bound

    The bound is given by a number of decimal digits.
    E.g., digits = 6 is a bound of 1000000

    Simulatenously palindromic, e.g. 585=1001001001_2
    """
    vals = [1, 3, 5, 7, 9] # all valid answers
    pivots = [''] + [str(n) for n in xrange(0,10)]
    pivots += ['0'*n for n in xrange(2,digits-1)]
    for rhs in xrange(1,10**(digits//2),2): # 2 because only odds
        rhsstr = str(rhs)
        lhs = reverse(rhsstr)
        for p in pivots:
            palstr = lhs+p+rhsstr
            if len(palstr) <= digits and ispali(bin(int(palstr))[2:]):
                vals.append(int(palstr))
    return sum(vals)

Good times. Also good times: writing this post in emacs, and using Eric Finster’s wplatex python package to post it. Well, once I got plastex working in my Ubuntu install. Something strange was going on with an import statement, but I hacked around it a little.

Problem 198 – Ambiguous Numbers

December 10, 2009 by

In problem 198, a number, x, is defined to be ambiguous for the denominator d if x has two best approximations by rationals with denominators no bigger than d. And then an x is ambiguous if it is ambiguous for some denominator. The problem is to determine how many ambiguous rationals there are in a certain interval, with a certain bound on the denominator.

Yes, I’m getting to this problem a little out of order. Jaime told me to. I was talking about Neighbors in Farey Sequences, on another blog, and he said this problem would fit well with what I was doing. In fact, if you don’t know about Farey sequences, you should probably go to that post (or Wikipedia) and read up, because I’m going to assume you know about them.

Once you have read a little about Farey sequences, you realize a number is ambiguous for a denominator if the number is the midpoint of neighbors in the appropriate Farey sequence. Moreover, the midpoint of Farey neighbors is always reduced, so you don’t have to worry about finding gcds. Also, Farey sequences can be easily generated inductively, by inserting the “mediant” between neighbors, when appropriate.

So a first go at this problem would be to use that idea. My first solutions are always ones that are completely stupid, but I’m pretty convinced are correct. So they’re slow, but good for testing optimizations or other solutions (for small inputs, anyway). Below is my third revision, approximately two orders of magnitude quicker than my first solution, an improvement gained by being slightly careful how far you iterate and what order you do things in. Oh, and I’m representing rationals as pairs of integers, to remain in the world of integer arithmetic. I worry about floating point operations when my denominators are allowed to be as big as 10^8.

def midpoint(lhs, rhs):
    """ Midpoint of lhs and rhs """
    return (lhs[0]*rhs[1] + lhs[1]*rhs[0], 2*lhs[1]*rhs[1])

def solve(r, D):
    """ Find ambiguous rationals

    Find ambiguous p/q with:
    1. p/q < 1/r (i.e. r*p < q)
    2. q <= D
    """
    seq = [(0,1), (1,1)] # F_1, the first Farey sequence
    den = 1
    ret = r < 2 and 1 or 0 # if 1/2 is included or not
    while den < D:
        # loop invariant: seq is the Farey sequence F_{den}
        # at least up to 1/r

        idx = 0 # our location in F_den
        # iterate to the end of F_den or until seq[idx] > 1/r
        while idx < len(seq)-1 and r*seq[idx][0] < seq[idx][1]:
            nextidx = idx + 1
            mediantden = seq[idx][1] + seq[idx+1][1]
            if(mediantden == den + 1): # new element of F_{den+1}
                mediant = (seq[idx][0]+seq[idx+1][0],mediantden)
                seq.insert(idx + 1, mediant)
                nextidx = idx + 2
                # possible new ambiguous numbers
                midleft = midpoint(seq[idx], mediant)
                midright = midpoint(mediant,seq[idx+2])
                if r*midleft[0] < midleft[1] and midleft[1] <= D:
                    ret += 1
                if r*midright[0] < midright[1] and midright[1] <= D:
                    ret += 1
            idx = nextidx
        den += 1
    return ret

This is still quite slow. I’m pretty sure it’d basically never solve the stated problem.

I spent some time thinking about ways to change the algorithm around and get the necessary speed improvements. I had spent some time thinking about neighbors in Farey sequences, and the extended Euclidean algorithm. I thought I could translate it into code that would iterate over reduced fractions, find neighbors for each, and count ambiguous numbers that way. It was a little messy, and then I realized I was over counting. And there seemed to be some boundary cases to worry about. I would still like to get a nice solution going along these lines.

But yesterday I came up with a recursive solution that runs quickly enough to solve the stated problem in about 10 minutes on my computer. That’s still 10 (or 100 :)) times longer than I usually shoot for on these problems, but I’ll take it (for now). It’s a pretty clean solution, I think. The idea is to loop through all the intervals [h/k, H/K] between neighbors in Farey sequences where (1) the interval contains points smaller than the specified bound, and (2) the product kK isn’t too big, for the specified denominator bound. Once you’ve got one such interval you can split the interval into two, at the mediant, and then run the same things on those two smaller intervals.

While the solution is essentially recursive, for whatever reason I ended coding it up as a loop. Here’s what I came up with:

def solve(r, D):
    """ assumes D >= 2 """

    # pairs will be my list of integrals to check
    # an element is (h,k), (H,K), neighbors in a Farey sequence
    # along with the denominator, 2*k*K, of the midpoint
    pairs = [ [(0,1),(1,1), 2] ]
    ret = 0
    while len(pairs):
        pair = pairs.pop()
        midden = pair[2] # we already know this is <= D
        midnum = pair[0][0]*pair[1][1] + pair[0][1]*pair[1][0]
        if r * midnum < midden: # midpoint is small enough, ambiguous
            ret += 1
        # construct the mediant
        mednum = pair[0][0] + pair[1][0]
        medden = pair[0][1] + pair[1][1]
        mediant = (mednum, medden)
        # and the denominators of the next midpoints
        leftden = 2 * pair[0][1] * medden
        rightden = 2 * pair[1][1] * medden
        # include each subinterval in the list of pairs
        # as long as the midpoint denominator isn't too big
        if leftden <= D:
            pairs.append( [ pair[0], mediant, leftden ] )
        if rightden <= D:
            pairs.append( [ mediant, pair[1], rightden ] )
    return ret

You can optimize this a bit. If you re-organize a little, don’t bother with tuples, and store some computed values, you can cut the running time nearly in half. That’s as much of an improvement as I’ve been able to make so far though. I look forward to doing better… but I do like the cleanliness of this solution. It’s 18 lines of code, and with some syntactic sugar (multiple assignments) could be reduced to or below 15, without harming readability much (perhaps depending on how readable you think it is now).

Problem 35 – Circular Primes

November 1, 2009 by

In problem 35 we are asked to find the number of “circular primes” below a bound. A prime is circular if all of the numbers obtained by cyclic permutation of the original number are also prime.

In my first solution, I’ll use my PrimeOracle from a while ago:

from PrimeOracle import PrimeOracle

def iscircularprime(num, isprime):
    if not isprime[num]:
        return False
    numstr = str(num)
    cyclics = [ int(numstr[i:] + numstr[:i]) for i in range(0,len(numstr)) ]
    return reduce(lambda a,b:a and b, map(lambda n: isprime[n], cyclics))

def solve(bound):
    """ Find how many circular primes there are below bound """
    isprime = PrimeOracle(bound)
    return len(set([n for n in range(2,bound) if iscircularprime(n, isprime)]))

This is pretty straightforward. I like how we can make all of the cyclic permutations using that one line (line 7), using the slice notation. And testing if all elements in a list are “True”, using reduce, is nice. I had hoped that maybe instead of the lambda expression I wrote, I could just pass in “and” as the first argument to that reduce. No dice though. Anybody know a nice way? I’ve also wondered about passing in the identity function in that first spot (in a filter, not a reduce), but don’t know how (without a lambda function or so).

Anyway, I think it’s good to have a solution like this that is easy to write and easy to read (hopefully) and, therefore, easy to check for correctness. Then I like to try to make it faster. The code above checks each cyclic permutation of each number many times, which is pointless.

Here’s a slightly faster solution, which doesn’t invoke the PrimeOracle:

def solve(bound):
    """ Find how many circular primes there are below bound """
    # isprime[n] if 2n+3 is prime, eventually
    # so to test if m is prime, check isprime[(m-3)/2]
    isprime = [True for n in range(0, bound//2)]
    for val in xrange(3, bound, 2):
        if isprime[(val-3)//2]:
            for mul in xrange(val*val, bound, 2*val):
                isprime[(mul-3)//2] = False

    ret = 1 # 2 isn't included in isprime, but it is circular
            # we've assumed bound > 2

    # let isprime[n] mean 2n+3 is a prime we haven't checked for circular primality

    And = lambda a,b: a and b # since we'll use it twice

    for val in xrange(3, bound, 2):
        if isprime[(val-3)//2]:
            valstr = str(val)
            circs = [ int(valstr[i:] + valstr[:i])
                      for i in range(0,len(valstr))]
            if reduce(And, map(lambda d:d%2, map(int, valstr))): # no even digits
                if reduce(And, map(lambda n:isprime[(n-3)//2], circs)): # all circs are prime
                    ret += len(set(circs))
            for n in circs:
                isprime[(n-3)//2] = False

    return ret

In line 25, where we increment ret, my original thought was to add len(valstr), since that’s how many cyclic permutations there are. However, these cyclic permutations need not be distinct, so I switched to making a set and getting its length.

What do you think of lines 23 and 24, an if within an if, and not much else? Would it be better (more readable?) to have simply “if () and ()”, even though that’d look best on two lines anyway? I guess python wouldn’t have to make one more nested context for that inner if, if we used 1 instead of 2…

It has, after writing all of this, occured to me that there’s some issues if “bound” isn’t a power of 10. Should we consider only primes such that all cyclic permutations are also primes less than the bound? We’d also have to think about how big to make our isprime array. So it would probably be better to re-write the above code accepting the bound as a number of digits.

Problem 34 – Factorials of Digits

October 31, 2009 by

In problem 34 we are supposed to find the sum of all numbers that are the sum of the factorials of their digits.

My first thought was to just loop through all the possible numbers, compute sums of factorials of digits, and see what numbers worked. To do this, we have to know a stopping point. Well, 9! is 362880, so the largest digit factorial sum (DFS) that could be made with n digits is 362880n. An n digit number is always at least as big as 10^{n-1}, so if 362880n<10^{n-1}, then no number with n digits (or more) will be equal to its DFS. It turns out that 8 digits is too many. So we could stop looking once we got to 7\cdot 9!=2540160. We could clearly lower this bound some more, but let’s go with it for now.

Suppose we’ve got a function “fact” that returns the factorial of whatever number you give it. Then we can compute DFS using the following function:

def dfs(n):
    """ Compute the sum of the factorials of the digits of n """
    return sum(map(fact, map(int, str(n))))

I think that’s fun. Convert n to a string, then change the string to a list of digits (as integers), apply “fact” to each, and then sum. I’m sure I’ve used similar lines in other problems here.

Now, given this, we could, in theory, solve this problem with one more line:

print sum([n for n in xrange(10,7*fact(9)) if dfs(n) == n])

This certainly has the advantage of being readable. But I think it could be improved a bit, as it’s fairly slow.

The problem with this code is that we can quickly ignore lots of numbers, without computing the digit sum for each. If a number starts with a bunch of 1s, the DFS will only be so big. On the flip side, there’s no point in considering any 5 digit numbers with a 9 in them. Another issue with the code above, I expect, is that the dfs calculation used is not particularly fast.

Here’s some more efficient code:

fact = [1,1,2,6,24,120,720,5040,40320,362880]
digs = range(0,10)

def solvedigs(k, head = [], headval = 0, sumhead = 0):
    """ All solutions with k digits, starting with head

    assumes:
      len(head) <= k
    """
    if len(head) == k: # we have a k digit number
        if headval == sumhead:
            return headval
        else:
            return 0
    rem = k-len(head) # number of digits remaining to append
    smallval = headval*(10**rem) # smallest value that could be made
    smallsum = sumhead + rem*1 # smallest digit sum possible
    largeval = smallval + 10**rem - 1 # largest value that could be made
    largesum = sumhead + rem*fact[9] # largest digit sum possible
    if(largeval < smallsum or smallval > largesum):
        return 0
    return sum([solvedigs(k, head + [d], 10*headval + d, sumhead + fact[d])
                for d in digs
                if len(head) or d>0]) # no 0 leading digits!

if __name__ == "__main__":
    print sum([solvedigs(k) for k in xrange(2,8)])

I know I cheated with the global variable fact, but it is quicker than computing factorials every time we want to know them. And with problems of this size, having global variables really shouldn’t be an issue.

In the code above I pass around the digits of the number I’m considering, the value that digit string represents, and the sum of the factorials of those digits. Never involving strings (more precisely, string to/from integer conversions) seems like a good idea for efficiency in this problem.

Problem 33 – Cancelling Fractions

October 30, 2009 by

49/98 = 4/8 if you cancel the 9s. This happens for few others fractions n/d where 0<n<d<100, and we’ll also ignore (n'\cdot 10)/(d'\cdot 10) for digits n',d'. Now take all these fractions that cancel like this, consider their product as a reduced fraction, and find the denominator. This is the challenge of problem 33.

Here’s what I came up with:

def gcd(a,b):
    if b == 0:
        return a
    return gcd(b,a%b)

def solve():
    fracs = [ ]
    conv = lambda t:int("%s%s" % t)
    for bothdig in xrange(1,10): # the digit that cancels
        for numdig in xrange(1,10): # digit left in numerator
            for dendig in xrange(1,10): # digit left in denom
            nums = map(conv, [(numdig,bothdig),(bothdig,numdig)]) # consider bother orderings
            dens = map(conv, [(dendig,bothdig),(bothdig,dendig)])
            for v in [(n,d) for n in nums for d in dens if n<d]:
                if v[0]*dendig == v[1]*numdig: # v[0]/v[1] = numdig/dendig
                    fracs += [ v ]
    p = reduce(lambda l,r:(l[0]*r[0],l[1]*r[1]),fracs) # multiply together solutions
    return p[1] / gcd(*p) # reduce the fraction, return denom

Certainly room for improvement, but it works.