It's Not Magic
Writings of a techie wizard
Thu, 15 Dec 2011
Another Brief Nerd Interlude
Years ago, Doug McIlroy, the inventor of the Unix pipe, published a paper on techniques for computing the terms of power series. The paper talks about a number of key concepts in programming, such as "lazy" evaluation, that were not well supported by most programming languages at the time, which is why McIlroy spent a good portion of the paper describing an implementation of his techniques in a new language designed by Rob Pike.
I came across this paper recently and realized that Python's generators would be a perfect fit for representing power series. They support all the key techniques McIlroy described, particularly "lazy" evaluation (a generator doesn't compute any specific term of its series until it is asked for it in sequence). You can see the Python implementation I came up with on github here. A particularly neat feature is that you can recursively include a Python generator in itself; this allows the recursive nature of many power series to be directly represented in the code. For example, here's the exponential series:
def _exp(): for term in integral(EXP, Fraction(1, 1)): yield term EXP = PowerSeries(_exp)
No monkey business with factorials; just a generator that recursively integrates itself. This same trick also works for implementing operations on power series; for example, any series can be exponentiated by a method similar to the above. The reciprocal and inverse operations on series use similar tricks, which basically make the code look just like the mathematical descriptions of those operations in McIlroy's paper.
Once I got the Python implementation working smoothly, I began checking online to see what other recent implementations of these techniques existed, and found that McIlroy posted an implementation in Haskell on the web in 2007. All of the key operations are one-liners. This is possible because Haskell has built-in support for expressing these operations declaratively, instead of having to define functions and use for loops and so on. So in a sense, my Python implementation is a case of Greenspun's Tenth Rule. But it's still fun.
Open Source Projects
Old Open Source Projects
Copyright © 2011-2013
by Peter A. Donis
All Rights Reserved