Tag: python

Project Euler’s problem #11 statement goes:

In the 20×20 grid below, four numbers along a diagonal line have been marked in red.

08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08
49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00
81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65
52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91
22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80
24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50
32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70
67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21
24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72
21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95
78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92
16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57
86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58
19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40
04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66
88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69
04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36
20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16
20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54
01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48

The product of these numbers is [pmath]26 * 63 * 78 * 14 = 1788696[/pmath].

What is the greatest product of four adjacent numbers in any direction (up, down, left, right, or diagonally) in the 20×20 grid?

This one is remarkably easy but also was quite fun. I think it’s because it reminds me of the kind of work we’d do during our Algorithms classes during my first year in college. And so this one goes to my Algorithms teacher, Ricardo Vargas Dornelles—best teacher I’ve ever had too.

import operator
numbers = [[8,2,22,97,38,15,0,40,0,75,4,5,7,78,52,12,50,77,91,8],
[49,49,99,40,17,81,18,57,60,87,17,40,98,43,69,48,4,56,62,0],
[81,49,31,73,55,79,14,29,93,71,40,67,53,88,30,3,49,13,36,65],
[52,70,95,23,4,60,11,42,69,24,68,56,1,32,56,71,37,2,36,91],
[22,31,16,71,51,67,63,89,41,92,36,54,22,40,40,28,66,33,13,80],
[24,47,32,60,99,3,45,2,44,75,33,53,78,36,84,20,35,17,12,50],
[32,98,81,28,64,23,67,10,26,38,40,67,59,54,70,66,18,38,64,70],
[67,26,20,68,2,62,12,20,95,63,94,39,63,8,40,91,66,49,94,21],
[24,55,58,5,66,73,99,26,97,17,78,78,96,83,14,88,34,89,63,72],
[21,36,23,9,75,0,76,44,20,45,35,14,0,61,33,97,34,31,33,95],
[78,17,53,28,22,75,31,67,15,94,3,80,4,62,16,14,9,53,56,92],
[16,39,5,42,96,35,31,47,55,58,88,24,0,17,54,24,36,29,85,57],
[86,56,0,48,35,71,89,7,5,44,44,37,44,60,21,58,51,54,17,58],
[19,80,81,68,5,94,47,69,28,73,92,13,86,52,17,77,4,89,55,40],
[4,52,8,83,97,35,99,16,7,97,57,32,16,26,26,79,33,27,98,66],
[88,36,68,87,57,62,20,72,3,46,33,67,46,55,12,32,63,93,53,69],
[4,42,16,73,38,25,39,11,24,94,72,18,8,46,29,32,40,62,76,36],
[20,69,36,41,72,30,23,88,34,62,99,69,82,67,59,85,74,4,36,16],
[20,73,35,29,78,31,90,1,74,31,49,71,48,86,81,16,23,57,5,54],
[1,70,54,71,83,51,54,69,16,92,33,48,61,43,52,1,89,19,67,48]]

maxp = 0
for row in range(0,16):
        for col in range(0,16):
                maxp = max(maxp, reduce(lambda x,y: x*y, 
                                        [n for n in numbers[row][col:col+4]]))
                maxp = max(maxp, numbers[row][col] * numbers[row + 1][col + 1] * 
                                 numbers[row+2][col+2] * numbers[row+3][col+3])
                maxp = max(maxp, reduce(lambda x,y: x*y, 
                                        [n[col] for n in numbers[row:row+4]]))
                if col > 3:
                        maxp = max(maxp, numbers[row][col] * 
                                         numbers[row + 1][col - 1] * 
                                         numbers[row+2][col-2] * 
                                         numbers[row+3][col-3])

print maxp

As well, this runs in 0.020s, so not bad at all.

I decided to take on Project Euler’s problem #10. Its statement goes like this:

The sum of the primes below 10 is [pmath]2 + 3 + 5 + 7 = 17[/pmath].

Find the sum of all the primes below two million.

My first attempt used brute force and testing primality using only previously found primes:

primes = []

def is_prime(n):
        if not (n < 2 or any(n % x == 0 
           for x in filter(lambda x: x < math.ceil(n ** 2), primes))):
                primes.append(n)
                return True
        return False

sum_primes = 0
n = 1   
while n < 2000000:
    if (is_prime(n)):
        sum_primes += n
    n += 1  
print sum_primes

It worked if you consider an execution time of over 7 minutes reasonable. Clearly a bad solution. So I decided to rethink it a bit and tried again. This time, I decided to use a little creative thinking. For every number I tested for primality, I took the time to mark its multiples as not primes (when I first test, say, 3, I already know that [pmath]6, 9, 12, … n * 3 < 2000000[/pmath] will not be a prime numbers, so I don’t have to test their primality later.)

max_primes = 2000000
numbers = [0] * max_primes

def is_prime(n):
        if n <= 2:
                return True
        if numbers[n] == 1:
                return False
        c = 2
        m = 0
        while True:
                m = n * c
                if m < max_primes:
                        numbers[m] = 1
                else:
                        break
                c+= 1

        #if any(n % x == 0 for x in xrange(2, int(n ** 0.5) + 1)):
        i = 2
        while i < int(n** 0.5 + 1):
                if n % i == 0:
                        return False
                i += 1

        numbers[n] = 1
        return True

def sum_primes():
        soma = 0
        for i in range (2, max_primes):
                if is_prime(i):
                        soma += i
        return soma

print sum_primes()

This is a much better algorithm and it gave the correct result in 54 seconds. Project Euler has a “one-minute rule”, which state that all its problems should be solvable “on a modestly powered computer in less than one minute,” which means the solution applies. Still, it’s a brute force solution and I wanted a better one.

So today I picked my copy of Concrete Mathematics and read about primes. Honestly, I’ve been reading more about primes lately than I even thought I would.

Anyways, Knuth—yes, I know the book has three authors, but I can’t help thinking of it as a Knuth book—talks about several strategies to test primality, but it also mentions sieve algorithms that are used to generate lists of prime numbers. Exactly what I wanted!

def prime_list(limit):
        if limit < 2:
                return []
        sieve_size = limit / 2
        sieve = [True] * sieve_size

        for i in range(0, int(limit ** 0.5) / 2):
                if not sieve[i]:
                        continue
                for j in range((i * (i + 3) * 2) + 3, sieve_size, (i * 2) + 3):
                        sieve[j] = False

        primes = [2]
        primes.extend([(i * 2) + 3 for i in range(0, sieve_size) if sieve[i]])

        return primes

print reduce(lambda x,y: x+y, prime_list(2000000))

Et voilà! Correct result in 0.456s! The code is based on the description of the Sieve of Eratosthenes found in Concrete Mathematics. Also interesting to note, my idea in the second algorithm above is actually the basis of this sieve algorithm, I was just thinking “backwards.”

So the other night I was a bit bored and decided to do something to pass the time. I first came across Project Euler a while ago, but had never gone further than problem #1. Boredom is a great motivator and I went through problems #2 thru #9 last night and I decided to post my solutions in search of better ones. Feel free to comment with your suggestions.

Project Euler’s Problem #9 statement is —

A Pythagorean triplet is a set of three natural numbers, $latex a < b < c$, for which,

[latex]a² + b² = c²[/latex]

For example, $latex 3² + 4² = 9 + 16 = 25 = 5²$.

There exists exactly one Pythagorean triplet for which [pmath]a + b + c = 1000[/pmath].

Find the product [pmath]abc[/pmath].

Butt ugly, super slow, brute-force solution:

from sys import exit

def is_triplet(a, b, c):
        if a < b and b < c:
                return a ** 2 + b ** 2 == c ** 2

for a in range(0, 1000):
        for b in range(a + 1, 1000):
                for c in range(b + 1, 1000):
                        if a + b + c == 1000 and is_triplet(a, b, c):
                                print a * b * c
                                exit(0)

I’ll have to rethink this, as it’s really inefficient. As it is, it runs in 18.034s!

Updated: take a look at another take at this problem and algorithm.

Updated 2010/08/04: Gustavo Niemeyer pointed an obvious optimization to the algorithm above (written in C)—the innermost loop is unnecessary. I rewrote it in Python to see the difference:

from sys import exit

def is_triplet(a, b, c):
        return (a ** 2 + b ** 2 == c ** 2)

for a in range(1, 1000):
        for b in range(a + 1, (1000 - a) / 2):
                c = 1000 - a - b
                if is_triplet(a, b, c):
                        print "%d * %d * %d = %d" % (a, b, c, a * b * c)

From 18s to 0.076s. As well, following Eduardo Habkost’s suggestion, I used psyco and execution time went down to 0.023s.