2009-02-11

Project Euler: Problem 10

I just finished up Problem 10; given my earlier work on Problem 7, it was trivial to adapt it and arrive at the following program.


"""Solves Problem 10 for Project Euler."""
import math

def is_prime(candidate, known_primes):
    """Determines whether candidate is prime by trial division using \
    known_primes.

    For this function to work, known_primes *must* be accurate.
    """
    last_possible = math.sqrt(candidate)
    for current_prime in known_primes:
        if current_prime > last_possible:
            break
        if not candidate % current_prime:
            return False
    return True

def primes_generator(upper_bound):
    """A generator for all the primes < upper_bound."""
    candidates = xrange(2, upper_bound)
    primes = []
    for n in candidates:
        if is_prime(n, primes):
            primes.append(n)
            yield n

def problem_10():
    """Sum the primes less than 2 million."""
    return sum(primes_generator(2000000))

if __name__ == '__main__':
    print problem_10()

Back to flipping out...

2009-02-10

Project Euler: Problem 9

Problem 9 deals with one of the more interesting things I learned in high school geometry: Pythagorean triples. In high school, I just memorized the 2 most common ones (3, 4, 5 and 5, 12, 13) and thought to myself: Wouldn't it be cool if I could generate all of these? But that was in the before time; Wikipedia didn't exist and my text book wasn't cool enough to dwell on them. At any rate, now I know Euclid's formula for generating Pythagorean triples: m ** 2 - n ** 2, 2 * m * n, m ** 2 + n ** 2

Here's the script I used to solve the actual problem:


"""Solves Problem 9 from Project Euler."""

import operator

def triples(upper_bound):
    """Generator for Pythagorean Triples (represented as tuples).

    Uses Euclid's formula to generate Pythagorean Triples
    (see http://en.wikipedia.org/wiki/Pythagorean_triple#Generating_a_triple).
    """
    for m in xrange(2, upper_bound):
        for n in xrange(1, m):
            yield m ** 2 - n ** 2, 2 * m * n, m ** 2 + n ** 2
# Only uncomment if the triples we get from the original Euclid's are insufficient.
# for k in xrange(1, upper_bound):
#     yield k * (m ** 2 - n ** 2), k * (2 * m * n), k * (m ** 2 + n ** 2)

def problem_9():
    """Finds the product of the Pythagorean Triple where a + b + c = 1000."""
    for triple in triples(1000):
        if sum(triple) == 1000:
            return reduce(operator.mul, triple)
    return 0

if __name__ == '__main__':
    print problem_9()

Back to flipping out...

2009-02-09

Project Euler: Problem 8

Just finished up with Problem 8. Brute-forcing it was pretty straightforward, so I decided to play about with some of the more functional aspects of Python. Enter reduce. Here's the original version I used to solve the problem.


"""Solves Problem 8 from Project Euler."""

def problem_8(num_in_question):
    """Finds and returns the greatest product of 5 consecutive digits \
    of num_in_question."""
    to_process = str(num_in_question)
    offset = 0
    highest_product = 0
    last_possible_start = len(to_process) - 5
    while (offset < last_possible_start):
        digits = [int(digit) for digit in to_process[offset:offset + 5]]
        product = 1
        for n in digits:
            product *= n
 
        if product > highest_product:
            highest_product = product
 
        offset += 1
 
    return highest_product

if __name__ == '__main__':
    print problem_8("73167176531330624919225119674426574742355349194934\
96983520312774506326239578318016984801869478851843\
85861560789112949495459501737958331952853208805511\
12540698747158523863050715693290963295227443043557\
66896648950445244523161731856403098711121722383113\
62229893423380308135336276614282806444486645238749\
30358907296290491560440772390713810515859307960866\
70172427121883998797908792274921901699720888093776\
65727333001053367881220235421809751254540594752243\
52584907711670556013604839586446706324415722155397\
53697817977846174064955149290862569321978468622482\
83972241375657056057490261407972968652414535100474\
82166370484403199890008895243450658541227588666881\
16427171479924442928230863465674813919123162824586\
17866458359124566529476545682848912883142607690042\
24219022671055626321111109370544217506941658960408\
07198403850962455444362981230987879927244284909188\
84580156166097919133875499200524063689912560717606\
05886116467109405077541002256983155200055935729725\
71636269561882670428252483600823257530420752963450")

And here's the same problem_8 function using reduce:


def problem_8(num_in_question):
    """Finds and returns the greatest product of 5 consecutive digits \
    of num_in_question."""
    to_process = str(num_in_question)
    offset = 0
    highest_product = 0
    last_possible_start = len(to_process) - 5
    while (offset < last_possible_start):
        digits = [int(digit) for digit in to_process[offset:offset + 5]]
        product = reduce(operator.mul, digits)
 
        if product > highest_product:
            highest_product = product
 
        offset += 1
 
    return highest_product

Pretty similar: 1 less line of code, 1 more line of imports, almost identical performance. I guess it all comes down to taste. One note: if you're doing functional programming and need to use a function supported by the operator module, that's the recommended way of doing it. Since it's part of the standard library, it's more obvious what's going on than a comparable lambda, plus they're implemented in C to give better performance. But since we're getting all functional, we might as well do it all the way. Here's another version:


def problem_8(num_in_question):
    """Finds and returns the greatest product of 5 consecutive digits \
    of num_in_question.

    This function expects num_in_question to be a string so we can
    slice it into 5-digit sequences.
    """
    SEQUENCE_LENGTH = 5
    sequences = [num_in_question[offset:offset + SEQUENCE_LENGTH] \
        for offset in range(len(num_in_question) - SEQUENCE_LENGTH)]
    nums = []
    for sequence in sequences:
        nums.append([int(num) for num in sequence])

    return max([reduce(operator.mul, num_list) for num_list in nums])

Back to flipping out...

Sneak Attack: Graph of NP-Complete Problems

Mmm… NP-completeness….

Back to flipping out...

2009-02-06

Project Euler: Problem 7, redux

Remember how I mentioned the Sieve wasn't the most performant solution? Here's a much faster solution. On my system the time dropped from 11.661 seconds to .332 seconds. Also, it makes use of one of my favorite features in Python, so far: generators.


"""Solves Problem 7 from Project Euler."""

import math
import sys

def is_prime(candidate, known_primes):
    """Determines whether candidate is prime by trial division using \
    known_primes.

    For this function to work, known_primes *must* be accurate.
    """
    last_possible = math.sqrt(candidate)
    for current_prime in known_primes:
        if current_prime > last_possible:
            break
        if not candidate % current_prime:
            return False
    return True

def primes_generator():
    """A generator for all the primes <= sys.maxint."""
    candidates = xrange(2, sys.maxint)
    primes = []
    for n in candidates:
        if is_prime(n, primes):
            primes.append(n)
            yield n

def problem_7(n):
    """Finds the nth prime number."""
    primes = primes_generator()
    i = 0
    while i < n - 1:
        i += 1
        primes.next()

    return primes.next()

if __name__ == '__main__':
    print problem_7(10001)

Back to flipping out...

Project Euler: Problem 7

My new math trick of the day? Upper and lower bounds on the nth prime:

n * ln(n) + n * ln(ln(n - 1)) is less than p[n] is less than n * ln(n) + n * ln(ln(n)) for n greater than or equal to 6

Using this little nugget, I can give myself an upper bound on the numbers I need to test for primality and unleash the Sieve of Eratosthenes on Problem 7. BTW, I realize the Sieve isn't the most performant way to do most of these tests; I just think it's a really elegant solution for finding primes and it isn't slow enough to be an issue for most numbers of this size. At any rate, here's my script:


"""Solves Problem 7 from Project Euler."""

import math

FIRST_SIX_PRIMES = {1:2, 2:3, 3:5, 4:7, 5:11, 6:13}

def e_sieve(upper_bound):
    """Uses the Sieve of Eratosthenes to get a list of the primes up to max."""
    primes = []
    candidates = range(2, upper_bound)
    while candidates:
        head = candidates[0]
        primes.append(head)
        candidates = [n for n in candidates[1:] if n % head]

    return primes

def problem_7(n):
    """Finds the nth prime number.

    Thanks to
    http://en.wikipedia.org/wiki/Prime_number_theorem#Approximations_for_the_nth_prime_number
    we know that it will be less than n * ln(n) + n * ln(ln(n))
    (for n >= 6).
    """
    if n < 6:
        return FIRST_SIX_PRIMES[n]

    upper_bound = int(n * math.log(n) + n * math.log(math.log(n)))

    primes = e_sieve(upper_bound)
    return primes[n - 1]

if __name__ == '__main__':
    print problem_7(10001)

Back to flipping out...

2009-02-05

Project Euler: Problem 6

I just finished up Problem 6. This one was really straightforward. The interesting part was learning a new math trick: Square Pyramidal Numbers. Here's the most brute-force way (I thought of) to solve it:

"Naïve" Solution


"""Solves Problem 6 from Project Euler."""

def sum_of_squares(upper_bound):
    """Sums the squares of all the natural numbers from 1 to \
    upper_bound (inclusive)."""
    return sum(n ** 2 for n in range(1, upper_bound + 1))

def square_of_sums(upper_bound):
    """Sums all the numbers from 1 to upper_bound."""
    return (sum(range(1, upper_bound + 1)) ** 2)

def problem_6(n):
    """Finds the difference between the square of the sums and the \
    sum of the squares for the first n natural numbers."""
    return square_of_sums(n) - sum_of_squares(n)

if __name__ == '__main__':
    print problem_6(100)

Even though this solution has the worst Big O, it wasn't noticeably slower than any of the others for numbers the size of the ones here.

I used a slightly more sophisticated approach because I remembered an insight for summing the first N natural numbers, discovered by none other than Gauss himself (especially entertaining to use Gauss to solve problems on Project Euler, I know): n * (n + 1)) / 2. I fell back on brute force for the sum of squares part, though.

Solution I Used


"""Solves Problem 6 from Project Euler."""

def sum_of_squares(upper_bound):
    """Sums the squares of all the natural numbers from 1 to \
    upper_bound (inclusive)."""
    return sum(n ** 2 for n in range(1, upper_bound + 1))

def square_of_sums(upper_bound):
    """Sums all the numbers from 1 to upper_bound."""
    # Use Gauss' insight about n/2 * (n + 1).
    return (upper_bound * (upper_bound + 1) / 2) ** 2

def problem_6(n):
    """Finds the difference between the square of the sums and the \
    sum of the squares for the first n natural numbers."""
    return square_of_sums(n) - sum_of_squares(n)

if __name__ == '__main__':
    print problem_6(100)

Since there was a math trick for solving the summing problem, I suspected there was one for the sum of squares problem, as well. A quick Google revealed: Square Pyramidal Numbers. I wish I had known about this when I was doing stupid brain teasers; I always knew I was wasting too much time counting squares.

Least Brutish Solution


"""Solves Problem 6 from Project Euler."""

def sum_of_squares(upper_bound):
    """Sums the squares of all the natural numbers from 1 to \
    upper_bound (inclusive)."""
    # Use the rule for square pyramidal numbers
    # (http://en.wikipedia.org/wiki/Square_pyramidal_number).
    return ((2 * upper_bound ** 3) + (3 * upper_bound ** 2) + upper_bound) / 6

def square_of_sums(upper_bound):
    """Sums all the numbers from 1 to upper_bound."""
    # Use Gauss' insight about n/2 * (n + 1).
    return (upper_bound * (upper_bound + 1) / 2) ** 2

def problem_6(n):
    """Finds the difference between the square of the sums and the \
    sum of the squares for the first n natural numbers."""
    return square_of_sums(n) - sum_of_squares(n)

if __name__ == '__main__':
    print problem_6(100)

All in all, I'm really glad I started working these problems. I've usually learned/remembered at least one interesting math trick or Python trick during each one I solved, or at least had good programming practice (KISS) driven home.

Back to flipping out...