research!rsc

These are the most elegant programs I've ever seen.

A power series is an infinite-degree polynomial. For example,

Computation and manipulation of power series has a long and distinguished history as a test of so-called stream processing systems, which manipulate (arbitrarily long finite prefixes of) infinite streams. In the 1970s, Gilles Kahn and David MacQueen used power series as an unpublished test of their coroutine-based stream-processing system. Hal Abelson and Gerald Sussman devote Section 3.5 of their well-known 1985 textbook Structure and Interpretation of Computer Programs to stream processing, although only a single exercise mentions power series. In the late 1980s, Doug McIlroy explored power series processing in the context of Rob Pike's concurrent programming language Newsqueak. He published the work in his 1990 paper “Squinting at Power Series” (the Newsqueak interpreter is named squint).

Stream processing requires lazy evaluation to keep computations finite. All of these implementations built frameworks for lazy evaluation atop other building blocks: concurrent processes (Kahn and MacQueen, McIlroy) or first-class functions (Abelson and Sussman). It's only natural, then, to consider what the implementations would look like in a language with explicit support and encouragement for lazy evaluation, a language like Haskell. McIlroy revisited the topic in the context of Haskell in his 1998 Functional Pearl, “Power Series, Power Serious” (gzipped PS).

Symbolic manipulations can represent a power series as just a list of integer coefficients. Mathematically, this is equivalent to writing the polynormial in Horner form, which repeatedly factors x out of the remainder of the polynomial: F = f₀ + xF₁, or in Haskell, (f:fs) (= f + x fs).

In this form, e^x is [1, 1, 1%2, 1%6, 1%24, 1%120, ...].

As a simple example, McIlroy's addition and multiplication implementations mimic the usual rules:

(f₀ + xF₁) + (g₀ + xG₁) = (f₀ + g₀) + x(F₁ + G₁)

  (f:fs) + (g:gs) = f+g : fs+gs

(f₀ + xF₁) (g₀ + xG₁) = f₀g₀ + x(f₀G₁ + g₀F₁) + x²(F₁G₁)

{- (f:fs) * (g:gs) = (f*g : 0) + (0 : f.*gs + g.*fs) + (0 : 0 : fs*gs) -}
   (f:fs) * (g:gs) = f*g : (f.*gs + g.*fs + (0 : fs*gs))
       c .* (f:fs) = c.*f : c.*fs

The commented-out definition is a more literal, but less efficient, translation of the equation.

Haskell's pattern matching and operator overloading are responsible for the elegance of the presentation, but it is lazy evaluation that is responsible for the elegance of the computation. The underlying lazy computation mechanism (whether provided directly, as in Haskell, or via other primitives, as in Scheme and Newsqueak) takes care of all the bookkeeping necessary to compute only the terms needed.

McIlroy implements derivative and integral operators as well. Both rely on a helper function to keep track of the leading multiplicative constant:

deriv (f:fs) = (deriv1 fs 1) where
    deriv1 (g:gs) n = n*g : (deriv1 gs (n+1))

integral fs = 0 : (int1 fs 1) where
    int1 (g:gs) n = g/n : (int1 gs (n+1))

These definitions enable elegant definitions of the power series for exp, sin, and cos:

expx = 1 + (integral expx)

sinx = integral cosx
cosx = 1 - (integral sinx)

To me, these three equations programs are beautifully elegant, almost magical.

When learning recursion in an introductory programming course (or induction in an introductory math course), it is common to feel like there's no actual work going on: one case just got rewritten in terms of others. The base cases, of course, provide the foundation on which the others are built. Learning to determine which base cases are necessary given the recursive steps is the key to being comfortable with recursion. Only then can you see that the program or proof is structurally sound.

The recursion in the definitions above gives me the same introductary queasiness, because it is hard to see the base cases. Upon closer inspection, the base case is in the expansion of integral, which begins with a constant zero term, making it possible to determine the first term in each of those series without further recursion. Having the first term makes it possible to find the second term, and so on.

For me, understanding why they work only enhances the beauty of these programs.

McIlroy continues the theme in his 2000 ICFP paper “The Music of Streams” (gzipped PS), which contains even more examples but fewer programs.

I've extracted runnable code from the paper. McIlroy maintains equivalent code. McIlroy's version uses these shorter versions of integral and deriv, which will delight Haskell aficionados:

int fs = 0 : zipWith (/) fs [1..]
diff (_:ft) = zipWith (*) ft [1..]

For the numerical computing aficionados, there is a one-line implementation of power series reversion:

revert (0:ft) = rs where rs = 0 : 1/(ft#rs)

(For comparison, Donald Knuth devotes Section 4.7 of Volume 2, about 10 pages, to the manipulation of power series. The equivalent reversion implementation takes about half a page.)

Posted by rsc on Friday, February 22, 2008 | 3 comments |
Labels: algorithms, code, coroutines, math, notation

At first glance, Duff's Device is one of the most mysterious pieces of C code you'll ever see:

    void
    send(short *to, short *from, int count)
    {
        int n=(count+7)/8;
        switch(count%8){
        case 0: do{ *to = *from++;
        case 7:     *to = *from++;
        case 6:     *to = *from++;
        case 5:     *to = *from++;
        case 4:     *to = *from++;
        case 3:     *to = *from++;
        case 2:     *to = *from++;
        case 1:     *to = *from++;
                }while(--n>0);
        }
    }

It's an 8x-unrolled while loop interlaced with a switch statement. The switch jumps into the middle of the while loop to handle the part of the count that isn't a multiple of 8. Once inside the while loop, the switch logic never executes again.

Tom Duff first described Duff's Device in a November 1983 email to Ron Gomes, Dennis Ritchie, and Rob Pike. It spread to a larger group when he revised the note and posted it to net.lang.c in May 1984. The 1984 message is the less frequently reproduced but is the one in which Duff gave it a name: “I think I'll name it after myself — ‘Duff's Device’ has a nice ring to it.” Bjarne Stroustroup used a variant as an example in The C++ Programming Language. Eventually Duff posted the 1983 message along with commentary to Usenet in 1988, in response to a heated debate about the merits of Duff's Device as an idiom.

Duff's Device is rarely worth reusing, but in the specific case that it was written for, it was an apt tool for the job. The more amazing thing about it is that it's valid C. Since a switch statement is just a computed goto and the case labels are just labels, switching into the middle of a while is just as legal as goto'ing there. In Duff's description, he wrote:

It amazes me that after 10 years of writing C there are still little corners that I haven't explored fully. (Actually, I have another revolting way to use switches to implement interrupt driven state machines but it's too horrid to go into.)

In 2000, Simon Tatham wrote about a C coroutine implementation based on this idea, to make callback-based event-driven programming easier. (He used it in his Windows SSH client, PuTTY.) A switch-based implementation has the advantage of being portable and relatively simple to emit using preprocessor macros, though it doesn't allow the use of other switch statements in the code.

A year ago, I asked Duff if that kind of coroutine implementation is what he had in mind, and he confirmed that it was, pointing at a blog comment and adding:

I only ever did this twice: once for a coroutinish application in an SDLC driver (for UNIX to IBM OS/360 communication) in the 1970s and once for loop unwinding in 1983. Both were desperate situations. I have never felt the need to stoop so low since then.

(Today, of course, there are few excuses not to use threads, or at the very least, coroutines with real stacks. If you are stuck with callbacks and events, something like Max Krohn's tame preprocessor is a cleaner solution than macros.)

Posted by rsc on Wednesday, January 30, 2008 | 4 comments |
Labels: coroutines