Archive for the 'Programming' Category

Firefox 3 is a big step forward

March 13, 2008

Firefox has released Beta 4 of Firefox 3 this week. I have periodically checked the new releases as they’ve been posted, but it was only this version that impressed me enough to make it my default browser. The feel of the browser is much faster than previous versions. As some of the developers have posted, the memory use is much improved. I like the new bookmark management system Places, though I can’t really claim to “get” it yet. I can’t really point to any of these as the main reason for this version becoming my default when the others did not, and I will sorely miss Google Browser Sync, until they get that extension to work with v3.0. But the overall user experience with 3.0B4 is wonderful.


Scientific packages in Python I just discovered

February 20, 2008

I like to think that I’m abreast of the best of the scientific packages in Python, but I just discovered several cool packages today:

  • SymPy, a symbolic math package in Python, which actually is part of the Sage program that I like so much;
  •  sfepy, a finite element analysis package in Python (also symfe, for symbolic finite element analysis).

Testing whether a number is a square integer

February 13, 2008

In an earlier post I mention that I didn’t know of a good method to determine whether a number was an integer squared or not. I spent some time looking into this and here’s what I came up with.

The naive solution

squares = Set(i*i for i in xrange(1,BIG)); def is_square(n): return n in squares

doesn’t really work for large values of BIG, since it has to compute all of those elements.

Googling around led me to the beast that is the integer square root. The wikipedia page led me to the following page, which had the following Python implementation:

def isqrt(n):
    xn = 1
    xn1 = (xn + n/xn)/2
    while abs(xn1 - xn) > 1:
        xn = xn1
        xn1 = (xn + n/xn)/2
    while xn1*xn1 > n:
        xn1 -= 1
    return xn1

Thus one can simply use:

def is_square(n): return isqrt(n)**2 == n


D’oh moment from Project Euler

February 12, 2008

I don’t think that the caliber of my posts has necessarily been so high that people will be shocked by overt stupidity from me, but I just discovered some in myself, and stupidity is one of the things that is better shared.

One of the things I’ve had to do very often in Project Euler is figure out whether some long number is an integer or not. So I’ve defined the following function, that I’ve used over and over again:

def is_int(f): return int(f) == f

It actually works pretty well in any number of cases, until the floating point number f starts getting big. However, that doesn’t mean it should be used all of the time.

One place where I’ve used this function incorrectly is figuring out whether a square root is integral, for example, in looking for Pythagorean triplets. It’s much better to test whether a**2+b**2 is square, than whether sqrt(a**2+b**2) is integral. However, since for really huge numbers this requires generating a long list of square integers, I don’t feel too stupid about this. (Does anyone know of a clever way to test whether an integer is a square??)

What I do feel stupid about is in checking whether a rational fraction is integral. I caught myself doing this:


which is staggeringly stupid. I was flipping through some of the forum solutions to Project Euler, and noticed one of the Pythonistas was using the much more intelligent

a % b == 0

One of those things that seems stunningly obvious once you see it, but, for whatever reason, I didn’t see it. (Hangs head in shame.)

Thoughts on Pascal’s triangle in Python

February 11, 2008

For a Project Euler problem I had to work on a generator for Pascal’s Triangle. Ultimately, it didn’t emerge as the right way to solve the problem, but I arrived at an interesting implementation that I’d like to share.

The naive implementation is:

def gen_pascal():
    last = []
    while True:
        n = len(last)
        row = [1]*(n+1)
        for i in range(1,n):
            row[i] = last[i-1]+last[i]
        last = row
        yield row

This is pretty similar to any number of implementations you can find on the web. Works quite nicely, too, except it requires 2*n storage, which gets big if you’re interested in really large triangles.So I started thinking about automatic ways of generating only a row at a time, and yielding an iterator for that row rather than a list, to reduce the memory usage. Again, the obvious way to generate the elements are:

def fact(n):
    if n <= 0: return 1
    return prod(xrange(1,n+1))def combos(n,m): return fact(n)/fact(n-m)/fact(m)

def combos_slower(n,m):
    if m == 0 or n == m: return 1
    return combos(n-1,m-1)+combos(n-1,m)

def gen_pascal_row(n):
    for r in range(0,n+1):
        yield combos(n,r)

So this version doesn’t use any memory, but actually takes a great deal of computation (roughly N^3??) per new element, which seems a bit wasteful, and starts to bog down even sooner than the naive version.

I played around with an alternate way to generate the combinations, in part inspired by the naive version:

 def combos_slower(n,m):
    if m == 0 or n == m: return 1
    return combos(n-1,m-1)+combos(n-1,m)

but, as the name implies, it was even slower still.

Ultimately I settled on this version:

def gen_pascal_row_fast(n):
    el = 1
    yield el
    for m in range(1,n+1):
        el = (el*(n-m+1))/m
        yield el

which only takes 4 operations to make a new element, and no memory, which I think is pretty cool. Of course, I don’t kid myself that I’m the first one to think of this, but I still wanted to share it, especially since, as it turns out, I can’t use it for the problem I’m trying to solve in Project Euler.

Sage: great free mathematics software

February 1, 2008

Sage is an open source effort to replace the commercial math codes like Mathematica, Matlab, Magma with a free product of equal quality. It uses Python, my favorite language, as the glue that holds together a variety of existing open source math programs (Pari, Maxima, Numpy, Scipy, Blas, Lapack, etc.).

Such a product is a real revelation to me. I rarely have need for a symbolic math program, so that even though I have access to Maxima and Mathematica, I use them rarely enough so that when I do need them, I’ve completely forgotten the required syntax. But when I need the program, I really need it, which requires a long slog through the documentation. However, I use Python almost every minute of every day. But Python has its problems as well. Getting users to install all the various widgets required to run one of the pieces of software I write is itself very time consuming. Sage has solved all of this, as far as I can tell, by simply distributing all of the source code with the program. This makes for a rather time consuming build procedure, but one that worked for me the first time I tried it.

I’ve already mentioned how much I love Project Euler, and it’s gratifying to see Python making such a respectable showing on the statistics page there. However, I found that my versions of some very critical programs, like testing whether a number is prime, or computing the number of integer partitions, were enough slower than fast implementations, that my programs often took much longer than those written in Mathematica, for example. Sage solves all of these problems for me quite ably, and lets me use Python to write the code that drives the functions.

The graphics are worth noting as well. Sage has integrated the best Python plotting package (matplotlib) into a Firefox-driven duplicate of the Mathematica notebook into something that is both simple to use and elegant. They even have a free web portal for the program, where you can sign up for a free account and try things out.

Very slick and very, very well done. In all seriousness, I haven’t been this excited since I discovered Python for the first time.

End Links:

Miller-Rabin primality test

January 20, 2008

There’s a great snippet of python code for the Miller-Rabin primality test available. I had seen other versions before (here, here), but this was the first one that I could fully follow and use. Quite useful for one or two of the Project Euler puzzles where a Sieve of Eratosthenes took too much memory.

Python Lessons Learned From Project Euler

January 2, 2008

Much of my winter vacation has been spent playing around with the problems on Project Euler. I love crosswords and other puzzles, and I’ve posted before about how I love the MacHeist game that’s currently going on. But Project Euler is different from all the puzzles I’ve solved before.

The idea is to write short little programs to solve math problems that Euler himself was interested in (and, since Euler was interested in a huge variety of problems, there’s little risk of running out of topics). For example, Problem 1 is:

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.Find the sum of all the multiples of 3 or 5 below 1000.
This is typical of most of the PE problems: they describe a problem, give a solution to an easy version of the problem, and then pose the hard version of the problem. The easy version is useful for working out stupid bugs in ones code. Python is ideal for solving these problems (and, indeed, is the 2nd most popular language for solving Project Euler problems, after C/C++). The arbitrary length integers make even large integer problems fairly easy to tackle, and the language has good support for string manipulations, which are also quite useful.One of the qualities that sets PE apart from other puzzles is that I feel that I’m enriching my math and programming skills by solving the problems. Indeed, I find that I’m learning or reinforcing some fairly important lessons. Here’s what I’ve learned in solving the first 70 problems, in no particular order:

  • Use sets or dicts for membership testing; the sets in particular are amazingly versatile: they’re not much harder to use than a normal list, but you immediately get uniqueness, and you can always convert back by taking list(set);
  • range() is fine for small lists, but takes up a large amount of memory for billion-element lists. xrange() is trivial to use in its place;
  • It’s always good to store your intermediate results in a list such that list.sort() returns your final answer without sifting through pages of data;
  • Always solve a small version of the problem before tackling a large problem; make the solution automatic enough that by changing a single number, you can go from the easy version to the hard version.

Google Chart API!

December 6, 2007

This is freakin’ awesome. Google has come out with a Chart API, for dynamically creating charts over the internet. What is so great about this is that anyone could have thought of this: I create hundreds of plots each week, mostly with the excellent Matplotlib library in Python. Did I think of this? Nope. Did Google? Yep. Awesome.

Python, Haskell, Ruby Smackdown

November 30, 2007

Antonio Cangiano writes about how the new and improved Ruby 1.9 is actually faster than Python for a fibonacci program; Haskell Hacking writes about how Haskell blows them both away. Though I’m a Python guy, I’m delighted to see news like this come out. First, I have a great deal of regard for both Ruby and Haskell. These are remarkable languages, well-designed, easy to read, and Python can benefit a great deal from the places where they exceed it. Secondly, they’re proof of concepts that Python can be faster in these areas, and can point to ways that the cPython interpreter itself can optimize itself. It’s nice to see dynamic language do as well as these do. Haskell, in particular, since it has been my belief that functional programming and lazy evaluation are the way to go to make dynamic languages faster.

Followup: There’s a rebuttal from another Haskell blogger here, who says:

we learn the shocking fact that static, compiled, finite-precision Haskell is faster than dynamic, interpreted, arbitrary-precision Python and Ruby at a heavily numeric and recursive task … Color me not impressed.

So maybe it isn’t as big a deal that Haskell beat the other two. Still, I’m impressed with almost everything that I see from that language.