research!rsc

Some of the earliest mathematics was invented to design calendars. Since the Earth rotates around its axis about 365.25 times per revolution around the sun, some years in a solar calendar need to be 365 days long and others 366 days long in order to keep the calendar in sync with the sun. The key challenge in designing a calendar is to find a pattern of short and long years that achieve the desired approximation. The Julian calendar uses the simple pattern of one leap year every four years. The Islamic calendar, which syncs against the moon, needs to have years approximating 354.37 days per year. To do this, it inserts 11 leap years over a 30 year cycle.

Bresenham's algorithm chooses which pixels to use when drawing a bitmap representation of a line. The basic idea, for a mostly horizontal line that angles upward, is to move left to right, one pixel at a time, moving upward when the distance between the actual y and the pixel-approximate y has grown too large.

drawline(dx, dy)    /* mostly horizontal line, upward line */
    x = 0
    y = 0
    error = 0
    while (x < dx || y < dy)
        drawpoint(x, y);
        x++;
        if (x*(dy/dx) - y > 0.5)
            y++;
    drawpoint(x, y);

For example, drawing a line 30 pixels long and 12 pixels high produces

In their 2004 paper “Line Drawing, Leap Years, and Euclid” (PDF), Mitchell Harris and Edward Reingold show that the problem of choosing a leap year pattern is a generalized version of Bresenham's line drawing problem. This is surprising, but makes intuitive sense when you realize that both are concerned with computing an integer approximation to a rational fraction as evenly as possible. Leap year computations are basically drawing a line with slope 1/365.25 (or 1/364.37): each pixel is a day, and its y coordinate determines which year it belongs to.

Harris and Reingold go on to show that Euclid's algorithm can compute leap year patterns! Euclid's algorithm usually computes the greatest common divisor of two numbers:

gcd(u, v)
    while (u != v)
        if (u > v)
            u = u mod v;
        else
            v = v mod u;
    return u;

The connection is that if you're trying to draw a line dx by dy but dx and dy have some common divisor d, then the Bresenham algorithm will simply emit the pattern for a line dx/d by dy/d, d times. The example shown above can be viewed as two 15x6 lines, three 10x4 lines, or six 5x2 lines:

Euclid's algorithm gives the largest possible d, and you can adapt the algorithm to build the associated repeating pattern while it computes d.

In Volume 2, Donald Knuth calls it “the granddaddy of all algorithms, because it is the oldest nontrivial algorithm that has survived to the present day.” It also seems to have some famous children.

Further reading:

Bresenham described his algorithm in his 1965 paper “Algorithm for computer control of a digital plotter”

Godfried Toussaint's 2005 paper “The Euclidean algorithm generates traditional musical rhythms” shows yet another pattern connection, this time to rhythms.

The Gregorian calendar is an evolved version of the Julian calendar and can't be viewed as a Bresenham line. However, Jeffrey Shallit's 1994 paper “Pierce expansions and rules for the determination of leap years” shows how to compute leap years using Pierce expansions, and this model can generate the Gregorian calendar too.

Edward Reingold and Nashum Dershowitz are the experts on calendar algorithms. Their pair of papers 1990 Calendrical Calculations and 1993 Calendrical Calculations, II: Three Historical Calendars gives algorithms (and Lisp code) for essentially any date-related computation you might want to do. They expanded upon these papers in their books Calendrical Calculations: The Millennium Edition and Calendrical Tabulations: 1900-2200.

Posted by rsc on Friday, February 29, 2008 | 1 comment |
Labels: algorithms, graphics, math, theory

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

(Today's late post is brought to you by Continental Airlines and the letters E, W, and R.)

ITA Software's flagship product, QPX, is an airline ticket search engine. You tell it where you want to fly and when, and it comes back with possible trips at many different prices. QPX powers Orbitz, CheapTickets, Hotwire, and a half dozen airline company sites.

The airlines want to charge as much as possible for each ticket they sell, but still fill each plane. But the airlines can't just ask possible customers how much money they have to spend, and they also can't explicitly set different prices based on who is buying the ticket. So they attach restrictions to less expensive fares, hoping that business travelers won't accept the restrictions and will pay for the more expensive tickets. The most well-known example is that trips that include a Saturday night are cheaper. There are also changes in ticket prices based on the rest of the trip: a one-way ticket from Boston to Newark might cost $349, but a one-way ticket from Boston to Raleigh via the very same Newark flight only $145*. Of course, if you make a habit of buying multi-leg tickets and skipping the connecting flights, the airlines might blacklist you (or so a travel agent once told me).

These restrictions make QPX's job harder than you might think. In fact, the general case is undecidable!

ITA's Chief Scientist, Carl de Marcken, gives a talk titled “Computational Complexity of Airline Travel&rdquo (HTML, also PDF), in which he proves that:

1. Even if the airlines publish only a single fare (with every ticket a single one way PU), and all the airports in a passenger's itinerary are fixed, so that the only remaining choice is what flights to use between the consecutive airports, deciding whether there is a valid ticket is NP-hard.

2. Fixing the flights and priceable units, but leaving open the choice of fares to pay for each flight, deciding whether there is a valid ticket is NP-hard.

3. Removing bounds on the size of solutions, the general question of whether there is a ticket for a query is EXPSPACE-hard. That is, air travel planning is at least as hard (it might be harder) as deciding whether a computer program that can use space exponentially bigger than the input halts. EXPSPACE-hard problems are (thought to be) much harder than NP-complete problems like the traveling salesman problem, and even much harder than PSPACE-complete problems like playing games optimally. There is no practical hope that computers will ever be able to solve EXPSPACE-hard problems perfectly, even if quantum computing becomes a reality.

4. The final result shows that just finding out whether there is a valid solution for a query is actually harder than EXPSPACE-complete: it is unsolvable (undecidable). The question of whether a valid ticket exists can not be solved for all databases and all queries no matter how long a computer thinks. However the full proof of this result is considerably more complex than the EXPSPACE-hard proof without offering any greater understanding of the problem.

The NP-hardness proofs reduce 3-color, the EXPSPACE-completeness proof encodes a Turing machine in ticket restrictions, and the undecidability proof encodes Diophantine equations. Here's a 3x3-bit binary multiplier. (Be sure to read the early slides too, which give necessary details about ticket pricing.)

* As of February 13, Continental flights on March 12, 2008 had those prices. Actual fares may no longer match.

Posted by rsc on Wednesday, February 13, 2008 | 0 comments |
Labels: algorithms, theory, undecidability

Starting in the late 1970s, Ken Thompson proved by example that computers could completely analyze chess endgames, which involve only a small number of pieces. Thompson's key insight was that when there are only a few pieces on the board, there might be very many possible move sequences, but there are relatively few distinct board positions. An exhaustive search over board positions can run efficiently even though an exhaustive search over move sequences would not. (The same idea underlies Thompson's efficient parallel NFA simulation for regular expression matching.)

Thompson's first endgame analyses ran on boards with only four pieces, but over time he took advantage of improvements in processor speeds and storage sizes to analyze five- and six-piece endgames as well. Thompson described the programs in the ICCA Journal: his 1986 article “Retrograde Analysis of Certain Endgames” (PDF; 1MB) introduces the technique and covers five-piece endgames, and his 1996 article “6-Piece Endgames” (PDF; 1MB) describes refinements and optimizations for six pieces.

Writing about a 262-move optimal ‘rook and knight versus two knights’ game Thompson found in the 6-piece endgames, former championship chess player Tim Krabbé explained:

Playing over these moves is an eerie experience. They are not human; a grandmaster does not understand them any better than someone who has learned chess yesterday. The knights jump, the kings orbit, the sun goes down, and every move is the truth. It's like being revealed the Meaning of Life, but it's in Estonian.

On his web page, Thompson links to the online 6-piece endgame database as “Play Chess with God.” The CD art for the five-piece databases used a similar theme.

Further reading: With Joe Condon, Thompson built the chess playing computer Belle, which won the ACM North American Computer Chess Championship five times. Belle was later confiscated by U. S. Customs on its way to a contest in Moscow, on suspicion of being of military use. Thompson was quoted as saying that Belle's only military use would be “to drop it out of an airplane. You might kill somebody that way.”

Thompson's work in computer chess was featured in a Computer History Museum exhibit, which included a video of Thompson telling his story. (Beware that in the transcript, ‘endgames’ is usually written as ‘n games.’)

Posted by rsc on Friday, January 18, 2008 | 3 comments |
Labels: algorithms, Bell Labs, brute force, chess

Quicksort is a worst-case O(n²) algorithm, but if you assume some randomness in either the input or in the decisions made by the algorithm itself, the worst case becomes exceedingly unlikely and the expected runtime becomes O(n log n). As a result, most people treat quicksort as an O(n log n) algorithm. The only wrinkle in this logic is that most quicksort implementations are entirely deterministic, so whether they run in O(n log n) time depends entirely on whether the input data is “random enough.” So there must exist inputs that force these quicksorts into quadratic behavior.

Doug McIlroy's 1999 paper “A Killer Adversary for Quicksort” describes a simple way to find these deadly inputs. The basic idea is that since the C library qsort calls a user-supplied comparison function as the only way it inspects the input, the comparison function can decide what the input is “on the fly,” making decisions that are consistent with previous comparisons but bad for qsort. See the paper for the elegant details.

The paper graphs the 64-element input that forces Digital's qsort to go quadratic:

McIlroy provides a simple C program illustrating the technique. The GNU C library quicksort goes quadratic on exactly the same kind of input as the Digital quicksort.

(A note of caution for those playing along at home: the GNU C library qsort function is actually mergesort if the C library thinks the computer has enough free memory; to force quicksort, call the undeclared _quicksort function and compile with gcc -static.)

Posted by rsc on Wednesday, January 02, 2008 | 3 comments |
Labels: algorithms