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...

2009-02-04

Project Euler: Problem 5

At last! A Project Euler problem I didn't brute-force. This problem was a fairly straight-forward LCM problem. I just reused my factorize function from before and implemented the algorithm to find LCM by prime factorization.


"""Solves Problem 5 of Project Euler."""

def factorize(to_factor):
    """Use trial division to factorize to_factor and return all the resulting \
    factors."""
    factors = []
    divisor = 2
    while (divisor < to_factor):
        if not to_factor % divisor:
            to_factor /= divisor
            factors.append(divisor)
            # Note we don't bump the divisor here; if we did, we'd have
            # non-prime factors.
        elif divisor == 2:
            divisor += 1
        else:
            # Trivial optimization: skip even numbers that aren't 2.
            divisor += 2
    if not to_factor % divisor:
        # Don't forget the last factor
        factors.append(to_factor)
    return factors

def lcm(numbers):
    """Finds the Least Common Multiple of numbers."""
    highest_degree_factors = {}
    for number in numbers:
        degrees_by_factor = {}
        for factor in factorize(number):
            # Translate the raw list of factors into a dictionary of degrees
            # keyed on the factor.
            current_degree = degrees_by_factor.setdefault(factor, 0)
            degrees_by_factor[factor] = 1 + current_degree

        # Update the top-level dict so it really is tracking the highest
        # degrees.
        for k, v in degrees_by_factor.iteritems():
            highest_degree_factors.setdefault(k, v)
            if highest_degree_factors[k] < v:
               highest_degree_factors[k] = v 

    running_product = 1
    for factor, degree in highest_degree_factors.iteritems():
        running_product *= factor ** degree

    return running_product

if __name__ == '__main__':
    print lcm(range(1,20))

Back to flipping out...

Project Euler: Problem 4

I really enjoyed solving this one; palindromes have intrigued me since learned about them in grade school. Rather than repeat my mistake of getting too cute from Problem 3, I decided to brute-force this one. That approach did lead me to learn some nifty Python tricks though:

  • [::-1] for reversing a string
  • the new reversed feature
  • xrange for lazy range generation; I'd seen this one, but this is the first time I'd written it

Without further ado, here's my solution to Problem 4.


"""Solves Problem 4 from Project Euler.
"""

def is_palindrome(to_test):
    """Determine whether to_test is a palindrome."""
    # Uses extended slice syntax to reverse a sequence.
    return to_test == to_test[::-1]

def problem_4():
    """Function to solve problem 4."""
    # I arrived at the values in the range by finding the smallest and largest
    # products of 2 3-digit numbers.
    for n in reversed(xrange(10000, 998001)):
        if is_palindrome(str(n)):
            for factor_1 in reversed(xrange(1000)):
                for factor_2 in reversed(xrange(1000)):
                    if n == (factor_2 * factor_1):
                        return n

    return 255  # Return a sentinel value to indicate failure.

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

Back to flipping out...

2009-02-02

Project Euler: Problem 3

I just finished playing around with Problem 3 from Project Euler. My rough idea for the problem was to build a list of the primes up to the number in question and start testing them, since I already knew the basics of the Sieve of Eratosthenes and I didn't know any good algorithms for factoring. After I finished, I looked up some of those algorithms (see Wikipedia's (incomplete) list), and I probably wouldn't have bothered with anything that complex anyway. Here's my original "solution" (I use scare quotes since I know it has bugs, and I knew it at the time).


"""Solves Problem 3 from Project Euler."""

import math

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 find_highest_prime_factor(to_factor):
    """Find the highest prime factor of to_factor."""
    highest_prime = None
    primes = e_sieve(int(math.sqrt(to_factor)))
    for prime in primes:
        if not to_factor % prime:
            quotient = to_factor / prime
            if quotient in primes:
                return quotient
            highest_prime = prime

    if highest_prime:
        return highest_prime

    return to_factor

if __name__ == "__main__":
    print find_highest_prime_factor(600851475143)

After I coded up my solution and verified I had the right answer (after what seemed like forever but was really ~15 minutes), I just knew there had to be a faster way. A little work with cProfile quickly revealed I was spending all my time in the Sieve. It turns out that even a "naïve" algorithm for factoring is better than generating the primes for a number this size. My new, "naïve" solution ran in 44 ms.


"""Solves Problem 3 from Project Euler."""

def factorize(to_factor):
    """Use trial division to factorize to_factor and return all the resulting \
    factors."""
    factors = []
    divisor = 2
    while (divisor < to_factor):
        if not to_factor % divisor:
            to_factor /= divisor
            factors.append(divisor)
            # Note we don't bump the divisor here; if we did, we'd have
            # non-prime factors.
        elif divisor == 2:
            divisor += 1
        else:
            # Trivial optimization: skip even numbers that aren't 2.
            divisor += 2
    if not to_factor % divisor:
        # Don't forget the last factor
        factors.append(to_factor)
    return factors

if __name__ == "__main__":
    print max(factorize(600851475143))

As an added bonus, this code has fewer bugs. One of these days, I'll learn to follow KISS. One of these days….

Back to flipping out...

2009-01-31

Project Euler: Problem 2

Here's my solution to Problem 2.

Back to flipping out...

Project Euler: Problem 1

As part of my ongoing efforts to learn Python, I've decided to start working through the Project Euler problem sets. If you don't want to see spoilers, then leave now. Here's the Python script I whipped out for Problem One.

Back to flipping out...

2009-01-30

Sneak Attack—Send IMs from Quicksilver

If you're using Adium and Quicksilver, and you hate having to mouse over to the contacts list and double-click, this one's for you.

Back to flipping out...