Developing for Developers

New blog: Papers in Computer Science

MSDNArchive — Sun, 01 Mar 2009 01:20:00 GMT

Hey all - I apologize for the (extremely) long period of no updates, I've been prioritizing other things. I've been accepted this Fall to begin my Ph.D. at University of California, Berkeley. Since I won't be at Microsoft any longer, I've started a new blog with a new focus and more regular updates called Papers in Computer Science at Wordpress, and I've posted my first post to it:

http://papersincomputerscience.wordpress.com/

My summary from the introduction post: "Papers in Computer Science (originally called A CS Paper a Day) is a blog for abstracts and brief summaries and discussions about published computer science articles. I’ll cover a variety of areas, focus on some of the most interesting and important papers in each one, and I’ll try to explain them to a general computer science audience without specialist background. Because most of these are available online, many of them without cost, readers can consult the original paper for more information if they’re interested."

I hope anyone interested in reading more of my stuff will take a look at Papers in Computer Science and leave comments. Thank you all for your patronage!

P-complete and the limits of parallelization

MSDNArchive — Sat, 08 Sep 2007 06:11:00 GMT

We're entering an era where CPU clock speeds will soon cease to scale upwards and instead CPU manufacturers are planning to put more and more independent cores on a chip. Intel plans to release an 80-core chip within 5 years. Consequently the research community is setting their eye on the manifold barriers to writing correct and efficient highly-parallel programs. Will we be able to continue to get the same boosts in computational performance that we have been getting through raw clock speed?

You might think it's just a matter of building effective tools for producing programs that can take advantage of highly parallel platforms, and training programmers to take advantage of them. But there are complexity theory results to show that there are some problems that, under widely accepted assumptions, fundamentally cannot be parallelized, and this is the theory of P-completeness. For these problems, it may be literally impossible to get any significant benefit out of throwing additional cores at them - such problems are called inherently sequential. More to the point, however long it takes your computer to solve these problems right now, that's about how fast it's going to stay for a long time.

The complexity class NC (Nick's Class, named after Nick Pippenger by Steven Cook) is the set of problems that can be solved in polylogarithmic time (O(log^k n) time for some integer k) by a machine with polynomially many processors. For example, given a list of n numbers, the typical way to add them up on a single processor would be with a loop, as in this C# snippet:

    public static int AddList(List<int> list) {
        int sum = 0;
        foreach (int i in list) {
            sum += i;
        }
        return sum;
    }

But now suppose we had n/2 processors. Then each processor can add just two of the numbers, producing a list of n/2 numbers with the same sum. We add these in pairs as well, continuing until just one number remains, the sum. This fan-in algorithm requires only O(log n) time, placing the problem in NC (or technically the decision problem version of it).

We can use similar ideas to parallelize many algorithms. For example, the best known sequential algorithms for matrix multiplication requires time O(n^2.376) for an n × n matrix, but given O(n³) processors, we can solve it in O(log n) time: for each entry, we assign n processors to compute the products that are added to form that entry, then use the fan-in procedure above to add them together. This can still be done in O(log n) time. There are efficient parallel algorithms for other common problems too, like sorting, string matching, and algorithms in graph theory and computational geometry.

In real life of course we don't have such a polynomial number of processors, so this is in some ways an impractical definition; all it really tells us is that if we have many processors, we can solve the problem dramatically faster. Practical parallel algorithms pay more attention to speedup - how many times faster the algorithm is with k processors than with one. Effective speedup requires a good way of dynamically decomposing the problem. But speedup is difficult to characterize formally in a machine-independent way, so we'll stick with the definition of NC.

The famous P = NP problem asks whether all problems whose solutions can be verified in polynomial time, can also be solved in polynomial time. In this problem, the problems in P are considered "feasible to solve", and most computer scientists believe P ≠ NP - that is, that some problems which we can feasibly verify cannot be solved feasibly. We believe this in part because we have shown that many important problems are NP-hard. To say a problem is NP-hard means there is a polynomial time algorithm that can take any problem in NP and convert it into an instance of that problem; consequently, if an NP-hard problem can be solved in polynomial time, so can any problem in NP.

But when you have a large number of processors, the question turns around: we no longer consider P feasible, because any algorithm using just one CPU would face a dramatic slowdown compared to algorithms that take advantage of all of them. Now NC algorithms are considered feasible, and we ask: can all problems with efficient sequential solutions be highly parallelized? That is, is NC = P?

Just as with P = NP, nobody knows. However, just as there are NP-hard problems, there are P-hard problems. A P-hard problem is a problem where there is an NC algorithm (running in polylogarithmic time on a polynomial number of processors) to convert any problem in P into an instance of that problem. If we had an NC algorithm to solve a P-hard problem, it would imply that NC = P, providing a constructive and mechanical way of parallelizing any problem we can solve sequentially in polynomial time. Just as most researchers believe P ≠ NP, most researchers believe NC ≠ P, which implies that no P-hard problem has an NC algorithm to solve it.

One of the most important P-complete problems is linear programming: given a linear target function in terms of some variables, and a number of linear inequalities expressing constraints between those variables, find the optimal (maximum) value of the target function. It finds uses in operations research, microeconomics, and business management, and it would be great if we could exploit large number of cores to tackle larger instances of these problems. But no one has found an efficient parallel algorithm for this problem.

Another important example is the circuit-value problem: given a digital logic circuit and its inputs, compute the output of the circuit. There is a trivial linear time solution in the number of gates: continually get the next gate that you know both inputs of and compute its output, until you're done. It would seem that certain gates don't influence each other and be computed in parallel. But because the wires can express complex interdependencies between values and can create long chains of dependencies, there's no clear way to parallelize a general circuit. Besides obvious applications to circuit design and simulation, many problems can be solved by families of circuits, and if we can evaluate them quickly we could solve these problems. But not only has no one has found a general way of parallelizing such evaluation, they've shown that even if we restrict the circuits so that wires can't cross, or such that only the inputs can be negated, the problem remains P-complete, and so difficult to parallelize (these are the planar and monotone circuit-value problems, respectively).

The first problem ever shown to be NP-complete, as well as an important problem in practice, is the boolean satisfiability problem: given a formula using AND, OR, and NOT, with a number of variables set to either true or false, determine if there is some way of setting the variables to make the whole formula evaluate to true. The P-complete version of this problem is called Horn satifiability.

Any boolean formula can be reduced efficiently to a form where it represents a three-layer circuit: a literal is either a variable or a negated variable; a clause is the OR of some literals; and an entire formula is the AND of some clauses. A "Horn clause" is a clause where all of the variables are negated except possibly one. This provides an efficient solution to the problem, since Horn clauses can be rewritten as implications (not(x1) or not(x2) or x3 is the same thing as (x1 and x2) implies x3). Formulas made up of Horn clauses appear in the context of logical deduction systems, such as the resolution systems used to model problems in AI. Again, no one has found an effective way to parallelize finding solutions for these formulas.

In the world of P = NP, people have long been studying ways to "get around" the infeasibility of finding a polynomial-time solution for NP-hard problems, such as approximation algorithms, special cases, fixed-parameter tractable solutions, and so on. We could ask the same questions about the NC = P question: could we find an approximate solution much more quickly by using all our cores, instead of running a classical algorithm to find an exact solution on just one? Are there special cases of these problems which are parallelizable? But comparatively little investigation has been done into this area, although some examples exist (like Karger and Motwani's "An NC Algorithm For Minimum Cuts").

Finally, I should mention that P-complete is relevant not only to the study of what problems can be efficiently parallelized, but also to what problems can be solved in a small amount of space. I'll leave this for another time though.

Robin's theorem

MSDNArchive — Tue, 17 Jul 2007 00:12:00 GMT

Most computer scientists are familiar with the P = NP problem, which asks essentially whether we can verify more problems in polynomial time than we can solve. So fundamentally does complexity theory hinge on this result that the Clay Mathematics Institute has labelled it one of their seven Millennium Prize Problems, and will award $1,000,000 to the researchers that solve it.

But there are six other Millennium Prize problems that are often less accessible to computer scientists, requiring background in areas such as topology and complex analysis. One of the most important of these is the Riemann Hypothesis, which is normally stated in terms of complex zeros of the Riemann Zeta function. Here, we will describe a recent characterization of the Riemann Hypothesis more suited to a computer scientist's background based on growth of integer functions.

Consider the sum-of-divisors function σ(n) that sums up the divisors of the positive integer n. For example, σ(15) = 1 + 3 + 5 + 15 = 24. If you've studied asymptotic notation, a natural question to ask is: how fast does σ(n) grow?

It's clear that since 1 and n always divides n, σ(n) ≥ n + 1, so σ(n) is Ω(n). But σ(n) cannot be O(n), because σ(n)/n takes on arbitrarily large values. To see this, consider n equal to k primorial, the product of the first k primes. Because σ(n)/n is multiplicative, σ(n)/n will equal the product of σ(p)/p = (1+p)/p = (1 + 1/p) for each prime factor. But if you expand out this product, you get a sum including 1/p for each prime factor, and the sums of the reciprocals of the primes, also known as the harmonic series of primes, diverges. But the harmonic series of primes grows quite slowly, and is O(ln ln n).

From this we might conjecture that σ(n) grows at a rate of n ln ln n. Gronwall's theorem, proven in 1913, shows that in fact:

lim sup_n→∞σ(n) / (n ln ln n) = e^γ.

Here, e and γ are both constants (the latter is the Euler-Mascheroni constant), so e^γ is a constant, equal to about 1.781. This is somewhat stronger than the statement that σ(n) is O(n ln ln n) - the lim sup actually says that for any ε > 0, some tail of the sequence σ(n) / (n ln ln n) must be bounded above by e^γ + ε.

But instead of a lim sup, what we'd really like is to have a hard bound; we'd like to say:

σ(n) / (n ln ln n) < e^γ for all n.

This is false, and a quick search identifies 27 small counterexamples: 2, 3, 4, 5, 6, 8, 9, 10, 12, 16, 18, 20, 24, 30, 36, 48, 60, 72, 84, 120, 180, 240, 360, 720, 840, 2520, and 5040. However, no search to date has located any larger counterexamples. From this we might conjecture:

σ(n) / (n ln ln n) < e^γ for all n > 5040.

In 1984, Guy Robin showed that in fact, this statement is true if and only if the Riemann hypothesis is true. This is Robin's theorem. Thus, the Riemann hypothesis could, in theory, be proven false by exhibiting a single positive integer n > 5040 such that σ(n) > e^γ(n ln ln n). Because most mathematicians believe the Riemann hypothesis is true (just as most computer scientists believe P ≠ NP) it's more likely that the above statement is true, and a proof of this fact would win you $1,000,000.

Cache-oblivious data structures

MSDNArchive — Tue, 12 Jun 2007 23:54:00 GMT

In most data structure and algorithms classes, the model used for basic analysis is the traditional RAM model: we assume that we have a large, random-access array of memory, and count the number of simple reads/writes needed to perform the algorithm. For example, selection sort takes about n(n-1)/2 reads and 2n writes to sort an array of n numbers. This model was relatively accurate up through the 1980s, when CPUs were slow enough and RAM small enough that reading and writing the RAM directly didn't create a significant bottleneck. Nowadays, however, typical desktop CPUs possess deep cache hierarchies - at least a sizable L1 and L2 cache - to prevent runtime from being dominated by RAM accesses. Data structures and algorithms that efficiently exploit the cache can perform dramatically more quickly than those that don't, and our analysis needs to take this into account.

For example, suppose you have a list of bytes stored in both an array and a linked list. Provided the linked list nodes are allocated at random positions, iterating over the array will be dramatically faster, exhibiting optimal locality of reference, while iterating over the linked list will trigger a cache miss for every element. We can formalize this: assume a cache line - the size of the block read into the cache at one time - has size B bytes. Then traversing the linked list requires n cache misses, while traversing the array requires precisely ceiling(n/B) cache misses, a speedup of B times.

It is possible to design data structures that specifically take advantage of the cache line size B. For example, in a previous post I described an unrolled linked list data structure where each node fills up a cache line, achieving similar cache performance to arrays while still allowing efficient insertions and deletions. Such a data structure is termed "cache-aware", because it explicitly takes knowledge about the cache model into account in its construction.

Unfortunately, this approach has a compelling disadvantage: one has to tune the data structure based on the cache line size B. This is a problem for many reasons: first of all, programs, even those compiled to native machine code, are often run on many different processors with different cache line sizes. A Pentium 4 can run the same machine code as a Pentium II, but has a very different cache model. But even when running on a fixed machine, these simple cache-tuned data structures can't take advantage of multilevel caches, because each one has a different size and a different block size. One can attempt to build a data structure for a specific cache hierarchy, but this is even less robust against a change of machine.

In 1999, in his masters thesis, Harold Prokop came up with an interesting solution to this problem: the idea of a cache-oblivious algorithm. This does not mean that the algorithm does not take advantage of the cache; to the contrary, it does so quite effectively. What it means is that the algorithm does not need to know the cache line size; it works effectively for all cache line sizes B simultaneously. This allows the algorithm to robustly exploit caches across many machines without machine-specific tuning. More importantly, it allows the algorithm to exploit multilevel caches effectively: because it works for all B, it applies between each two adjacent levels of cache, and so every level is well-utilized.

In fact, the multilevel advantage extends beyond simply taking advantage of both the L1 and L2 caches - such algorithms perform well even when the storage hierarchy is expanded to encompass much slower media such as slow memory, Flash, compressed memory, disks, or even network resources, without knowledge of block sizes or access patterns. For example, cache-oblivious B-trees can be used to implement huge database indexes that simultaneously minimize disk reads and cache misses.

Note that just because the algorithm doesn't depend on B doesn't mean the analysis of the algorithm doesn't depend on B. Let's take an example: iterating over a simple array. As noted earlier, this requires about n/B cache misses, which is optimal. But neither the array's structure nor the iteration algorithm explicitly take B into account. Consequently, it works for all B. In a multilevel cache, all data prefetched into every cache is used. This is the theoretical benefit of cache-oblivious algorithms: we can reason about the algorithm using a very simple two-level cache model, and it automatically generalizes to a complex, deep cache hierarchy with no additional work.

While cache-oblivious algorithms are clearly useful, at first it's not clear that there even exist any other than simple array iteration. Thankfully, extensive recent research has revealed cache-oblivious data structures and algorithms for a multitude of practical problems: searching (binary trees), sorting, associative arrays, FFT computation, matrix multiplication and transposition, priority queues, shortest path, query processing, orthogonal range searching, and more emerging every year. Practical comparison studies have shown a significant performance gain on real processors compared to traditional algorithms, although carefully-tuned cache-aware algorithms still enjoy an advantage on the specific machines they are tuned for (see e.g. [1][2]).

Many of the cache-oblivious data structures and algorithms that have been published are relatively complex, but here I'll describe a simple one just to give you a feel for it. This discussion is based on parts of this paper.

Consider the following simple search problem: we have a static list of records, and we wish to find the record with a particular key. Traditionally, this problem is solved with either an array and binary search, or a binary search tree. Both of these approaches exhibit dismal cache behavior. Database applications rely on B-trees, which group several records into each block. This greatly improves performance, but requires knowledge of the block size, and does not work across all levels of the cache hierarchy.

The key is the van Emde Boas layout, named after the van Emde Boas tree data structure conceived in 1977 by Peter van Emde Boas. Suppose for simplicity that the number of elements in a power of 2. We create a complete binary tree containing our records (where by "complete" we mean that all leaves are at the same level). This tree will have height h = log₂n.

Take a look at the first h/2 levels of the tree. This subtree has 2^h/2 leaves, each itself the root of a tree of height h/2. In the van Emde Boas layout, we first lay out the subtree of height h/2 rooted at the root, followed by the subtrees rooted at each leaf of this tree from left to right. Each subtree is recursively laid out using the van Emde Boas layout. An example diagram is shown in Figure 1 of the paper.

The analysis proceeds like this: at each step of the recursion, the size of the subtrees being laid out is the square root of the size of the subtrees laid out at the previous step. Consequently, as we recurse, at some point we will be laying out subtrees between sqrt(B) and B in size (let's call it C). Each of these subtrees can be retrieved in a single memory transfer, and this layout partitions the tree into subtrees of this size.

The search procedure works like this: we access the root, which pulls in the first log₂C levels of the tree. We have only cache hits until we access a leaf of the tree, which pulls in the next log₂C levels of the tree rooted at that leaf. This continues until we reach the bottom of the tree. Because only 1 out of every log₂C accesses causes a cache miss, the total number of misses is log₂n/log₂C = log_Cn. Because C is at least sqrt(B), log₂C is at least (log₂B)/2, and the total number of misses is at most 2log_Bn.

We made one small assumption above: that each subtree is aligned with the cache block boundaries; if this is not true, each subtree can require two accesses, bringing the final maximum to 4log_Bn. Probabilistic analysis shows that with a random starting point in memory and sufficiently large block size, the real constant is closer to 2.

This is the same asymptotic performance as B-trees, and is faster than ordinary binary trees by a factor of about log₂n/2log_Bn = (log₂B)/2. For example, for a disk block size of 2048 elements, this would be a speedup of about 6.5 times. The advantage is more than theoretical: to quote the paper, "preliminary experiments have shown, surprisingly, that CO [cache-oblivious] B-trees can outperform traditional B-trees, sometimes by factors of more than 2".

But as practical as the research is in cache-oblivious algorithms, many applications and libraries have yet to take advantage of them. The data structures supplied with .NET, Java, Lisp, and so on are not cache-oblivious. We need to start putting this research into practice and reaping the benefits. I hope this gives you some insight into cache-oblivious algorithms, and I hope you will consider taking advantage of it the next time you develop a product whose performance critically depends on memory access patterns.

K-nearest neighbor spatial search

MSDNArchive — Fri, 08 Jun 2007 03:23:00 GMT

I apologize to everyone for the hiatus - I realise a post is more than a little overdue and will try to be more regular in the future. I appreciate the support that I've received for the blog and continuing it.

Consider the following problem: you're starting a website for a chain of restaurants, and you want to have a "Locations" feature where a person can enter their address and get a list of the 10 outlets closest to their location, in order of distance from their location (as the crow flies, for simplicity). Moreover, each time they hit the "next" button, they should get the next 10 locations in order of distance.

More abstractly, say you have a list of locations on a 2D map, each specified as a point with an x and y coordinate. You can do any sort of preprocessing you like of this list of points. Then I ask you to efficiently answer: given some point, which is not necessarily one of the points in the list, produce a list of the k points closest to that point, in order of their distance from that point.

This problem, called the k-nearest neighbor problem, is very well-studied and has a variety of solutions. Here we describe one very simple but effective incremental algorithm first proposed by Hjaltason and Samet in their "Ranking in Spatial Databases" in 1995 (Citeseer). A spatial database is a database specialized to perform this type of query over its (potentially very large) data set.

Consider the following very simple but inefficient algorithm: we take each point, one at a time, and add it to a priority queue with its key equal to the distance from the search point. We then extract minimum repeatedly to obtain the points in order of their distance from the search point. Once we've extracted the first k, we're done, and with a good data structure this requires only O(k log n) time. Unfortunately, even with the best priority queue data structures based on Fibonacci heaps, adding all n points requires amortized linear (O(n)) time and space, which for a very large number of points as one would expect to find in a spatial database is prohibitive.

One approach to narrowing down the points we have to examine is to split the space into larger regions. The idea is that points close to the search point should be found in regions close to the search point. In more detail, we first divide the space into two rectangular regions, each containing some of the points. Each of these regions is, in turn, subdivided into two rectangular regions, and we continue in this manner until all the resulting regions have some small constant number of points in them. The result is a tree of regions, where each node is a subregion of its parent node, and each point lies in one of the leaves. It really doesn't matter how we choose to divide up the space - the authors used PMR quadtrees, but you could use kd-trees, simple quadtrees, or whatever you like. All that matters is that in the end each leaf region contains only a small number of points.

Creating the regions requires examining all the points, but the key is that this part does not depend on the search point, and so only has to be done once in preprocessing. The regions may require rebalancing as points are added or removed, but this can be done efficiently as well.

However, there's a trick to this: if the search point lies near the boundary of a region, its closest points may actually lie in a different, adjacent region. How can we correctly locate the nearest point without spuriously gathering all the points in all the regions around the point?

Here's the key insight: our priority queue will contain two types of objects, single points and regions. The key of each point is set to its distance from the search point, while the key of each region is set to the minimum distance of that region from the search point. Each time we extract a point from the queue, we output it; each time we extract a region from the queue, we expand it by adding its subregions to the queue, or if it's a leaf region, the points in that region. Because the region's key is set to the minimum distance of that region from the search point, any point that is outputted before the region is expanded is closer to the search point than any point in the region. This ensures correctness.

For example, consider the search problem shown to the right. There are 9 points, and we wish to find the 5 closest points to the big red search point in order. The space is recursively subdivided into regions as shown in the tree.

For the purpose of discussion we will refer to each region by the set of points it contains. Initially, the queue contains only the large region containing all points.

Current queue in order: {1,2,3,4,5,6,7,8,9}

This contains the red search point, so its minimum distance from the search point (its key) is zero. We remove it from the queue and add its two children {1,3,4,8} and {2,5,6,7,9}. The region {2,5,6,7,9} again contains the red search point and so has key zero. The minimum distance to the region {1,3,4,8} is shown by the red line extending out of the red point and to the left and is larger than zero.

Current queue in order: {2,5,6,7,9}, {1,3,4,8}

We extract and expand {2,5,6,7,9}, getting leaf region {2} and region {5,6,7,9}. The distance to region {2} is shown by the red line extending upward from the red point, which is clearly longer than the distance to region {1,3,4,8}. The region {5,6,7,9} has distance zero.

Current queue in order: {5,6,7,9}, {1,3,4,8}, {2}

We extract and expand {5,6,7,9} to get the regions {5,9} and {6,7}. The distance to region {6,7} is shown by the red line extending right from the red point, the longest distance so far encountered. The region {5,9} contains the red point and has distance zero.

Current queue in order: {5,9}, {1,3,4,8}, {2}, {6,7}

Next, we extract and expand region {5,9}. This is a leaf region, so we add the points 5, 9 to the queue. Point 5 is only slightly closer to the red point than region {2}, but farther than region {1,3,4,8}. Point 9 is farther than everything so far encountered.

Current queue in order: {1,3,4,8}, 5, {2}, {6,7}, 9

We extract and expand {1,3,4,8} to get {1,3,4}, at the same distance as {1,3,4,8}, and {8}, the distance of which is shown by the diagonal red line.

Current queue in order: {1,3,4}, {8}, 5, {2}, {6,7}, 9

We extract and expand {1,3,4}, which is a leaf region, and add the three points. Point 4 is closest, point 3 lies between {6,7} and 9 in distance, and point 1 is slightly closer than point 9.

Current queue in order: 4, {8}, 5, {2}, {6,7}, 3, 1, 9

We extract and output point 4, the closest point. We expand region {8} and add point 8, which is slightly farther than point 3.

Current queue in order: 5, {2}, {6,7}, 3, 8, 1, 9

We extract and output point 5. We expand region {2} and add point 2, which is slightly closer than point 8.

Current queue in order: {6,7}, 3, 2, 8, 1, 9

We extract and expand region {6,7}, adding points 6 and 7. Point 6 is farthest yet and point 7 even farther.

Current queue in order: 3, 2, 8, 1, 9, 6, 7

We've now expanded all regions and the remaining points in the queue are the remaining points not yet outputted in order.

The nice thing about this particular algorithm is that it's fully incremental: we could have stopped at any time once we had enough points, and pick up later when we want more. This makes it quite appealing in database queries, where you don't really know how many points you need in advance, because you might want to filter the results according to additional attributes. In practice, in a typical spatial database, one only needs to expand a small fraction of the regions in order to find a small number of closest points. The idea can be extended into higher-dimensional regions, spatial objects other than points such as areas and line segments, and even more general problems where one can easily compute the minimum key of a subtree of objects.

I hope you found this informative and I appreciate your comments.

Functional list processing in C# 2.0 with anonymous delegates

MSDNArchive — Sat, 01 Jul 2006 01:16:00 GMT

One of the benefits of functional languages is their great flexibility in list manipulation, which enables them to express certain computations concisely that would require one or more verbose loops in procedural languages. Many of the features that functional programmers take for granted such as first class functions, persistent data structures, and garbage collection were included in Lisp partly for the purpose of convenient list processing (Lisp stands for "LISt Processing").

Although C# is not a functional language, in C# 2.0 and the accompanying .NET Framework 2.0, a number of constructs and library calls have been added that allow generic List objects to be manipulated in a more functional way. The most important of these is anonymous delegates, which are the same construct that functional programmers call closures. We start out with the following method, which removes all elements from a list that match a given condition, passed in as a delegate:

    delegate bool Predicate<T>(T arg);

    static void RemoveAll<T>(List<T> list, Predicate<T> condition)
    {
        int i=0;
        while (i < list.Count)
        {
            if (condition(list[i]))
                list.RemoveAt(i);
            else
                i++;
        }
    }

The use of the word "predicate" here comes from the mathematical predicate, which is just a function that yields a boolean true/false value.

Now say we have a list of integers, and we want to remove all the multiples of x, where x is a given number. In C# 1.0, using RemoveAll for this is so no easy task, requiring a temporary object:

    class RemoveMultiplesOfHelper
    {
        private int x;
        public RemoveMultiplesOfHelper(int x)
        {
            this.x = x;
        }
        public bool IsMultiple(int n)
        {
            return (n % x) == 0;
        }
    }

    static void RemoveMultiplesOf(List<int> list, int x)
    {
        RemoveAll<int>(list, new RemoveMultiplesOfHelper(x).IsMultiple);
    }

The advantages of reusing RemoveAll are important: it's tricky to iterate through a list while removing elements from it, and programmers often introduce subtle bugs while attempting to do so. But this code has serious problems too: it's a hassle to write, it's verbose, the "temporary" class isn't very reusable and is difficult to name, and the code is spread all over the place, making the flow of execution difficult to follow. Because of these issues, many procedural programmers opt to reinvent the wheel instead.

Anonymous delegates solve all these problems. An anonymous delegate allows you to write code right in the middle of a method which has access to all the local variables in that method, and then package it up and pass it to another method to execute later. Here's our example using anonymous delegates:

    static void RemoveMultiplesOf(List<int> list, int x)
    {
        RemoveAll<int>(list, delegate(int y) { return (y % x) == 0; });
    }

In this much briefer implementation, the delegate "captures" the local value of x at the instant of its creation, then later RemoveAll uses it in computations.

Another handy trick you can do with anonymous delegates is take a method signature that doesn't quite fit an existing delegate and translate it into one that does. For example, suppose you're sorting a List using the overload of the Sort method that takes a System.Comparison delegate. Now say you're sorting a list of strings, and you have a string comparison method that takes a boolean determining whether the comparison is case-sensitive:

    delegate int Comparison<T>(T left, T right);
    void Sort<T>(Comparison<T> comparison);
    int StringCompare(string left, string right, bool caseSensitive);

At first you appear to be in quite a pickle, as the signature of the StringCompare method does not exactly match the delegate signature Comparison<string>. Anonymous delegates make it easy to overcome this problem:

    static void SortStrings(List<string> a, bool caseSensitive)
    {
        a.Sort(delegate(string left, string right) {
            return StringCompare(left, right, caseSensitive);
        });
    }

We simultaneously match the desired delegate signature and push the decision of whether to sort case-sensitively or not to the caller, where it belongs, all with very little new code.

The advantages of anonymous delegates make functions like our Sort and RemoveAll considerably more useful. Consequently, a number of such functions were included with the generic List class in .NET 2.0:

RemoveAll(Predicate<T>): Like our RemoveAll implementation above, removes all elements from a list matching a predicate.
Sort(Comparison<T>): Mentioned above, this is a Sort method that takes a delegate. Sorting with anonymous delegates is very handy for unusual sorting orders; for example, the following sorts a list of integers in reverse order:
```
list.Sort(delegate(int x, int y) { return y.CompareTo(x); });
```
Find(Predicate<T>): Finds the first element in a list satisfying a predicate. Handy for replacing linear searches. For example, this single line of code identifies the first file in the current directory modified in the last hour, or yields null if none:
```
new List<string>(Directory.GetFiles(Directory.GetCurrentDirectory())).Find(
    delegate(string path) { return File.GetLastWriteTime(path) >=
                                   DateTime.Now - new TimeSpan(1, 0, 0); });
```
Exists(Predicate<T>): Determines whether a list contains an element satisfying a predicate. Besides its added readability in situations where you only care about the existence of a value, it also has another use: since "Find" returns the default value for the type T to indicate "not found", you need to use Exists for predicates that might match null or zero, such as:
```
bool containsEven = list.Exists(delegate(int x) { return (x % 2) == 0; });
```
TrueForAll(Predicate<T>): Similar to Exists, but determines whether the predicate is satisfied by all elements of the list. For example, this line of code might run some tests and see if they all pass:
```
if (tests.TrueForAll(delegate(Test t) { t.Run(); return t.SuccessCode == 0; })
```
FindAll(Predicate<T>): Creates a new list containing all elements satisfying a predicate. In functional languages this is sometimes called select or filter, as we're "filtering out" the parts of the list not satisfying the predicate. Here's an example which, in a single line, gets a list of all processes using more than k threads:
```
List<Process> x = new List<Process>(Process.GetProcesses()).FindAll(
    delegate(Process p) { return p.Threads.Count > k; });
```
FindIndex(Predicate<T>): Just like Find, except that instead of returning the first value that satisfies the predicate, it returns its index in the list. For example, this expression yields the index of the first alphabetical character in a string:
```
new List(s.ToCharArray()).FindIndex(Char.IsLetter)
```
FindLastIndex(Predicate<T>): Same as FindIndex, except that it finds the last occurrence satisfying the predicate.
ForEach(Predicate<T>): Simply executes the given delegate on each member of the list. Although just using a foreach loop is cleaner than using ForEach with an anonymous delegate, ForEach is still handy if you already have a named delegate, as in this example that deletes all the files in the current directory:
```
new List<string>(Directory.GetFiles(Directory.GetCurrentDirectory())).ForEach(File.Delete);
```
ConvertAll(Predicate<T>): Perhaps the most powerful of the delegate-based List methods, ConvertAll feeds each element of the list through a conversion function, producing a new list. In functional languages this is often called map because, like a mathematical function, it "maps" one set to another. Here's a simple example which produces a list of strings that are the lowercased versions of the strings in the original input list:
```
list.ConvertAll(delegate(string s) { return s.ToLower(); });
```

What's more, besides being individually useful, complex transformations can be achieved by chaining these methods together in clever ways. For example, here's two lines of code that produce a list of all files larger than a given size in a given directory:

    static List<string> GetBigFiles(string directory, int bigLength)
    {
        List<string> paths = new List<string>(Directory.GetFiles(directory));
        return paths.ConvertAll<FileStream>( File.OpenRead )
                    .FindAll( delegate(FileStream f) { return f.Length >= bigLength; } )
                    .ConvertAll<string>( delegate(FileStream f) { return f.Name; } );
    }

The first ConvertAll opens all the files by using the static library method "File.OpenRead" as the converter. The FindAll filters out just the FileStream objects corresponding to large files. The final ConvertAll extracts the filenames from each stream. Each of these "list filters" is independently reusable. At first this way of programming may seem unusual or cumbersome, but just compare the size of the above method to a procedural implementation.

Don't forget also that you can write your own methods taking delegates to reap more benefits from anonymous delegates. A frequently useful application is for error handling and debugging code. Suppose you find yourself writing several methods that all look like this:

    void foo()
    {
        Debug.Write("entering foo()");
        try
        {
            // Do some stuff
        }
        catch (Exception e)
        {
            // Do some stuff
            throw;
        }
        Debug.Write("exiting foo()");
    }

It's not clear how to factor out the two "Do some stuff" parts, since they're nested right in the middle of the method. Once again, anonymous delegates are the answer. We can write a method like this taking delegates for the parts that get filled in:

    delegate void DoAction();
    delegate bool ProcessException(Exception e);
    void DebugInvoke(string name, DoAction action,
                     ProcessException processException)
    {
        Debug.Write("entering " + name);
        try
        {
            a();
        }
        catch (Exception e)
        {
            if (processException(e))
            {
                throw;
            }
        }
        Debug.Write("exiting " + name);
    }

Now we can mix-and-match try bodies and catch bodies at will, and get the tracing for free:

    void foo(int x, int y)
    {
        int z = 0;
        DebugInvoke("addFoo", delegate() { z = x + y; }, OverflowHandler);
        DebugInvoke("addBar", delegate() { z = z + x; },
                    delegate(Exception e) { Debug.Assert(false); return true; });
        DebugInvoke("addBaz", delegate() { y = z + z; }, OverflowHandler);
    }

    static bool OverflowHandler(Exception e)
    {
        if (e is OverflowException)
        {
            Debug.Write("Encountered overflow");
            return true;
        }
        return false;
    }

Unfortunately, anonymous delegate syntax, while much more concise than the alternative, is still a bit verbose, requiring not only the word "delegate" but full argument lists. One of the goals of the next version of C# is to encourage the use of these functional idioms, partly by introducing a new closure syntax that will take advantage of type inference to overcome this syntactic overhead.

Just remember to keep this rule of thumb in mind: if you're iterating over a list, and the body of the foreach is pretty small, chances are you don't need to. See if you can do the same thing more concisely and with less risk using a combination of the above methods and anonymous delegates (and amaze your coworkers). I hope you found this article helpful.

Factoring large numbers with quadratic sieve

MSDNArchive — Mon, 19 Jun 2006 22:13:00 GMT

Today I'm going to talk about how the quadratic sieve factoring algorithm works, giving a comprehensive description assuming knowledge of only basic university-level mathematics.

The foundation of the most popular public-key cryptography algorithm in use today, RSA, rests on the difficulty of factoring large integers. When keys are generated, efficient algorithms are used to generate two very large prime numbers and multiply them together. The person who generated the keys knows these two numbers, but everyone else only knows the product. The product contains enough information to encrypt a message to the person; the two primes allow the recipient to decrypt it. There is no known way to decrypt it without using the primes, but by factoring, we can extract the two prime factors from the product and break the encryption.

At the time that RSA was invented in 1977, factoring integers with as few as 80 decimal digits was intractable; all known algorithms were either too slow or required the number to have a special form. This made even small, 256-bit keys relatively secure. The first major breakthrough was quadratic sieve, a relatively simple factoring algorithm invented by Carl Pomerance in 1981, which can factor numbers up to 100 digits and more. It's still the best known method for numbers under 110 digits or so; for larger numbers, the general number field sieve (GNFS) is now used. However, the general number field sieve is extremely complicated, and requires extensive explanation and background for even the most basic implementation. However, GNFS is based on the same fundamental ideas as quadratic sieve, so if factoring the largest numbers in the world is your goal, this is the place to start.

We'll begin by addressing a few problems that at first glance have nothing to do with factoring, then assemble them into a working algorithm. I won't be motivating them first, but trust me - they're important.

Finding a subset of integers whose product is a square

Suppose I give you a set of integers and I ask you to find a subset of those integers whose product is a square, if one exists. For example, given the set {10, 24, 35, 52, 54, 78}, the product 24×52×78 is 97344 = 312². The brute-force solution, trying every subset, is too expensive because there are an exponential number of subsets.

We'll take a different approach based on prime factorizations and linear algebra. First, we factor each of the input numbers into prime factors; for now we will assume that these numbers are easy to factor. For the above example set, we get:

10 = 2 × 5
24 = 2³ × 3
35 = 5 × 7
52 = 2² × 13
54 = 2 × 3³
78 = 2 × 3 × 13

When you multiply two numbers written as prime factorizations, you simply add the exponents of the primes used. For example, the exponent of 2 in 24×52×78 is 6, because it's 3 in 24, 2 in 52, and 1 in 78. A number is a square if and only if all the exponents in its prime factorization are even. Suppose we write the above factorizations as vectors, where the kth entry corresponds to the exponent of the kth prime number. We get:

[1 0 1 0 0 0]
[3 1 0 0 0 0]
[0 0 1 1 0 0]
[2 0 0 0 0 1]
[1 3 0 0 0 0]
[1 1 0 0 0 1]

Now, multiplying numbers is as simple as adding vectors. If we add rows 2, 4, and 6, we get [6 2 0 0 0 2], which has all even exponents and so must be a square. In more familiar terms, we want the last bit of each entry in the sum to be zero. But in this case, we don't need to store all of the numbers above, only the last bit of each number. This gives us the following:

[1 0 1 0 0 0]
[1 1 0 0 0 0]
[0 0 1 1 0 0]
[0 0 0 0 0 1]
[1 1 0 0 0 0]
[1 1 0 0 0 1]

Moreover, since we're only interested in last bits, we can perform all our addition using one-bit integers with wraparound semantics (in other words, mod 2). If we add rows 2, 4, and 6 in this way, we get [0 0 0 0 0 0 0], the zero vector. In fact, all squares correspond to the zero vector.

Let's rephrase this as a matrix equation problem. If we transpose the above matrix, so that rows become columns, we get this:

[1 1 0 0 1 1]
[0 1 0 0 1 1]
[1 0 1 0 0 0]
[0 0 1 0 0 0]
[0 0 0 0 0 0]
[0 0 0 1 0 1]

Call this matrix A. If we multiply A by the vector [0 1 0 1 0 1], using one-bit integer arithmetic, we get the zero vector. This tells us precisely which numbers we need to multiply to get a square. So, our goal is to find a nonzero vector x such that Ax=0 (remember, all arithmetic here is with one-bit integers).

If you've had a course in linear algebra, this problem should look very familiar; it's the problem of finding the null space of a matrix, the set of vectors such that Ax=0. The problem can be solved using row reduction (Gaussian elimination). We row reduce the matrix, and then assign values to the free variables in a way that gives us a nonzero solution. The other variables will be determined by these values and the matrix. You probably studied this problem using rational numbers, not one-bit integers, but it turns out row reduction works just as well for these. For example, if we add row 1 to row 3 in the above matrix, we get the following:

[1 1 0 0 1 1]
[0 1 0 0 1 1]
[0 1 1 0 1 1]
[0 0 1 0 0 0]
[0 0 0 0 0 0]
[0 0 0 1 0 1]

Completing the row reduction, we eventually end up with this matrix:

[1 0 0 0 0 0]
[0 1 0 0 1 1]
[0 0 1 0 0 0]
[0 0 0 1 0 1]
[0 0 0 0 0 0]
[0 0 0 0 0 0]

If we turn this back into a system of equations and rearrange, we get this:

x₁ = 0
x₂ = −x₅ − x₆
x₃ = 0
x₄ = −x₆

Suppose we choose x₅=0, x₆=1. From the above equations, it follows that the first four vectors have the values 0, 1, 0, and 1 (remember, one-bit integer arithmetic). This gives us our final vector, [0 1 0 1 0 1]. If we were to choose x₅=1 and x₆=0 instead, we'd get a different solution: [0 1 0 0 1 0], corresponding to 24×54 = 1296 = 36².

Moreover, a theorem of linear algebra tells us precisely how many input numbers we need to guarantee that a square can be found: as long as we have more columns than rows, the null space is guaranteed to be nontrivial, so that we have a nonzero solution. In other words, we just need more numbers than prime factors used by those numbers. As this case shows, though, this isn't a necessary condition.

The one remaining problem with this method is that if one of the numbers in our set happens to have very large factors, our matrix will have a large number of rows, which requires a lot of storage and makes row reduction inefficient. To avoid this, we require that the input numbers are B-smooth, meaning that they only have small factors less than some integer B. This also makes them easy to factor.

Fermat's method: factoring using a difference of squares

You might be wondering what squares have to do with factoring. The connection is the very simple factorization method known as Fermat's method. Although not efficient in general, it embodies the same basic idea as quadratic sieve and works great for numbers with factors close to their square root.

The idea is to find two numbers a and b such that a²−b² = n, the number we wish to factor. If we can do this, simple algebra tells us that (a+b)(a−b) = n. If we're lucky, this is a nontrivial factorization of n; if we're not so lucky, one of them is 1 and the other is n.

The concept behind Fermat's algorithm is to search for an integer a such that a²−n is a square. If we find such an a, it follows that:

a²−(a²-n) = n

Hence we have a difference of squares equal to n. The search is a straightforward linear search: we begin with the ceiling of the square root of n, the smallest possible number such that a²−n is positive, and increment a until a²−n becomes a square. If this ever happens, we try to factor n as (a−sqrt(a²−n))(a+sqrt(a²−n)); if the factorization is trivial, we continue incrementing a.

Here's an example of Fermat's method from Wikipedia. Let n=5959; a starts out at 78. The numbers 78²−5959 and 79²−5959 are not squares, but 80²−5959=441=21². Hence (80-21)(80+21)=5959, and this gives the nontrivial factorization 59×101=5959.

The reason Fermat's method is slow is because simply performing a linear search of all possible a hoping that we'll hit one with a²−n square is a poor strategy − there just aren't that many squares around to hit. A better way of going about it is to proactively compute an a having this property (actually a similar property).

The key is to notice that if we take a number of a²−n values, none of which are squares themselves, and multiply them, we may get a square, say S. Let A be the product of the corresponding values of a. Basic algebra shows that A²−S is a multiple of n. Hence, (A−sqrt(S))(A+sqrt(S)) is a factorization of some multiple of n; in other words, at least one of these shares a factor with n. By computing the greatest common divisor (GCD) of each with n using Euclid's algorithm, we can identify this factor. Again, it may be trivial (just n itself); if so we try again with a different square.

All that remains is, given a list of a²−n values, to find a subset whose product is a square. But this is precisely an instance of the problem discussed in the last section. Unfortunately, recall that that the method we came up with there is not efficient for numbers with large factors; the matrix becomes too large. What do we do? We simply throw away numbers with large factors! Theoretical results show that there are a fairly large number of values in the sequence a²−n that are smooth (recall that smooth numbers have only small factors). This gives us a new factoring method that works pretty well up to a point.

For example, consider the number 90283. If we start a at 301 and increment it up to 360 while computing a²−n, we get the following values:

318, 921, 1526, 2133, 2742, 3353, 3966, 4581, 5198, 5817, 6438, 7061, 7686, 8313, 8942, 9573, 10206, 10841, 11478, 12117, 12758, 13401, 14046, 14693, 15342, 15993, 16646, 17301, 17958, 18617, 19278, 19941, 20606, 21273, 21942, 22613, 23286, 23961, 24638, 25317, 25998, 26681, 27366, 28053, 28742, 29433, 30126, 30821, 31518, 32217, 32918, 33621, 34326, 35033, 35742, 36453, 37166, 37881, 38598, 39317

None of these are squares (the first square occurs at a=398); however, if we factor each value we will discover that 7 of these values have no factor larger than 43:

6438, 10206, 16646, 19278, 19941, 30821, 35742

If we take these 7 values and feed them to the algorithm described in the last section, it finds a square: 19278×19941×30821×35742 = 423481541612104836 = 650754594². The corresponding original a were 331, 332, 348, and 355, and their product is 13576057680. Now, we can factor the number:

(13576057680−650754594)(13576057680+650754594) = 12925303086 × 14226812274 is a multiple of 90283
GCD(90283, 12925303086) = 137
GCD(90283, 14226812274) = 659
137 × 659 = 90283.

Making it faster: sieving for smooth numbers

The factorization algorithm above is considerably better than Fermat's algorithm, but if we try to scale up the size of number we factor, we quickly encounter a bottleneck: finding the smooth numbers in the sequence. Only 7 of the 60 values we computed in our last example were 43-smooth (actually we were lucky to get a square with so few vectors). As the size of the number that we're factoring grows, so does the size of the numbers in the sequence, and the proportion of smooth numbers rapidly shrinks. Although finding smooth numbers doesn't require completely factoring every number in the sequence (we only have to test primes up to the smoothness limit), it's still too expensive to test every number in the sequence this way.

The key is to observe that the prime factors of the sequence a²−n follow a predictable sequence. Let's take a look at the prime factorizations of the first ten or so numbers in our example sequence above:

318 = 2×3×53
921 = 3×307
1526 = 2×7×109
2133 = 3³×79
2742 = 2×3×457
3353 = 7×479
3966 = 2×3×661
4581 = 3²×509
5198 = 2×23×113
5817 = 3×7×277

The most obvious pattern is that every other number is even, beginning with the first one. This should be no surprise, since we're effectively adding 2a+1 to get each new number, which is always odd. Also, you'll notice that the first and second numbers are divisible by 3, as are the fourth and fifth, the seventh and eigth, and so on. If you look at the larger list, you'll notice similar patterns for larger primes; for example, the 3rd and 6th numbers are divisible by 7, and every 7th number after each of them as well. And, mysteriously, not one number in our entire sequence is divisible by 5!

So what's going on? The answer involves what number theorists call quadratic residues. A number a is called a quadratic residue mod p if there is some square S such that S−a is divisible by p. Half of all numbers are quadratic residues mod p, regardless of the value of p, and there's a simple formula for determining whether or not a particular number is: just take a, raise it to the power (p−1)/2, and then take the remainder after division by p. Then a is a quadratic residue mod p if and only if the answer is 1. Although this computation seems to involve very large values, in fact we can compute it quite quickly using exponentiation by squaring combined with frequent remainder operations.

This explains why none of our values are divisible by 5. If we compute 90283^(5-1)/2 mod 5, we get 4, which is not 1 (remember that 90283 is our original n to be factored). Thus, there is no square such that S−n is divisible by 5; but all numbers in our sequence have this form. In practice, this means we can compute just once ahead of time which factors may occur in the sequence (primes p such that n is a quadratic residue mod p), and ignore all others.

For our next mystery, why is it that given a number in the sequence divisible by p, every pth number after that is also divisible by p? Well, simple algebra shows that if a²−n=kp, then:

(a+p)²−n = (a²−n)+p(2a+p) = kp+p(2a+p).

But this doesn't explain why it always seems to be the case that there are exactly two different initial values of a such that a²−n is divisible by p (with the exception of p=2). For example, in our sequence above the 3rd and 6th values were divisible by 7. The answer again is quadratic residues: it can be shown that the modular equation x²≡y (mod p) has exactly two solutions (if it has any), and in fact there is an efficient algorithm for computing these two solutions called the Shanks-Tonelli algorithm. I won't go into it here since it requires some background in number theory, but for small primes it isn't really needed; it suffices to test the first p numbers to see which are divisible by p. For larger primes, it becomes important to avoid this expensive scan.

Recall the Sieve of Eratosthenes, an algorithm for locating prime numbers. It starts with a list of numbers, then crosses off all numbers not divisible by 2 except 2, then does the same for 3, 5, and so on until it's done. The numbers that remain must be prime. When attempting to find a list of prime numbers, this strategy is much more efficient than running even the most advanced primality test on each number individually.

We take a similar strategy here: we begin with a table of the original values in the sequence. We then visit all the numbers divisible by 2 and divide out a factor of 2. We do the same for each power of 2 up to the size of the sequence. We then do the same for every other prime up to our smoothness bound (43 in our example). In the end, the smooth numbers and only the smooth numbers will have become 1. Since we visit less and less list elements as the prime factor increases, the overall work is much less. For example, here's our original list from the above example:

We visit elements 1, 3, 5, and so on, dividing out 2. Here's the list after this first pass is complete:

159, 921, 763, 2133, 1371, 3353, 1983, 4581, 2599, 5817, 3219, 7061, 3843, 8313, 4471, 9573, 5103, 10841, 5739, 12117, 6379, 13401, 7023, 14693, 7671, 15993, 8323, 17301, 8979, 18617, 9639, 19941, 10303, 21273, 10971, 22613, 11643, 23961, 12319, 25317, 12999, 26681, 13683, 28053, 14371, 29433, 15063, 30821, 15759, 32217, 16459, 33621, 17163, 35033, 17871, 36453, 18583, 37881, 19299, 39317

Here it is after dividing out the prime factors 3, 5, 7, 11, 13, and 17:

53, 307, 109, 79, 457, 479, 661, 509, 2599, 277, 1073, 7061, 61, 163, 263, 3191, 1, 10841, 1913, 577, 6379, 1489, 2341, 2099, 2557, 1777, 1189, 5767, 2993, 18617, 1, 23, 10303, 1013, 1219, 22613, 3881, 163, 12319, 2813, 619, 26681, 4561, 1039, 2053, 9811, 5021, 37, 103, 10739, 16459, 1601, 1907, 35033, 851, 12151, 18583, 1403, 919, 39317

We see a couple 1s have already appeared; these are 17-smooth numbers. When we get all the way up through 43, we have:

53, 307, 109, 79, 457, 479, 661, 509, 113, 277, 1, 307, 61, 163, 263, 3191, 1, 293, 1913, 577, 6379, 1489, 2341, 2099, 2557, 1777, 1, 5767, 73, 18617, 1, 1, 10303, 1013, 53, 22613, 3881, 163, 12319, 97, 619, 26681, 4561, 1039, 2053, 9811, 5021, 1, 103, 10739, 16459, 1601, 1907, 35033, 1, 419, 18583, 61, 919, 39317

We see several numbers set to 53 or 61; these would be smooth if we raised our bound a little bit.

This sieving process is where quadratic sieve gets its name from. This drastically decreases the overall work needed to find a sufficient number of smooth numbers, making it practical for very large numbers. This basic implementation could probably handle numbers up to 50-60 digits.

Improvements and optimizations

Quadratic sieve admits a number of "bells and whistles" to dramatically improve its runtime in practice. We mention only a few of the most important ones here.

The simple row reduction method of Gaussian elimination is not able to accomodate the very large smoothness bounds needed to factor large numbers, which often range in the millions, mostly due to space limitations; such matrices, if stored explicitly, would require trillions of bits. However, this method is wasteful, because most of the entries in the matrix are zero (they must be; each number has no more than log₂n prime factors). Instead of using an actual two-dimensional array, we can just keep a list for each column that lists the positions of the 1 bits in that column. We then use a method well-suited to reducing sparse matrices such as the Lanczos algorithm. This still requires a fair amount of space; it's common to use block algorithms that work on small portions of the matrix at one time, storing the rest of the matrix on disk. The matrix step is notoriously difficult to parallelize and for large problems is often done on a single high-performance supercomputer.

The most expensive step by far is the sieving, which can require scanning billions of numbers to locate the needed smooth numbers. A common trick is to only track the approximate logarithm of each number, usually in fixed-point arithmetic. Then, when visiting each number, instead of performing an expensive division we only have to subtract. This introduces a bit of rounding error into the algorithm, but that's okay; by rounding consistently in the correct direction, we can ensure that we don't miss any smooth numbers and only capture a few spurious numbers that we can quickly check and reject. Because the logarithms of small primes are small, and require visiting more numbers than any others, primes like 2 and 3 are often dropped altogether.

Another problem is that a²−n grows fairly quickly; because smaller numbers are more likely to be smooth, we get diminishing returns as we scan higher in the sequence. To get around this, we scan values of not just the sequence a²−n but also a number of similar sequences such as (Ca+b)²−n for suitable constants C, b. This variation is called the multiple polynomial quadratic sieve, since each of these sequences can be seen as the values of polynomial in a.

Finally, although the matrix step does not admit simple parallelization due to many data dependencies, the sieving step is perfectly suited to massive parallelization. Each processor or machine simply takes a portion of the sequence to scan for smooth numbers by itself, returning the small quantity of smooth numbers that it discovers to a central processor. As soon as the central processor has accumulated enough smooth numbers, it asks all the workers to stop. In the multiple polynomial variant, it's common to assign some of the polynomials to each machine.

One peculiar idea for massively parallelizing the sieving step, invented by Adi Shamir, is to use not computers but a specially constructed sieving device based on light emitters and sensors that he calls TWINKLE. The concept is that we have a light for each prime number whose intensity is proportional to the logarithm of that prime. Each light turns on just two times every p cycles, corresponding to the two square roots of n mod p. A sensor senses the combined intensity of all the lights together, and if this is close enough to the logarithm of the current value, that value is a smooth number candidate.

Conclusion

I hope this gives you all some insight into the workings of one of the most powerful factoring algorithms available, and into how unexpected algorithms and mathematics can be applied to familiar problems. Please leave any feedback, whether clarifications, corrections, complaints, or encouragement. Thanks for reading.

Color quantization

MSDNArchive — Wed, 10 May 2006 03:57:00 GMT

If you've ever done work with Web graphics, I'm sure that at some point you reduced an image with thousands of colors, such as a screenshot or photograph, to an image with only 256 colors, for example to store the image in the GIF format. Remarkably, the reduced image is usually very visually similar to the original image, despite the massive loss of information. Most paint programs are capable of performing this transformation, and operating systems such as Windows will perform the transformation on-the-fly when the user is running a video mode with only 256 colors. But how does it work?

In the literature, this problem is called the color quantization problem -- generally, the word "quantization" refers to any process that maps several values to a single representative value. In this case, the values that we're mapping are colors. There are two main steps to color quantization:

Generate a good palette
Dither

Generate a good palette

First we'll talk about generating a good palette. The simplest approach is to use a fixed palette for all images. We can map pixels in the source image to pixels in the reduced image by simply choosing the closest color. What do we mean by "closest"? Well, we can identify each color with a point in three-dimensional space whose coordinates are the red, green, and blue values of the color. Then, we consider two colors to be close if their points are close in this three-dimensional space. This isn't ideal, but it's very simple.

The fixed palette approach has the advantage that all images using the same palette can be rendered on the same screen simultaneously; unfortunately, it severely compromises image quality. Making this even worse is that typical fixed palettes, such as the Windows standard palette and the web colors palette, are chosen as uniform palettes, meaning that the colors are evenly spaced throughout RGB space. This is a terrible way to construct a palette, as humans can distinguish much finer color differences in some areas of RGB space than others.

Below, we show a photograph of a rose taken by Stan Shebs, but modified to remove its blue channel. This means that all of its colors lie in a single plane, the plane blue = 0. That color plane is shown in the graph below. We also show a 16-color optimized palette found by Adobe Photoshop - each blue point is a palette color, and each region shows the section of the color space mapping to that palette color. As you can see, the palette is far from uniform, and there are many colors that the image doesn't even use.

From this picture, we can see that generating a good palette is essentially a point clustering problem. We want to identify clusters in the color space that contain a large number of points close together and assign a single palette entry to those colors.

One of the most popular algorithms for this problem is called the median cut algorithm. It was invented by Paul Heckburt in 1980 in his paper "Color Image Quantization for Frame Buffer Display" and is still in wide use today. The algorithm follows just a few basic steps:

Create a block or box surrounding the points; it should be just large enough to contain them all.
While the number of blocks is less than the desired palette size:
1. Find the block with the longest side out of all blocks.
2. Cut this block into two blocks along its longest side. Cut it in such a way that half of the points inside the block fall into each of the two new blocks (that is, we split it through the median, thus the name).
3. Shrink these two new boxes so that they just barely contain their points.
Average the points in each box to obtain the final set of palette colors.

Despite a number of more sophisticated modern algorithms, this strikingly simple approach is adequate for generating 256 color palettes for most images. You can read about a C++ implemention of the median cut algorithm at my wiki.

Dither

If the palette is large enough and of high enough quality, Reducing by nearest color is usually sufficient. If not, however, distinct artifacts will be visible, such as banding, where what were once smooth gradients become striped regions (see this image; look closely at the cat's neck). One effective technique for dealing with these palette limitations is called dithering, which mixes dots of different colors to simulate a larger number of colors and smooth out gradients. For example, if you mix a grid of red and blue dots, the result will appear purple. Dithering is also commonly used in printing, where for example a grayscale image can be simulated with great accuracy using a higher-resolution dithered black and white image.

The trick to dithering is to effectively simulate the original color in each area of the image while avoiding telltale artifacts formed by poor dithering patterns; see ordered dithering, commonly used in newspapers, as an example of a dithering scheme with clearly visible artifacts. One of the most popular dithering schemes to this day is the Floyd-Steinberg algorithm, a classic dithering algorithm falling into the general class of error diffusion dithering algorithms.

The basic concept behind error diffusion is to begin with nearest color matching, but whenever we can't match a pixel color exactly, the error between the original color and the palletized color is distributed into adjacent pixels not yet visited. In Floyd-Steinberg, the pixel to the right receives 7/16 of the error, the pixel below 5/16, the pixel below and to the left 3/16, and the pixel below and to the right 1/16. When these pixels are visited, this added error affects the palette entry selected for them, leading to an average colour over that portion of the image close to the actual average color of the pixels in that area.

What's great about Floyd-Steinberg is that it generalizes easily to reducing full 24-bit color images to palettes of any size. The only modification we need to make is that we perform the error distribution step for each of the three channels. The palette entry is still chosen by nearest color. You can read about a C implementation of Floyd-Steinberg dithering on my wiki.

This isn't the end of the story of course - there is ongoing research and all kinds of fancy variations on these approaches, many of which can produce much nicer results than these classical algorithms. But at least now the next time you see a GIF or 8-bit PNG file on the web, or use a computer in 256-color mode, you can appreciate the algorithms under the cover that helped make your images look good under these severe limitations.

How does JPEG actually work?

MSDNArchive — Thu, 13 Apr 2006 02:33:00 GMT

JPEG is an image encoding designed to compress photographs and similar images effectively, often 5 to 15 times over a raw bitmap format. It's a lossy format that exploits properties of human vision to eliminate information that is difficult to distinguish. You see JPEGs all the time, and if you've seen a JPEG compressed with a low quality level, you have a good idea of the kinds of visual artifacts that emerge when JPEG compression is pushed to its limits. But how does it really work, and how does this cause JPEG artifacts?

The JPEG standard defines a number of different encoding options, and there are a number of container file formats defined other than the usual JFIF, but I'm going to discuss a single popular JPEG encoding method used by many applications. This encoding combines three different completely independent methods of image compression: downsampling the chrominance channels, quantization of the discrete cosine transformation, and entropy coding.

Method 1: Downsampling the chrominance channels

The most obvious way to make an image file smaller is to make the image smaller - for example, a 100×100 bitmap is 1/4 the size of a 200×200 bitmap. You can expand the image back out to the original size when viewing it by simply upsampling it. This isn't a great compression method, since it causes clearly visible artifacts, usually blocky or blurry effects, especially around detailed, sharp areas.

A different approach is instead of downsampling the entire image, to only downscale certain channels of the image. For example, the eye's red-green sensors are more numerous and sensitive than its blue-yellow sensors, so if we shrink just the blue channel and later re-expand it, the effect is much less noticable. But the RGB color space is not our only choice - images can be encoded in many color systems, and for each of these the channels we choose to downsample will have different visual effects. This experiment by Eric Setton of Stanford University demonstrates the effect on image quality of downsampling various channels in various color systems.

The results are visually obvious: we are sensitive to red and green, but even more sensitive to intensity or luminance, the relative brightness of pixels compared to other nearby pixels. This is expressed by the Y component of the YCbCr system and the V component of the HSV system. Any downsampling of this information is visually disastrous. On the other hand, downsampling the chrominance, expressed by the Cb and Cr components of the YCbCr system, has almost no visual effect (Cb and Cr stand for "chrominance blue" and "chrominance red" - see a sample image broken up into YCbCr channels in the middle of this page). For this reason, typical JPEG files convert images to the YCbCr color space, and then downsample both the Cb and Cr channels by a factor of 2 in both width and height. Further processing of the image proceeds on each of the three channels independently.

This is the method in JPEG encoding largely responsible for the artifact known as "color bleeding", where along sharp edges the color on one side can "bleed" onto the other side. This is because the chrominance channels, which express the color of pixels, have had each block of 4 pixels averaged into a single color, and some of these blocks cross the sharp edge.

Method 2: Quantization of the discrete cosine transformation

Normally, we view an image as simply a grid of pixels, and if we alter a pixel or group of pixels dramatically, the effect will be visually obvious. A different approach to storing images is to come up with a set of "basis" images which can be added together in various combinations to form the image we want. JPEG encodes each 8-by-8 square of pixels using the following basis image set of 64 images:

We assign each of these images a relative weight determining how much influence it has on the image. This weight can be positive or negative (when negative, the basis image is subtracted). Taken together, these 64 weights can describe any 8×8 image at all, and they describe the image precisely. This example at Wikipedia shows an 8×8 image, its original pixel array, and its new 8×8 array of weights after being converted using the discrete cosine transform (which I'll discuss later).

But what's the point of transforming 64 pixel values into 64 weights? For one thing, the resulting data may exhibit more structure, and so be easier to compress. For example, with blocks from typical photographs, the basis images in the upper left above have large weights and the ones in the lower right have small weights.

But to store the weights, which are real numbers, we have to encode them somehow. The simplest encoding is to just convert them to integers between -k and k by scaling and rounding them. This works, but there's a tradeoff in our choice of k: if it's too big, the encoded image will take too much space, but if it's too small, the rounding might visibly distort the image, since it's effectively altering the weights.

A more effective approach is to use a different k for each basis image. This works well because the visual difference caused by variations in the weights of different basis images is different. The eye can detect small variations of the weights of the images in the upper-left, but can overlook much larger variations in the weights of the images in the lower-right, so we can get away with smaller k and larger rounding error for these weights.

In practice, rather than deal with expensive floating-point numbers, we round the initial weights to integers and later divide them by fixed factors called the quantization factors using integer division. During decoding, we multiply them by these factors again. The upper-left images have small quantization factors, so that little information is lost in this process, while the lower-right images have larger factors. This example shows a commonly used matrix of quantization factors, one for each basis image, and the result of quantizing the sample image using them. Observe how much simpler the matrix becomes - it now contains a large number of entries that are small or zero, making it much easier to compress. This example shows how the weight values expand back out during decoding - notice that although some of the weights are now different, the visual difference is almost undetectable.

Quantization is the primary source of JPEG artifacts. Because the images in the lower-right tend to have the largest quantization divisors, JPEG artifacts will tend to resemble combinations of these images. The matrix of quantization factors can be directly controlled by altering the JPEG's "quality level", which scales its values up or down.

Quantization compresses the image in two important ways: one, it limits the effective range of the weights, decreasing the number of bits required to represent them. Two, many of the weights become identical or zero, improving compression in the third step, entropy coding.

But this discussion leaves one important question unanswered: how do we take an arbitrary square of 8×8 pixels and decompose it into a combination of the basis images? As it turns out, we specifically chose the images above in a way to make this easy. These images are formed by cosine waves of various frequencies in the x and y directions, and a transformation called the discrete cosine transform from Fourier transformation theory allows us to rapidly decompose a signal into these parts in only O(n log n) time. Because of the importance of JPEGs, over the years very fast and clever specialised solutions for 8×8 blocks have been created in both software and hardware. This is also the reason we use small, constant-sized blocks instead of processing the image as a whole: the cosine transform step doesn't need to scale up to large bitmaps, so the transformation as a whole can run in linear (O(n)) time.

Method 3: Entropy coding

Unlike the other steps, entropy coding is lossless. If you've studied the simple compression techniques used in ZIP files and BMP files, the techniques of entropy coding should come as no surprise.

We start by turning our 8×8 grid of integers into a linear sequence of integers. The obvious row-major and column-major orders are not ideal here, since we expect there to be many zeros in the lower-right area. Instead, we start with the top-left corner and zig-zag in a diagonal pattern across the matrix, going back and forth until we reach the lower-right corner.

Next, we apply the same run-length encoding algorithm used by GIF files to condense sequences of identical integer values. In areas where the values are small, we can expect many such runs, particularly runs of zeros. We use a special code to say "all the remaining values are zero", which comes in handy.

Finally, we apply Huffman compression to the remaining sequence, which is a standard lossless compression technique that encodes frequently used values using less bits than infrequently used values.

Conclusion

And that's it! For each YCbCr channel, for each 8×8 block, we perform these transformations, and we store the results in a bitstream. This data is wrapped up in a JFIF container giving metadata about the image, and the file goes out ready for web browsers to view, completely unaware of all the data that has been lost. Just don't go sending JPEGs into space and expect aliens to see them the same way.

Finally, no discussion of the JPEG encoding would be complete without noting that JPEG is based on research ideas that are now very old and largely supplanted. The more modern JPEG 2000 encoding based on wavelets can be pushed much farther before visible artifacts appear, and includes a surprisingly efficient lossless mode. Of course, all this improvement is rather moot until common applications actually support it.

Encoding correctness in types

MSDNArchive — Mon, 27 Mar 2006 23:52:00 GMT

This article will discuss effective use of types to catch some common problems at compile time not normally found by the typechecker.

The typechecker is the most important tool in the programmer's arsenal. Types enable you to structure code in a more maintainable way and prevent entire classes of errors from ever occurring at runtime, many of which can be insidiously difficult to debug. If you've written anything extensive in assembly language, you can understand the frustration of dealing with misinterpretations of data leading to corruption of program state. Researchers like to say "well-typed programs never go wrong".

But what does "go wrong" mean? One answer is: "more than you think". By making effective use of the type system to encode invariants about your program, you can use the compiler to catch more errors at compile time that you currently catch at run time.

A simple example is SQL injection, where insertion of user-supplied data into a SQL statement string can enable attacks where the user modifies the nature of the statement with carefully chosen data. At first the solution seems simple: double your apostrophes to escape them, and then the value cannot be misinterpreted. This can be trickier than sounds, though, since in a complex program you might not be able to identify which values come from an untrusted source, directly or indirectly, which are safely produced from trusted sources, and which are untrusted data that has already been properly escaped. These three types of data travel from method to method, following a winding path that obscures their source. The compiler is of no help, since it sees them all as just strings.

The solution, called "tainting" in languages like Perl and Ruby, is to make them not "just strings", but to create two types of data: tainted data, which is unsanitized data derived from an untrusted source, and untainted data, which is either derived from a trusted source or already sanitized. For our example, suppose we're performing our queries with SqlDataReader. We create a wrapper class for it which has an identical interface except that wherever SqlDataReader takes a string, our wrapper takes a SafeSqlString. This new type is simply a wrapper for string that can only be generated in a few strictly limited ways:

static SafeSqlString AssumeSafe(string): Used for trusted data, just wraps it
static SafeSqlString Escape(string): Escapes apostrophes; used for untrusted data.
SafeSqlString operator+(SafeSqlString): Concatenates two SafeSqlStrings
static SafeSqlString Format(SafeSqlString format, params object[] args): Takes a safe format string and a list of arguments of safe types such as SafeSqlString, int, and double, and produces a new SafeSqlString.

What does this buy us? Now, as long as we use our new SqlDataReader wrapper for all database access, we are essentially immune to SQL injection: the only way it could sneak in is through AssumeSafe(), which is obviously dangerous and makes the code easy to review. As an added bonus, double escaping can never happen: the data only gets escaped once as it changes from string to SafeSqlString, and Escape() cannot be applied to SafeSqlString.

Because these types flow through the program, we can specify that a particular method always returns a SafeSqlString, or assumes that its arguments are SafeSqlStrings, or assumes that some private variable is always safe, and enforce this statically instead of by convention. This is important, because it allows us to use AssumeSafe() at a point where it's obvious that the string is safe, as opposed to five calls up the call stack where we really have no idea who else might be calling us.

What else can we encode in types? All sorts of things! Here are some more examples:

We could have a type PositiveInt for positive integers, which will throw a runtime error if constructed with a nonpositive integer. Now, if a method receives a PositiveInt, it can safely pass it on to another method expecting a positive integer without any additional runtime checks. This trick works with any object where you want to enforce a complex runtime property on an object but only want to check it once (or only when it might change). In unmanaged code, it's particularly useful for checking buffer/string lengths.
Some binary tree algorithms require that every node has zero or two children, never one. We can represent this using a subclass hierarchy, where the abstract base class is TreeNode and the subclasses are Leaf and InternalNode. Leaf contains no data at all, while InternalNode contains a piece of data and two (non-null) references to TreeNodes.
There are many other properties of data structures that can be encoded in types. For example, in a red-black tree, red nodes cannot have red children. We can enforce this statically by creating a subclass hierarchy where BlackNode objects contain references to the abstract base class Node, while RedNode objects contain only references to BlackNode objects. Similarly, the root must be black, so we just use a BlackNode reference for the root.
If you're writing a sorting algorithm, it might not be clear that the array you get back contains a rearrangement of the same elements as the array you passed in. There's no quick way to check this, even at runtime. On the other hand, if instead of an array you pass in an object whose only valid non-read-only operation is to exchange two elements, then this becomes statically obvious.

In short, wherever you find yourself wanting to keep track of a static property as it flows through a program, types are a great way of doing so - use them to catch errors sooner.

In programming languages with more sophisticated type systems, you can encode even more information in types. Taking this to the extreme are dependent types, a relatively new system of types based on values that are extremely difficult to typecheck but so powerful that they can express many correctness properties of real programs. A number of research languages such as Dependent ML and Epigram already implement dependent types. With the power of systems like these, maybe in a few years the primitive C-based type systems of C# and C++ will look a lot like assembly language does now. But before we can take advantage of new technology, let's start making better use of what we have - encode your program's properties using types.

LiteratePrograms wiki

MSDNArchive — Sat, 04 Mar 2006 04:54:00 GMT

Hi everybody. This is a short post, but I just wanted to tell you all about a new wiki I created called LiteratePrograms, based on an extension of the same MediaWiki software used by Wikipedia.

Some of you read my earlier post about literate programming and some of the advantages it has for easy-to-read exhibition of code. On LiteratePrograms, every article is a literate program: it is simultaneously a document and a piece of code which can be downloaded, compiled, and run. The "download code" tab at the top of every article runs the noweb literate programming tool which extracts code segments from the discussion, puts them back in a form that can be compiled, and zips them up in a package for download. Additionally, the wiki markup, editable by anyone, allows rich documentation including mathematical formulas and diagrams. Finally, all content is released under the liberal MIT/X11 non-copyleft license, so anybody can use it, even for commercial purposes. You can see a sample of LP in action at Insertion sort (C, simple).

Recently featured at Lambda the Ultimate (thanks Allan!), LP has attracted a lot of attention and has a couple dozen articles and users now. The users have some pretty diverse interests in languages and subject areas that I would never have anticipated (see Pi with Machin's formula (Python), Hello World (occam)).

I invite your suggestions for and participation in and the site - in particular, we're really lacking in Windows and .NET programs and programmers at the moment, and I know you guys can help. See How to write an article and take a look at the existing articles to learn about the markup syntax. Please ask me any questions you have. Thanks everyone.

Dependency tracking in builds

MSDNArchive — Wed, 15 Feb 2006 06:20:00 GMT

Happy Valentine's Day, everyone! Today I'm going to talk about the important problem of successfully tracking dependencies in builds in an automatic manner.

In short, there are two kinds of build systems: the ones where all dependencies are tracked successfully and incremental rebuilding is available at all times, and the ones where you have to rebuild everything whenever you change anything remotely significant. There is nothing inbetween - as soon as a build system fails to track dependencies in one circumstance, developers will become paranoid and insist on costly full rebuilds in the future. After all, the consequences of a partial rebuild are insidiously difficult to diagnose, resulting in broken builds, symbol mismatches, exceptions while loading assemblies, mixtures of interacting new and old code, and other nonsense that we don't want to waste our time debugging. This means that a dependency tracking system has to be very close to perfectly reliable, which can make it very difficult to design and implement successfully.

Although the problem is a very old one, with tools like GNU make dating back decades, I believe in most project environments dependency tracking is not solved to the level of reliability that is needed for this kind of trust. Behind all of these failed dependency-tracking systems is one cardinal sin:

Dependencies cannot be maintained by hand.

I'm not merely saying maintaining dependencies by hand is difficult or even error-prone, but that it's damn near impossible. In C/C++ in particular, adding a #include statement to a header file implicitly adds a dependency to every source file including that header file directly or indirectly. Even worse is removing a #include, in which case some of the files depending on it might no longer require it, but then some might, depending on what they directly or indirectly include. Waste all the time you want on this - you won't get it right. Even in C#, Java, or Visual Basic, which thankfully have no concept of header and source files and load an entire project for compilation at once, interproject dependencies can become just as complex in large solutions.

Even worse are systems which require you to explicitly specify not only the dependencies but the build order. Often hierarchically structured solutions of many projects will use pretty good dependency tracking within each project, but revert to this "build these things in this order" model for the whole tree. The result is that related modules in different parts of the tree cannot be built together incrementally, and highly parallellized building becomes impossible. Scalability demands that we create some level of abstraction for solution-level builds, but that doesn't mean we have to sacrifice incremental building altogether.

The alternative to maintaining something by hand is of course maintaining it automatically. The code is already quite explicit about what depends on what - whether in the form of #include statements, import statements, assembly references, or other means. The fundamental principle of avoiding redundancy tells us that we should be generating dependency information from these sources. For example, the following buried treasure from the GNU make manual details how to use the -M flag of gcc together with some (rather esoteric) script to automatically track header dependencies in C:

The practice we recommend for automatic prerequisite generation is to have one makefile corresponding to each source file. For each source file `name.c' there is a makefile `name.d' which lists what files the object file `name.o' depends on. That way only the source files that have changed need to be rescanned to produce the new prerequisites.

Here is the pattern rule to generate a file of prerequisites (i.e., a makefile) called `name.d' from a C source file called `name.c':
%.d: %.c
        @set -e; rm -f $@; \
         $(CC) -M $(CPPFLAGS) $< > $@.$$$$; \
         sed 's,$$*$\.o[ :]*,\1.o $@ : ,g' < $@.$$$$ > $@; \
         rm -f $@.$$$$
http://www.gnu.org/software/make/manual/html_mono/make.html#SEC51

The concept here is that the dependencies for each source file are stored in a separate file which itself depends on that source file. If the source file is modified, the ".d" file containing that source file's dependencies is rebuilt before it is loaded by make. In this way the dependencies can never be wrong or out of date for even a moment. This concept could be generalized to a pluggable framework that enables you to determine build-time dependencies of arbitrary files using arbitrary code.

We can play similar tricks to dynamically track the dependencies of projects in a solution. Although this is one level of abstraction higher, it's really just the same thing with source files replaced by projects. We analyze the projects to determine their dependencies at build-time, set up a build order, and then locally build each project using its own local dependency structure. This enables automatic dependency tracking to scale effectively. With a solution like this, you can simply sync up the latest changes from the repository and do an incremental build from root in minimum time.

Dependency tracking isn't just useful for builds either - it can often be used to evaluate the impact of a change. For example, if I change a header file, a good dependency tracking system should be able to name for me all the affected source files in the project, as well as any other project that might be affected. In the end I might decide to go with a less destabilizing change instead. They can also be incorporated into version control systems, where they can help to evaluate whether concurrent changes made by other people will break what you're working on right now. This is largely unexplored territory.

At Microsoft, a new build system called MSBuild has recently been developed and released with Visual Studio 2005. Currently used for almost all of the Visual Studio build, it will hopefully in time replace older ad hoc systems based on batch files across the company. One of its advantages is its incremental building support, described here:

http://msdn2.microsoft.com/en-us/library/ms171483.aspx

Although a sophisticated system, even MSBuild's support for incremental building is relatively limited, focusing on dependencies within projects between files that have direct input-to-output mappings.

Build scheduling algorithms

Let's talk a little bit about the technology behind build scheduling - a simple problem, but one with a number of interesting generalizations. To simplify the discussion, we'll refer to both inputs and targets/outputs as "nodes". Thus, source files, header files, binary outputs, object files, resource files, and so on are all nodes. We can describe this as a directed graph where each vertex is a node and each edge from A to B says that "B depends on A". In simple terms, this is just a drawing where we have a circle for each node and an arrow pointing from each node to nodes that depend on it. This is called the dependency graph, and is the fundamental data structure used in build scheduling (as well as other types of scheduling).

The simplest scheduling problem is sequential scheduling, where we have a one-processor machine and we need to figure out what order to build things in. In this case, we can use topological sorting. This just means that first we start by building all the things with no dependencies, and we proceed to build the things whose only dependencies are things we already built. There are efficient algorithms for determining this order in advance.

In large solutions, however (like say, Windows), building everything on one machine could take weeks or months. In this case, we'd like to use many machines, each with many processors, to build different parts of the product simultaneously. However, they still have to respect the build-time dependencies between components. This is where things get interesting, and we run into many of the same problems as scheduling instructions on microprocessors, and even run into some NP-complete problems. A couple common heuristics, the Decreasing-Time Algorithm and the Critical-Path Algorithm, are described here:

http://www.ctl.ua.edu/math103/scheduling/scheduling_algorithms.htm

Both of these require a time estimate that each step will take, which is not too difficult to generate automatically from metrics like lines of code. The concept behind the Critical-Path Algorithm is that if a long chain of dependencies is waiting on a particular node to be built, that's a good reason to build that node first, even if it doesn't take long to build itself. Even these heuristics are based on NP-complete problems such as the longest path problem, but in practice they can often be solved efficiently.

Note that one can't just throw arbitrary hardware at a scheduling problem; there are limits beyond which any additional hardware would simply sit idle, and there are even smaller limits for which very little benefit is derived. Part of the challenge of parallel scheduling is finding just the right number of systems to maximize cost-benefit.

There are even fancier variants of this where you have machines that build at different speeds, machines that can build some types of files but not others, and so on. I won't get into these here, but this discussion should make one thing clear: a developer doesn't want to waste precious development time figuring all this stuff out. Build order should always be determined by build technology that can be relied upon to be correct, to do efficient incremental building, and to scale up to large builds and large build clusters well, all automatically.

So why do we keep using hacked-up internal build systems made of Perl scripts and anarchy? I don't know, but if you decide to make something better, please come sell it to my group!

The nature of computing and infeasible machines

MSDNArchive — Tue, 07 Feb 2006 05:44:00 GMT

If asked what they do, some developers might say "I write code" or "I program computers". But more fundamentally, what is computing really about? After all, Euclid's algorithm for computing the greatest common denominator of two numbers in polynomial time was designed - and executed by hand - thousands of years ago, long before even early mechanical computing devices were available (besides the abacus). For centuries elaborate bureaucracies have set up processing units and communication channels to spread and delegate work, not unlike modern server farms. Even today, complexity theorists regularly study mythical machines such as nondeterministic machines, alternating machines, nonuniform circuits, oracle machines, and hypercomputers, none of which are practical to construct or even to simulate in general. So why design algorithms for machines we can't build?

In reality, computing has nothing to do with computers. Computing is fundamentally about solving problems using strictly defined processes under resource constraints. The theory of computing, by design, has little regard for what device you use to carry those processes out.

For example, suppose ten monks in medieval Europe are copying a sacred book of 100 pages, numbered 0 through 99, and just as they finish one of them accidentally knocks the pages all over the floor. How can they rapidly put the pages back in order? One simple process is this: each monk picks up 10 pages and places them in one of ten piles depending on the tens digit of the page number. Then, each monk takes one of the piles and sorts it by himself. The piles are put together and they are done - this could all be done in just a couple minutes, much more quickly than a single monk would take to sort the whole book. In computing terms, this is a parallelized version of radix sort, switching to selection sort or insertion sort for small sublists.

Indeed, good algorithms are often more important than the device you use; for example, a computer sorting a list of a billion integers using bubble sort would probably take more time than a room full of people using the parallelized radix sort described above. On the other hand, some devices can solve problems that others cannot; for example, the very weak deterministic finite automaton (DFA) can't even tell whether its input string contains the same number of zeros and ones, which is a trivial task for a Turing machine, an abstract machine modelling practical computers. This tells us something important: DFAs aren't good enough. In other words, a computer that at its best can only simulate a DFA is fundamentally incapable of solving certain problems. Similarly, theoretical hypercomputers can solve problems that Turing machines cannot.

Unfortunately, once machines become powerful enough, they quickly all become equivalent with respect to the problems they can solve. For example, deterministic, multitape, random access, probabilistic, nondeterministic, and quantum Turing machines all solve the same set of problems. Why go to the bother of defining all these different machines then? There are two reasons:

Certain devices make solutions to certain problems much easier to describe. For example, problems that are easy to solve by exhaustive search are also easy to express using algorithms that run on nondeterministic machines. Problems like primality testing and undirected graph connectivity are much easier to solve with known algorithms when we have the ability to generate random numbers.
Some devices can solve more problems than others using the same amount of resources, or solve the same problems using less resources. As a simple example, a machine with random access can search an ordered list in logarithmic time using binary search, whereas an ordinary Turing machine with only sequential access would require linear time. Shor's algorithm can factor integers on quantum computers in polynomial time, something we don't know how to do on classical machines.

We know that the "features" offered by new theoretical devices are useful in describing new algorithms for solving problems, since we designed them for that purpose, but what we'd really like to know is whether they really let us solve problems more efficiently than we possibly could without them. This question, on the "power" of devices, is the major stumbling block of modern complexity theory. Even the great unsolved question of whether P=NP can be phrased as a question of power: can nondeterminism allow us to solve problems in polynomial time that we could not otherwise? Another question, whether P=RP, asks whether we can solve more problems in poly time with the ability to generate random numbers than without it. Although also unsolved, many researchers believe that these classes are equal and randomness does not add any real new power. The unsolved question of whether L=P asks whether we can solve more problems in polynomial space and polynomial time than in logarithmic space and polynomial time. Occasionally, we show that two dissimilar devices are equally powerful, as in the famous results that PSPACE=NPSPACE, IP=PSPACE, and NL=coNL, but rarely ever have we unconditionally shown that one machine is truly more powerful than another (notable exceptions are the trivial L ≠ PSPACE and P ≠ EXPTIME; see time hierarchy theorem, space hierarchy theorem).

To come back to my original point, why do we design algorithms for machines we can't build? There are many reasons, but some of the most important are:

The algorithm may be more natural to describe on an abstract machine which we can then simulate on a real machine. By focusing our efforts on this simulation program, we can acquire the ability to easily solve any problem that the machine being simulated can solve.
Understanding what problems a machine can solve helps us learn about the machine and how powerful it is. This helps us determine whether we should invest money and effort in software and hardware simulating those features or whether we can adapt existing machines to solve the same problems just as efficiently.
Like the proof that the halting problem cannot be solved, they help us understand what problems are truly intractable, so that we can focus our effort on problems that we're more likely to crack successfully.

More philosophically, a popular saying among mathematicians is that if you were to wake up tomorrow and the universe were gone, you could still practice mathematics - abstract concepts are not dependent on any property of our physical universe, even if the specific representations and concepts that we use are. The same applies to computing, which at its core is a mathematical discipline: even when modern SDKs and classical architectures have become obsolete, the same abstract problem solving we now perform in C++, C#, or Java we will perform in yet another new programming language. And if you're ever stuck in purgatory for eternity, you can just quicksort a list of a billion numbers in your head.

Literate programming

MSDNArchive — Tue, 24 Jan 2006 00:53:00 GMT

Literate programming, invented in 1981 by the same Donald Knuth who wrote The Art of Computer Programming and the document language TeX, is a technique in which a program is written as a human-oriented document interspersing discussion and code. The code segments are arranged not according to execution order or the logical structure of the code, but in whatever order the author believes offers the best presentation to the reader. A tool can be used to automatically extract the code segments, put them back in the order expected by the compiler, and export them as a source file. Although it never really caught on for ordinary coding, it has the important ability to make code easier to understand for a reader, easing maintenance, transfer of ownership, and helping to prevent future bugs.

The idea is perhaps easiest to understand if we look at it not from the perspective of a coder, but from the perspective of a researcher writing a book. They write a section about an algorithm, showing first the core algorithm, then talking about issues, improvements, and variations. In each section they show a little more code. They might show a test driver to demonstrate how it might be invoked. Taken to the logical extreme, these code samples could, taken together, show every part of the program - all someone would need to do to run it is to copy the samples into a source file in the right order and build it. If we automate this extraction process, we now have a section of a book that simultaneously is a real program that can be run.

Some people find the idea of structuring a program in this seemingly random "presentation order" strange, but really we do it all the time. Suppose you arrive on a new team and take control of a module from someone else. At some point, you sit down with the previous owner, who proceeds to describe how the module works. If they do this effectively, they don't take an alphabetical file listing and go through each file in order - instead, they first focus on core parts of core functionality and expand outward to more peripheral functionality. You might ask a question about a piece of code you see, and they say "Don't worry about that - I'll get to that later when I talk about X." By controlling the presentation order, they enhance your comprehension by keeping code out of the way until you have the proper background to understand its purpose.

But what if this previous owner is long gone? They could have all the comments in the world, but you still won't have anything like the discussion above. Literate programming enables them to write a tutorial document that takes you through the code in the same friendly way as a sit-down with the owner, and yet never invokes redundant effort or becomes out-of-date, because it is the code as well.

An example

But there's no better way to understand literate programming than to see some. So, right here, I'm going to discuss an implementation in C# of Bloom filters, the topic of my very first post (read this first if you haven't already). This implementation is based on a paper by a former of professor of mine, Panagiotis Manolios, and one of his students, Peter Dillinger, called Bloom Filters in Probabilistic Verification. This paper solves a frustrating implementation issue: Bloom filters require a potentially unlimited number of independent hash functions, and there's no convenient way of coming up with these hash functions on the fly. Their concept is to use only two independent hash functions to generate all of the hash values required. The insertion method below, based on Algorithm 2 in their paper, demonstrates this:

<<BloomFilter Insert method>>=
    public void Insert(TKey obj)
    {
        int h1 = obj.GetHashCode();
        int h2 = obj.GetHashCode2();
        unchecked
        {
            h1 = (int)(((uint)h1) % table.Count);
            h2 = (int)(((uint)h2) % table.Count);
        }
        <<perform insertion>>
    }

The syntax <<perform insertion>>, used by the noweb tool, indicates a chunk of code that I haven't talked about yet. The label at the top is the name of this chunk of code. We use unchecked arithmetic for the casts because the hash code methods can return negative values, which we wish to convert to large positive values via the cast to uint. The narrowing cast back to int is safe because table.Count is in an int.

Above, TKey is a generic type parameter on the class specifying the type of values stored in the Bloom filter. Keys must implement the IDoubleHashable interface, which adds a second hash function GetHashCode2() independent from the standard GetHashCode() supplied by Object:

<<IDoubleHashable definition>>=
    public interface IDoubleHashable
    {
        int GetHashCode2();
    }

<<BloomFilter class definition>>=
public class BloomFilter<TKey> where TKey : IDoubleHashable
{
    <<BloomFilter implementation>>
}

We store the bits of the BloomFilter in a standard BitArray collection, which compactly stores the bits of the table and handles all the bit arithmetic for us:

<<BloomFilter member variables>>+=
        private BitArray table;

This is in the System.Collections namespace, so we need a using statement:

<<BloomFilter using statements>>+=
using System.Collections;

The += noweb syntax above is used to add a piece of code to a code section. This is useful for sections we wish to describe sequentially a part at a time, or for elements like the above where order is unimportant.

Now, we're prepared to describe the insertion. We set a number of roughly evenly spaced bits, where the first hash establishes the initial index, the second hash establishes the spacing, and the spacing is slightly increased as we go along. All arithmetic is done modulo the table size:

<<Perform insertion>>=
        for(int i=0; i<numHashes; i++) 
        {
            table[h1] = true;
            unchecked
            {
                h1 = (int)((uint)(h1 + h2) % table.Count);
                h2 = (int)((uint)(h2 + i) % table.Count);
            }
        }

The number of bits we set is determined by the fixed parameter numHashes. We'll describe how it's chosen a bit later on:

<<BloomFilter member variables>>+=
        private int numHashes;

The lookup procedure is analogous, checking all the same bits and returning true only if they're all set:

<<BloomFilter Contains method>>=
        public bool Contains(TKey obj)
        {
            int h1 = obj.GetHashCode();
            int h2 = obj.GetHashCode2();
            unchecked
            {
                h1 = (int)(((uint)h1) % table.Count);
                h2 = (int)(((uint)h2) % table.Count);
            }
            for (int i = 0; i < numHashes; i++) 
            {
                if (table[h1] == false)
                {
                    return false;
                }
                unchecked
                {
                    h1 = (int)((uint)(h1 + h2) % table.Count);
                    h2 = (int)((uint)(h2 + i) % table.Count);
                }
            }
            return true;
        }

The probability that this method will return true for a random value not in the set is given by:

p = (1 − e^−kn/m)^k

where the variables are defined as:

k: Called numHashes in the code, the number of table elements set by an insertion
m: Called table.Count in the code, the number of elements in the table
n: The number of distinct objects inserted so far

To initialize the object, we use the following method, which needs to know k and m:

<<BloomFilter initialization>>+=
        private void Initialize(int tableSize, int numHashes)
        {
            this.numHashes = numHashes;
            table = new BitArray(tableSize, false);
        }

However, a user is unlikely to know the best values for these (which is why this isn't a constructor). Instead, we address two common use cases:

In this case, the user knows m, the desired table size in bits, and n, an estimate of how many elements they will insert. We choose k in a way that minimizes the false positive probability p, using the value (ln 2)m/n rounded to the nearest integer:
```
<<BloomFilter initialization>>+=
        public BloomFilter(int tableSize, int numElements)
        {
            Initialize(tableSize, (int)Math.Round(ln2*numElements/tableSize));
        }
```
In this case, the user knows n, an estimate of how many elements they will insert, and p, the desired false positive probability. Bloom filter theory states that if k is chosen optimally (using the formula in part 1 above), we will have p = α^m/n, where the constant α = 1/2^{ln 2}. This means we must have m = (log_α p)n. I call α optimalBaseRate below:
```
<<BloomFilter initialization>>+=
    public BloomFilter(int numElements, double falsePositiveProb)
    {
        int tableSize = (int)Math.Ceiling(numElements * Math.Log(falsePositiveProb, optimalRateBase));
        Initialize(tableSize, (int)Math.Round(ln2 * tableSize / numElements));
    }
```

The constants used above are defined at the class level (we cannot use const because the Math calls are considered non-constant by the compiler):

<<BloomFilter constants>>+=
    private readonly double ln2 = Math.Log(2.0);
    private readonly double optimalRateBase = 1.0/Math.Pow(2, Math.Log(2.0));

We need the System namespace for Math:

<<BloomFilter using statements>>+=
using System;

Finally, we put it all together to complete the implementation of BloomFilter:

<<BloomFilter implementation>>=
<<BloomFilter constants>>
<<BloomFilter initialization>>
<<BloomFilter member variables>>
<<BloomFilter Insert method>>
<<BloomFilter Contains method>>

We structure the program as two source files, with the interface in its own file, and both entities in the BloomFilters namespace:

<<BloomFilter.cs>>=
<<BloomFilter using statements>>
namespace BloomFilters {
<<BloomFilter class definition>>
}

<<IDoubleHashable.cs>>=
namespace BloomFilters {
<<IDoubleHashable definition>>
}

Our implementation is now complete. My full version of this code includes double-hashable classes for integers and strings and a test driver, which I omit here for space. With a bit of added syntax, I could run this entire entry through the literate programming tool noweb, and it would produce both an HTML document much like the one you're reading now, along with two C# source files which I can then build into an assembly and use.

More information

For a longer, fancier example of literate programming, take a look at my Select implementation, which uses the literate programming tool noweb to write an implementation of the classical worst-case linear time selection algorithm using C++ and the document language LaTeX. Here we gain the additional advantage of a full-fledged document language, complete with the ability to typeset mathematical expressions and diagrams, which is impossible in ordinary source code.

Sometimes schemes such as Perl's pod, Javadocs, or C#'s XML comments are touted as being literate programming, but this is far from the case. These are simply isolated documentation of individual logical program units. While they are valuable and have little overhead, they do not present the program; they only document parts of it in an isolated way.

So when should you use literate programming? Certainly there is an overhead involved, but in any case where the ability to maintain a portion of a program is a critical factor, or where there is a high risk of an owner disappearing, or where code will be exposed to many external partners or customers, literate programming is a valuable tool for keeping the code accessible and easy to understand. Try rewriting a single confusing or troublesome portion of your program with it, update your build system to run the literate programming tools prior to compilation, and see if it makes a difference with your application.

Efficient selection and partial sorting based on quicksort

MSDNArchive — Wed, 18 Jan 2006 23:26:00 GMT

Most people who have studied algorithms remember quicksort, the ubiquitous sorting algorithm available in the standard library of nearly every programming language implementation in existence, including C, C++, Java, and the .NET Framework. Quicksort is a compelling example of how algorithms with poor worst-case behavior but good average-case or expected behavior can be highly practical. Much less well-known, however, are the simple variations on quicksort that can improve stack space usage and solve other important problems efficiently.

Let's briefly review quicksort. Quicksort is fundamentally based on the efficient partition operation, a linear-time operation which divides a list into two parts: the elements less than a value, and the elements greater than or equal to a value. The value is called the pivot value, and is usually chosen to be one of the values in the list. In efficient implementations of quicksort, partitioning is usually implemented as an in-place operation that just shuffles elements of the array around as it scans over it. Quicksort operates by simply choosing a pivot value, running partition on the array, and then recursively running quicksort on the two partitions:

 // Sort the elements of a between indexes left and right inclusive
 void quicksort(a, left, right) {
     if (right > left) {
         Choose a pivot element a[pivotIdx] with left <= pivotIdx <= right
         // Here, partition() returns the new index of the pivot element
         newPivotIdx = partition(a, left, right, pivotIdx);
         quicksort(a, left, newPivotIdx-1);
         quicksort(a, newPivotIdx+1, right);
     }
 }

Much research has gone into the problem of choosing a good pivot element, but simply using a typical pseudorandom number generator to select an array index is usually sufficient to get efficient performance on average. A common and effective optimization is to switch to insertion sort for small subarrays, for which insertion sort is considerably faster than quicksort (and pretty much every sorting algorithm). A useful variation on quicksort is the partially tail-recursive variant, where we only make a recursive call for the smaller of the two sublists, and use a loop to deal with the other one:

 // Sort the elements of a between indexes left and right inclusive
 void quicksort(a, left, right) { 
     while (right > left) {
         Choose a pivot element a[pivotIdx] with left <= pivotIdx <= right
         newPivotIdx = partition(a, left, right, pivotIdx);
         if (newPivotIdx - left > right - newPivotIdx) {
             quicksort(a, newPivotIdx+1, right);
             right = newPivotIdx - 1;
         } else {
             quicksort(a, left, newPivotIdx-1);
             left = newPivotIdx + 1;
         }
     }
 }

This variant has worst-case logarithmic space usage, and in practice makes about half as many calls, making it more competitive with in-place sorting algorithms like heapsort. It's especially useful on memory-limited systems, but since stack space often comes at a premium, it's also useful on desktops.

Now let's turn to how we can solve different problems using variants on this idea. Related to sorting is the selection problem, where we are asked to find the kth smallest element in a list. A common special case of this is finding the median. Obviously, we can just sort the list and index, but this isn't terribly efficient - this is only a good idea if we're going to be selecting many different elements and the list is static. Instead, notice this important property of the partition function: the pivot element is always placed in its final sorted position, while all elements which come before it in sorted order are before it, and all those that come after it in sorted order are after it. Consequently, if we know the sorted position of the element we're looking for, we already know which partition it's in - just check k against the index of the pivot after the partition operation:

// Get the value of the kth smallest element in the list
object quickselect(a, k) { 
     int left=0, right = a.Length - 1;
     while (right > left) {
         Choose a pivot element a[pivotIdx] with left <= pivotIdx <= right
         newPivotIdx = partition(a, left, right, pivotIdx);
         if (k < newPivotIdx) {
             right = newPivotIdx - 1;
         } else {
             left = newPivotIdx + 1;
         }
     }
     return a[k];
 }

This is nearly the same as the second quicksort example above, but we no longer need to make any recursive calls - this algorithm is in-place. Moreover, because the sizes of the lists we visit shrink geometrically on average, the total expected time is only O(n), or linear time. The algorithm is also highly efficient in practice, and is used in all implementations of C++'s nth_element call. There is a variation on this scheme with worst-case linear time that I won't discuss here, because it's not as efficient in practice, but it was one of the most important early theoretical results in algorithms.

A related problem is finding the k smallest or largest values in a list, or the "top k". This is important, for example, in cases where you want to display just the top results from a search based on some numeric metric of relevancy (or if you want to start your own pop music countdown). The simplest algorithm, running the first k iterations of selection sort, requires O(kn) time, which is quite expensive. Instead, we can exploit the properties of partition to avoid many of the recursive calls made by a complete quicksort:

 // Rearrange the elements so that the k smallest elements appear at the beginning
 void quicksortSmallestK(a, left, right, k) { 
     while (right > left) {
         Choose a pivot element a[pivotIdx] with left <= pivotIdx <= right
         newPivotIdx = partition(a, left, right, pivotIdx);
         if (k <= newPivotIdx) {
             right = newPivotIdx - 1;
         } else if (newPivotIdx - left > right - newPivotIdx) {
             quicksort(a, newPivotIdx+1, right);
             right = newPivotIdx - 1;
         } else {
             quicksort(a, left, newPivotIdx-1);
             left = newPivotIdx + 1;
         }
     }
 }

This algorithm uses O(n + klog k) expected time and O(log k) space, and again is quite efficient in practice. It can also be generalized to compute an arbitrary subarray of the sorted array (the jth through kth smallest values in the list in order).

When it comes to solving ordering problems, partitioning and quicksort-based algorithms are really your swiss army knife. Keep them in mind next time you encounter a problem like this.

Integer division by constants

MSDNArchive — Tue, 13 Dec 2005 04:07:00 GMT

Today I'm going to discuss a few clever algorithms for dividing values by constant values known at compile-time, without the use of the division instruction.

In compiler optimization theory, the basic problem we attempt to solve is to create algorithms which take a specific program as input and produce a better program, for example a program which runs faster, uses less code space, or even consumes less power. An efficient standard library is an important step in creating efficient programs, but it does not provide the advantage of optimizations, which can dynamically adjust their behaviour to suit each input program. To demonstrate, today I'm going to discuss a very narrow but suprisingly interesting compiler optimization: efficiently dividing an (unsigned) integer by a fixed constant, such as 37.

In the bad old days, the division instruction on x86 machines was extremely slow, and the use of clever bit tricks was very common to speed it up, especially in high performance areas like graphics and simulation. Today we're more concerned with embedded systems, where, to save on chip cost, division is often implemented as a software library routine. Such a routine is necessary if neither argument is known until runtime, but it is extremely slow, easily taking hundreds of clock cycles. In many common applications, we know the divisor at compile-time, and we can exploit this knowledge to craft a specialized division algorithm that drastically outperforms the generic one.

Most programmers are aware that a value can be divided by 2ⁿ, truncating the remainder, by merely shifting it right by n bits. This works because the positional binary numeral system multiplies the kth digit from the right by 2^k. This simple optimization is present in nearly every compiler; for this reason, overzealous programmers who convert their divisions by 4 into bit-shifts will usually discover that there is no improvement (at least if they measure the improvement, as they always should).

But what if we want to divide by a different number, say, 10? This would be quite handy to do efficiently, since division by 10 is often used in converting from binary to decimal or binary-coded decimal. Since 10 = 5×2, we really just have to shift right by 1 and then divide by 5. Dividing by 5 is the tricky part. There are several approaches; the following ideas are taken from the pre-released Chapter 10: Integer Division by Constants of the next edition of Hacker's Delight, an excellent book full of neat bit tricks like these.

First, supposing a multiplication instruction is available, let's think about how we can divide using multiplication. First notice that dividing x by y is the same as multiplying x by 1/y, and we know y at compile time, so we can find 1/y then. We use a fixed-point representation to perform the fractional arithmetic using integers. For example, in the decimal system, one can divide by 5 by simply multiplying by 2 and then moving the decimal point one place to the left (x/5 = 2x/10). In binary, however, 1/5 does not have a terminating decimal; instead we have:

1/5 = 0.001100110011...₂

Now, suppose we have an instruction that multiplies two 32-bit values and produces a 64-bit result. Just as x/5 = 2x/10, we also see that:

x/5 = (858993459.2)x / 2³².

We can approximate (858993459.2)x by (858993459.25)x, which is approximated in turn by 858993459*x + (x >> 2). We just do one multiply, one shift, and one 64-bit add, then take the high word (the high 32 bits). This is just a heuristic argument, but sure enough, if we exhaustively try this formula on every 32-bit value, it never fails. On many machines like the x86 this can be done with as little as 4 instructions.

But what about the machines not fortunate enough to have a 64-bit multiply result? Or what if you're writing in a high-level language like C where multiplication always has the same size operands and result? In this case the task is more difficult. The key is to notice that simulating an approximate multiply by 1/5 isn't so hard, because its bit pattern is very regular: 1/5 is about 0.00110011001100110011001100110011. Now, notice that:

0.00110011001100110011001100110011 = 0.0011001100110011 + (0.0011001100110011 >> 16)
0.0011001100110011 = 0.00110011 + (0.00110011 >> 8)
0.00110011 = 0.0011 + (0.0011 >> 4)
0.0011 = 1 >> 3 + 1 >> 4

Applying this logic in reverse, we devise the following algorithm for approximating x/5:

    q = (x >> 3) + (x >> 4);
    q = q + (q >> 4);
    q = q + (q >> 8);
    q = q + (q >> 16);

Unfortunately, this algorithm isn't correct off the bat like our last one, because we're discarding bits and accumulating rounding error all over the place. However, using multiplication, we can get an idea of how far off we are:

    r = x - 5*q;

By adding a carefully devined error adjustment based on this value, to be precise 11r/32, we come out with a correct algorithm that can be verified by exhaustive testing. See the PDF linked above if you're interested in the precise details of how we come up with this adjustment. Our complete algorithm reads:

    q = (x >> 3) + (x >> 4);
    q = q + (q >> 4);
    q = q + (q >> 8);
    q = q + (q >> 16);
    r = x - 5*q;
    return q + ((11*r) >> 5);

Of course, while being able to divide by 5 is great, ideally we'd want our compiler to come up with a clever algorithm for whatever constant we choose to divide things by. This is considerably more complicated, but there are compiler algorithms which will generate clever bit twiddling code like the above for any constant divisor and inline them right into the code performing the division. Especially for small divisors, the result is drastically more efficient than a call to a general software division function.

A much simpler version to understand in the general case is called exact division. Say you have a number that, for whatever reason, you already know is divisible by 22, and you want to divide it by 22. In this case, the solution is simple. We first shift right by 1; it remains to divide it by 11. Next, drawing from the theory of groups and modular arithmetic, we note for every odd number a less than 2³², there is necessarily a multiplicative inverse b such that ab = 1 (mod 2³²). In other words, multiplying a by b and taking the low 32 bits, which the hardware does already, will result in the number 1. Consequently, x/a = bx/(ab) = bx (since ab=1). There are simple efficient algorithms to find this inverse, such as the extended Euclidean algorithm. Here are the inverses of several small odd integers mod 2³²:

3	2863311531	11	3123612579
5	3435973837	13	3303820997
7	3067833783	15	4008636143
9	954437177	17	4042322161

In our example, to divide, say, 3916 by 11, we'd just compute 3916*3123612579, which when truncated to the low 32 bits by the hardware is 356.

If you're interested in this kind of stuff, check out Hacker's Delight - it covers not only tricks for arithmetic but for other more complex operations such as counting the number of bits set in a word (the population function). I hope you guys found this interesting. As always, please feel free to leave any comments or e-mail me.

Succinct data structures

MSDNArchive — Mon, 05 Dec 2005 21:43:00 GMT

Sorry for the long hiatus, everyone. Today I'm going to talk about succinct data structures, which are informally data structures that use very close to the absolute minimum possible space. This material is largely based on lectures by MIT professor Erik Demaine, as transcribed by 6.897 students Wayland Ni and Gautam Jayaraman (link). Also see additional notes from this year by students Yogo Zhou and Nick Harvey.

One big problem with classical data structures that is too often overlooked is how horribly inefficient their space usage is. For example, suppose you have a linked list of characters. Each node contains a character and a pointer; on a 32-bit platform, each pointer takes 4 bytes, so that's already 5n bytes for n elements. Additionally, it's unlikely with most standard allocators that the 3 bytes following the character will be used; these are normally reserved for padding for alignment. On top of this, unless some kind of memory pool is created for the nodes, most simple allocators will reserve 4 to 12 bytes of space per node for allocation metadata. We end up using something like 16 times as much space as the actual data itself requires. Similar problems appear in many simple pointer-based data structures, including binary search trees, adjacency lists, red-black trees, and so on, but the problem isn't limited to these: dynamic arrays of the sort implemented by Java and .NET's ArrayList classes can also waste linear space on reserved space for growing the array which may never be used.

A natural explanation to this state of affairs would be to say that all of this extra data is neccessary to facilitate fast operations, such as insertion, removal, and search, but in fact this isn't true: there are data structures which can perform many of the same operations using very close to the minimum possible space. What is the minimum possible space? For a list of data, it's simply the space taken by a static array, which is just the space taken by the data itself. But what about, say, a binary tree, or a set of integers?

This is a little trickier, but the idea is to use combinatorics to count the number of possible values of size n, and then take the logarithm of this. Because we need a distinct representation for each value (otherwise we couldn't tell which one we have), this will give us the average number of bits needed to represent an individual value. For example, suppose we have a binary tree with n nodes; we are ignoring data for now and just assuming each node has a left and right child reference. How many bits do we need to store this tree? Well, the number of rooted, ordered, binary trees with n nodes is described by C_n, the nth Catalan number, which can be defined by this formula:

C_n = (2n)! / ((n + 1)! n!)

Now, since Stirling's approximation tells us that log(n!) is nlg n plus some stuff asymptotically smaller than n (written o(n)), where lg is the log base 2, we can see that the logarithm of C_n is asymptotically approximated by:

log C_n ≈ (2n)lg(2n) − (n + 1)lg(n + 1) − nlg(n)
log C_n ≈ (2n)(lg(2) + lg(n)) − 2n lg(n)
log C_n ≈ 2n

In other words, we have to use at least 2 bits per node - this is the best we can possibly do. Surprisingly, there is a very simple encoding that does so. We traverse the tree in preorder, outputting 1 for each node and outputting 0 whenever we hit a null. Here's some C# code which implements the conversion back and forth between this compact bit string and a full in-memory tree structure, preserving both the tree's exact shape and its data:

        class BinaryTreeNode {
            public object data;
            public BinaryTreeNode left, right;
        }

        static void EncodeSuccinct(BinaryTreeNode node, ref BitArray bits, ref ArrayList data) {
            bits.Length++;
            bits[bits.Length - 1] = (node != null);
            if (node != null) {
                data.Add(node.data);
                EncodeSuccinct(node.left, ref bits, ref data);
                EncodeSuccinct(node.right, ref bits, ref data);
            }
        }

        static BinaryTreeNode DecodeSuccinct(BitArray bits, ArrayList data) {
            int bitsIndex = -1, dataIndex = -1;
            return DecodeSuccinctHelper(bits, ref bitsIndex, data, ref dataIndex);
        }

        static BinaryTreeNode DecodeSuccinctHelper(BitArray bits, ref int bitsIndex,
                                                   ArrayList data, ref int dataIndex) {
            BinaryTreeNode result;
            bitsIndex++;
            if (bits[bitsIndex])
            {
                result = new BinaryTreeNode();
                dataIndex++;
                result.data  = data[dataIndex];
                result.left  = DecodeSuccinctHelper(bits, ref bitsIndex, data, ref dataIndex);
                result.right = DecodeSuccinctHelper(bits, ref bitsIndex, data, ref dataIndex);
            } else {
                result = null;
            }
            return result;
        }

The original representation uses about 3 times as much space as the data, while the succinct representation uses only 2n + 1 extra bits, an overhead of 6% on a 32-bit machine. There is some constant overhead, but because this is asymptotically less than n, the data structure is still considered succinct. A number of applications are evident: trees can be compactly stored and transmitted across networks in this form, embedded and portable devices can incorporate pre-built tree data structures in the application's image for fast startup, and so on.

On the other hand, it would seem that you can't really do much with trees in this form - in limited memory environments, you might not be able to expand the tree back out to its full form. An alternative encoding with the same space requirement is to search the tree with breadth-first instead of depth-first search. If we do this, the resulting bit string has some useful properties. Define the Rank and Select functions as:

Rank(i) = the number of 1 bits between indexes 0 and i in the bit string
Select(i) = the index of the ith 1 bit in the bit string

The following can be proven: if we have the breadth-first search bit string for the tree, the children of the node represented by bit i are always represented by bits 2(Rank(i)) and 2(Rank(i)) + 1. The parent of index i is always located at index Select(i/2). Finally, the data for the node at index i is located at index Rank(i) in the data array. Using these lemmas, along with an efficient implementation of Rank and Select, we can traverse the tree and perform all the usual operations on it in this compact form. I won't get into implementing Rank and Select here, but suffice it to say that by adding only O(n/log n) additional bits (which again is asymptotically less than n) they can be implemented in constant time; if the hardware supports a find-first-zero and/or population (bit-count) operation, these can be quite efficient.

This same idea can be applied to all sorts of data structures: for example, there are m choose n (that is, m!/(n!(m-n)!)) different sets of integers between 1 and m, and there are data structures which use bits in the logarithm of this quantity to store such a set. The number of partitions of a set of objects is defined by the Bell numbers, and there are data structures using bits in the logarithm of these numbers. A variety of papers address succinct graph representations, and there are even schemes for representing graphs with specific properties that use no more bits than the logarithm of the number of graphs having those properties.

One of the most interesting succinct data structures, created by Franceschini and Grossi in 2003 in their Optimal Worst-Case Operations for Implicit Cache-Oblivious Search Trees, is a search tree which permits all the usual insert, delete, and search operations in O(log n) time, but with only constant space overhead! It stores the data in an array, and nearly all the information it needs to search it is contained in the order of the elements alone (similar to the heap used in heapsort, but with different operations). On top of this, the structure is cache-oblivious, meaning that it maximally utilizes every level of the memory hierarchy (within a constant factor) without any machine-specific tuning (I'll talk more about cache-oblivious data structures another time). I haven't read this paper yet, but for those who believe data structures are a "solved problem", this kind of result makes it clear that there are many wonders yet to behold.

Finally, I wanted to come back to the problem with ArrayList. Although the simple array-growing scheme used by typical dynamic arrays can waste linear space (sometimes as much unused space is reserved as is used for data), there are alternative dynamic array data structures which remain efficient but reserve only O(sqrt(n)) space, which is asymptotically optimal (there is a sqrt(n) lower bound). These are described in Brodnik et al.'s Resizeable Arrays in Optimal Time and Space.

Thanks for reading, everyone. Please ask any questions that you have, and don't hesitate to e-mail me if you want to.

Persistent data structures

MSDNArchive — Tue, 08 Nov 2005 22:41:00 GMT

When learning to program in a functional language such as Lisp, Ocaml, or Haskell, one of the most difficult aspects of the paradigmic shift is that data in these languages is almost entirely immutable, or read-only. In purely functional programs, there is no assignment, and data structures cannot be modified. But how are we to do anything useful with a data structure we can't modify? The answer is that instead of modifying it, we create new data structures that incorporate the old structure as a part.

The simplest example of this is the singly-linked list. Say you have a pointer to the head of a linked list A. I can take that list and add three new nodes to the beginning to get a new list B. Yet, due to the structure of the list, the old pointer to the first node of the list A will still point to the same data structure as before; you wouldn't even notice I prepended the new elements. We have retained access to both the previous version and updated version of the list simultaneously. Conversely, if I take a pointer to the third node of A, I get a new list C consisting of A with the first two elements removed (called a tail of the list), but existing pointers to A still see the whole list. A data structure with this special property, the ability to retain access to the old version during updates, is called a persistent data structure (not to be confused with a disk-based data structure, which might also be inferred from the word "persistent"). Why might this be useful? There are several reasons:

We can treat all lists as immutable, or read-only. As a result, a function can be sure that no one else is going to modify its list. It can safely pass it to untrusted functions, and annoying phenomena like aliasing and race conditions cannot arise.
Because the data is immutable, the compiler can make assumptions that allow it to produce faster code, and the runtime and garbage collector can also make simplifying assumptions.
It allows us to efficiently "undo" a sequence of operations - just retain a pointer to the old version you want to return to later, and when you're done with the new version, throw it away (assuming you have garbage collection). This is particularly useful for rolling back a failed operation.

On the other hand, other operations on the singly-linked list, such as node removal or insertion in the middle, would overwrite existing pointers and destroy the old version. We can make any data structure (even arrays) persistent by simply making a complete copy of it and then modifying it to perform each operation, but this approach is obviously very expensive in time and space. An active area of programming language research is to find new, more clever data structures that can perform more operations in an efficient, persistent way.

Following is a simple example, taken from Chris Okasaki's excellent and readable PhD thesis, Purely Functional Data Structures (you can also purchase a book version). Say we wish to implement a simple queue. Our queue will support three main operations: add an element to the front, add an element to the back, and remove an element from the front. A singly linked list supports only two of these persistently, adding or removing an element from the front. A simple solution to this problem is to maintain two singly-linked lists, call them L and R. The queue's contents are simply the contents of L followed by the the contents of the reverse of R. For example, if L={1,2} and R={4,3}, then our queue contains the values {1,2,3,4} in that order. To add to the beginning or end of the queue, we simply add to the beginning of either L or R. To remove an element from the front, we remove an element from the front of L. However, it may so happen that L is empty. If this happens, we take R, reverse it, and make it the new L. In other words, we create a new queue {L,R} where L is the reverse of the old R and R is the empty list.

But isn't this operation, reversing a potentially long list, expensive? Not in an amortized sense. Each element added to R will only get moved to L one time ever. Therefore, a sequence of n insertions and n removals will still take only O(n) time, which is an amortized time of O(1) per operation. We say that we can "charge" each of the original n insertions a cost of O(1) to "pay for" the O(n) reversal cost.

Binary search trees and even self-balancing trees like red-black trees can also be made fully persistent: when we add a new node N, we create a new version of N's parent, call it P, with N as a child. We must also create a new version of P's parent with P as a child, and so on up to the root. For a balanced tree, O(log n) space is required for such an insertion operation - more than the O(1) required for an in-place update, but much less than the O(n) required for a complete copy-and-update. These are very useful for building map data structures that need the advantages of persistence.

One key observation to building more clever persistent data structures is that, really, we can modify the old version, as long as it behaves the same as it did before. For example, if we have a binary search tree, we can rebalance the tree without changing its contents; it only becomes more efficient to access. Objects holding existing pointers to the tree won't notice that it's been rebalanced as long as they access it only through its abstract interface and can't touch its internal representation. This can be exploited to maintain an efficient representation as more and more versions are created. On the other hand, this also sacrifices the advantage of avoiding race conditions; we must ensure the data structure is not accessed in an intermediate state where it would not behave the same.

Persistent data structures aren't strictly limited to functional languages. For example, the MOO and ColdMUD virtual environment languages use immutable data structures for built-in string, list, and map types, but have mutable objects. Nevertheless, implementations of persistent data structures are, today, largely limited to functional languages. But they're useful in every language - perhaps one day we will see them in standard libraries. Meanwhile, take a look at Okasaki's thesis, and when thinking about your next design, give a thought to how persistent data structures could simplify your own code.

Software transactional memory

MSDNArchive — Fri, 21 Oct 2005 02:29:00 GMT

Software transactional memory (STM) is a scheme for concurrent programming with multiple threads that uses transactions similar to those used in databases. Today I'll discuss what STM is, how it works, some implementations, and why you should care.

It's well-known that if two or more threads blindly try to access the same data at the same time, a variety of undesirable race conditions can result, situations where the outcome of the program will depend on the order in which the threads access the data. Usually, only one of these orders is desirable. In multithreading, the predominant means of synchronization are locks, which in their simplest form enforce that only one thread can access a particular resource at one time. There are several problems with this approach:

If the locks are too coarse-grained (that is, the code sections they guard are too large), the locked code sections can become a contended bottleneck, leading to decreased concurrency and, on multiprocessor machines, performance closer to that of a single processor machine.
If the locks are too fine-grained, the code can end up spending more time locking, unlocking, and context-switching than it does doing real work.
Locks cause several unpredictable concurrency issues of their own, such as deadlock, livelock, priority inversion, and so on.
In the event of an unexpected error, a lock can be left permanently locked, although C#'s "lock block" helps to prevent this.
Last but not least, locks disturb the flow of the program and force the programmer to think about what concurrent threads could be doing - which programmers happen to be quite terrible at doing. This error-prone process leads to many of the most difficult to reproduce and difficult to fix bugs that we face in our daily lives.

Some of these problems are addressed to some extent by more advanced types of locks such as read-write locks, semaphores, monitors, write barriers, and events, but when you add complexity you add overhead and decrease performance, so we continue to use normal locks more often than anything else.

On the other hand, if you've ever written a SELECT or UPDATE command in SQL, you probably didn't stop to think: what if someone updates the rows I'm querying or updating concurrently? Shouldn't I be locking this table before I mess with it? Should I lock just the rows I'm working with or the whole table? In most applications, query concurrency is simply transparent; the database engine takes care of everything for you. Instead of providing locks, it exposes to the programmer the simpler concept of a transaction. At its heart, a transaction is simply a way of performing a group of operations in such a way that they appear to happen atomically, all at a single instant. The intermediate states cannot be read or updated by any other transaction. This enables you to keep tables consistent with their constraints at all times - one operation can violate a constraint, but both before and after each transaction, it will hold.

In typical implementations of transactions, you have the option of either commiting or aborting a transaction at any time. When you commit a transaction, you agree that the transaction went as planned and make its changes permanently in the database. If you abort a transaction on the other hand, this means that something went wrong and you want to roll back, or reverse, any partial changes it's already made. Good installers use a similar rollback scheme: if some installation step fails and recovery is not possible, it will erase all the cruft it's already put on your machine so that it's just like it was before. If it fails to do this, the machine may become unusable or future installation attempts may fail.

Let's think about how transactions might be useful in concurrent programming. Say you're writing a threadsafe doubly-linked list class. You want to insert a node into the middle of the list. In a single-threaded application, you would just write something like this:

    public void InsertAfter(Node node, T value)
    {
        Node newNode = new Node(value);
        newNode.prev = node;
        newNode.next = node.next;
        newNode.prev.next = newNode;
        newNode.next.prev = newNode;
    }

Unfortunately, the linked list becomes inconsistent between the last two lines (newNode is in the list, but newNode.next.prev is not newNode), and it's impossible to reorder these statements in such a way that the linked list will remain consistent between every pair of statements. You could put a lock around every method accessing the list, but then two threads can't access different parts of the list at the same time, which can be a serious performance problem if it's a large "global" list accessed by many threads at once. If multiple locks became involved, you would also have to be careful to avoid deadlock.

Let's take a different approach based on transactions. We introduce the concept of a version for each object (in our example, each node would be an object). The version is just a counter that starts at zero. Each time an object is modified, its version is incremented. Each time the object is read, the version is (atomically) read along with it. Now we're ready to implement the software transactional memory scheme of Robert Ennals:

Optimistic reads: Every time we read an object, we remember its version at that instant. When we commit the transaction, we verify that the version for every object we read is still the same as it was when we read it. If any of them have changed, we abort and roll back the transaction, undoing any writes it did, and then try the operation over again. If none of them have changed, we know that the concurrent writes would not have altered the transaction's behavior in any case, so we can commit. The assumption here is that the concurrent writes aren't going to tamper with the objects we're reading very often, just other objects we don't care about.
Private revocable writes: Whenever we want to write to a node, we first lock it. The lock prevents others from reading or writing it simultaneously, and the transaction will continue to hold this lock until it commits or aborts. Next, we make a private copy of the node. Any changes we want to make we make to this private copy. By ensuring that each object is accessed through only one reference (which may itself be referenced by multiple references), we can atomically swap out the old copy for the new, updated copy when we commit the transaction. To abort the transaction, we simply release all write locks and throw away all the private copies. Notice that if a transaction locks a node for writing and later aborts, transactions which read that node can still go on to succeed.
Deadlock detection: When we're automatically locking things all over the place, deadlock is sure to happen sooner or later. To deal with this, we periodically interrupt the program (much as a garbage collector does) and scan for a cycle of objects that are deadlocked, each waiting on a lock held by the next. We then force one of them to abort, releasing its locks and breaking the cycle. With this system in place, deadlocks can never occur. Databases with locking use a similar system to prevent deadlock.
Edge cases: If a transaction which has read inconsistent state tries to commit, it will fail, but it's also possible that such a transaction could go into an infinite loop or have a fatal exception because it wasn't written to deal with this situation (nor should it be). Again, the runtime is responsible for detecting and aborting such transactions.

Now, all this copying, swapping, and retrying of aborted operations looks expensive. However, the truth is that so much concurrency is gained from a transaction-based system that they are often much more efficient than lock-based systems even for a small number of processors, as you can see in the Performance section of Ennals' paper. The reason for this is that simple locking is often too pessimistic - it locks an entire data structure when really different threads could be safely operating on different parts of it concurrently, as long as you have a roll back mechanism ready in case a conflict arises.

So it's really a win-win - you can write your program almost as though it were single-threaded, and it runs faster! Even on a single processor, the penalty is not so severe as to outweigh the benefits to the programmer (who is, after all, more expensive than both software and hardware). If we return to our linked list example, we can simply take the single-threaded solution, put each method inside a transaction, and we're done.

This all sounds nice in theory, but few of us are qualified to write and use a production-quality STM library after a quick read of a paper. Luckily, there are already some good implementations floating around that you can try:

Robert Ennal's Lightweight Transaction Library, a C library written for UNIX and recently re-released under the modified BSD license (which allows commercial use as well as commercial use of derivative works such as ports). Someone should really port this to a .NET library.
Microsoft Research's C# Software Transactional Memory (SXM). (restricted to noncommercial)

Here are some more papers on the topic:

The original paper on software transactional memory.
Microsoft Research's Composable Memory Transactions.
Software Transactional Memory for Dynamic-Sized Data Structures, which introduced the obstruction-free DSTM implementation.
The Wikipedia article (which needs some expansion).
A small bibliography.

Give them a try - in 5 or 10 years, locks might just be a footnote in a Systems textbook.

Non-nullable types

MSDNArchive — Mon, 10 Oct 2005 21:39:00 GMT

If you write programs in C, C++, Java, or C#, you've gotten used to having the null value around. The null value is a special reserved reference (or pointer) value indicating that a reference does not refer to any object. It's useful for constructing a variety of data structures, but it's also a notoriously common source of bugs, because occasionally values that "should never" be null turn out to be, unexpectedly. On modern machines this type of runtime error can be detected quickly if it occurs, typically by using the value zero for null and unmapping the first page of memory, but it would be much nicer to prevent this type at error at compile-time.

Impossible, you say? Not at all. In fact, functional languages in the ML family such as Standard ML, Ocaml, and Haskell have always made a distinction between nullable types (types that can have a reserved null value, such as "int option") and non-nullable types (such as "int"). There's no reason for every type to include, by default, a feature you may not use. There's nothing magic about this; it doesn't make a program any more difficult to type check if you say some values can be null and other values can't.

We could make a similar distinction in a language such as C#. In fact, someone already has. One of the most exciting research projects at Microsoft is an extension of C# called Spec#. Spec# has a lot of truly cutting-edge features that I won't get into right now, but the one simple feature that I do want to mention is non-nullable types. These allow you to affix an exclamation mark (!) onto any reference type, yielding a new type that cannot be null but is otherwise the same. For example, consider this simple string comparison method:

        public int StringCompare(string left, string right)
        {
            for(int i=0; i<Math.Min(left.Length, right.Length); i++)
            {
                if (left[i] != right[i])
                {
                    return left[i] - right[i];
                }
            }
            return right.Length - left.Length;
        }
    
        public int Foo(string name)
        {
            return StringCompare(name, name + "bar");
        }

If I enter this into a Spec# buffer, it underlines the references to "left" and "right" with red lines to indicate that I may be dereferencing a null value, which is a compile-time error (not just a warning). The solution is to change the signature of StringCompare to the following, indicating that it only accepts non-null strings:

        public int StringCompare(string! left, string! right)

If I do this, however, I get an error in Foo() that it's passing a possibly-null string where a non-null string is expected (a type error). Supposing I want Foo to accept null values, I can add the appropriate null check and the error goes away:

        public int Foo(string name)
        {
            if (name == null) return 0;
            return StringCompare(name, name + "bar");
        }

Even more importantly, Spec# ships with improved signatures for a large number of .NET Framework methods. You no longer have to worry about NullReferenceExceptions occurring when you pass null to a Framework method that doesn't accept it - you find out at compile-time. They generated most of these signatures automatically by reflecting on the CLR and looking for an argument being tested for null followed by a throw of NullReferenceException, but ideally the CLR team would include the proper types from the beginning.

Non-nullable types are also present in other extension languages and tools. Cyclone, the safe C dialect, adds non-nullable pointers to C (along with several other types of pointers, see Pointers). ESC/Java adds non-nullable pointers, along with many other extensions, to Java (see non-null pragma). Microsoft has a well-known internal static checking tool called PREfix which allows the specification of non-null types in C++ code using the SAL_notnull SpecStrings annotation. This feature is so easy to add and so useful that just about every new language does so.

So why have we continued to battle null values for all these years when there's an easy-to-implement, well-known compile-time solution? Got me, but I wouldn't mind seeing this feature in C# 3.0.

Secret sharing

MSDNArchive — Mon, 03 Oct 2005 01:31:00 GMT

One of the most difficult problems in cryptographic key management is keeping a secret key safe from both compromise and loss. If you don't make enough backups, the key might be destroyed in a hardware failure or natural disaster. But if any backup is compromised, the key is compromised.

Rather than invent new tools, one might try to solve the problem with encryption. Just encrypting the key normally is not helpful, since then you have the same problem with the key needed to decrypt it. Instead, if you encrypt your key using multiple keys, then store these keys in different locations, you need all the keys to decrypt it. This protects against compromise, but if any key is lost, so is the original key. To fix this, you'd have to encrypt the original key many times, each time using just some of the keys. For example, suppose you had 5 keys and you wanted to be able to decrypt using any 3 of them. Then you'd need to perform all the following encryptions:

E₁(E₂(E₃(K)))
E₁(E₂(E₄(K)))
E₁(E₂(E₅(K)))
E₁(E₃(E₄(K)))
E₁(E₃(E₅(K)))
E₁(E₄(E₅(K)))
E₂(E₃(E₄(K)))
E₂(E₃(E₅(K)))
E₂(E₄(E₅(K)))
E₃(E₄(E₅(K)))

The number of encryptions needed grows very quickly. It also takes quite a bit of time and space to store all these multiply-encrypted keys. Let's look at some other solutions.

One solution is to simply break the key file into two parts, make a copy of each, and store the four parts in four locations. Then, even if any one is compromised, the key is not revealed. Unfortunately, this scheme is not that helpful, because just knowing half of someone's private key is enough to drastically decrease the time needed for a brute-force search to crack the remainder of the key.

Here's another simple solution: create a one-time pad (a file of random data) the same size as the private key, then XOR them together. Put the result in one place and the pad in another, and destroy the original key. We call each file a "share" or "shadow". Later, you can XOR the two shares together to retrieve the original key. Since the bytes of each file are random, you can be assured that neither piece will reveal any information about the key by itself. We can extend this scheme to n pieces by just creating n - 1 one-time pads, adding them together, and subtracting this from the key data to obtain the nth share. Unfortunately, if even one piece is lost, so is the key. This is called an (n,n) secret sharing scheme, because we have n shares and need n of them to reconstruct the secret.

Here's another more useful scheme: let the private key be represented as the number S. Choose a random line in the plane passing through the point (0, S). Next, choose n random points on that line as your n shares. Now, if we know any two of the shares, we know the line, and so can find S. If we only know one share, then we just know one point, and for every T there is a line passing through (0, T) and that point; hence, the secret could be anything at all. This is called an (n,2) secret sharing scheme, because we have n shares and need 2 of them to reconstruct the secret.

Particularly useful for individuals is the (3,2) scheme - you can keep one share on your person on a floppy or USB pen drive at all times, one share on your hard drive, and one share backed up at a secure remote location. When you need to use the key, you recreate it briefly using the shares on your hard drive and your person, then destroy it. If you lose the floppy or your hard disk crashes, you can recreate the key from the remaining shares (and make yourself 3 new shares). If any one share is compromised, the key is not compromised. If you know of a compromise, you can recreate the key, create 3 new shares, and destroy the old ones, preventing the attacker from getting a second share. The new shares are unrelated to the old ones because a different random line through (0,S) is used.

In 1979, Shamir generalized these ideas even further to create a general (n,k) secret sharing scheme for any n, k. The idea is to choose a random degree-k polynomial passing through (0,S), where S is the secret. Then choose n random points on this curve. Because there is exactly one degree-k polynomial passing through any given set of k points, any k of these points define the curve and allow us to retrieve the secret. If we have only k - 1 points, then we don't learn anything about S, because there is a polynomial passing through (0,T) and those k - 1 points for every T. To make the notion of choosing random polynomials and random points precise, Shamir works with polynomials modulo a prime number p, choosing polynomial coefficients and point coordinates randomly from the integers in [0, p - 1].

Secret sharing has more applications than just protecting private keys. It's also useful in situations where we want to entrust a secret to a group of people but want to ensure that several of them must cooperate to access the secret. We don't need to trust them all individually, nor do they all need to cooperate, so the secret is still available if members go missing or refuse to cooperate.

Finally, I wouldn't just tell you all this without providing some working code you can use. You can download an excellent implementation of the scheme for UNIX called secsplit (direct download). I don't know of any free .NET implementation, but I'm sure that porting this 600-line program into any environment wouldn't be difficult.

Some more resources on secret sharing:

Lectures notes from Cornell's CS 513
Wikipedia article
Adi Shamir's "How to Share a Secret"
RSA FAQ

Custom building and code generators in Visual Studio 2005

MSDNArchive — Tue, 13 Sep 2005 23:15:00 GMT

I'm a fervent fan of using code generator tools wherever possible to make your life easier. Although they come with issues related to effective building, diagnostics, and debugging, the amount of value they add to your application is immense: they can eliminate entire classes of potential bugs, save you a great deal of effort and time, make your module much easier to extend and maintain, and even yield runtime performance gains. Among the most frequently used code generator tools are the lexer and parser generators GNU Flex and Bison, based on the classic generators lex and yacc. Although I'll get into the details of how to use these tools effectively at a later time, today what I want to show you is a practical example of how to use the new custom build features of Visual Studio 2005 to effectively incorporate a code generator into your automatic build process.

Here is a very simple Bison grammar for evaluating arithmetic expressions involving addition, subtraction, and multiplication:

/* example.y */
%{
#define YYSTYPE int
%}
%token PLUS MINUS STAR LPAREN RPAREN NUMBER NEWLINE
%left PLUS MINUS
%left STAR

%%
line : /* empty */
     | line expr NEWLINE      { printf("%d\n", $2); }
expr : LPAREN expr RPAREN     { $$ = $2; }
     | expr PLUS expr         { $$ = $1 + $3; }
     | expr MINUS expr        { $$ = $1 - $3; }
     | expr STAR expr         { $$ = $1 * $3; }
     | NUMBER                 { $$ = $1; }
     ;

%%

int yyerror (char const *msg) {
    printf("Error: %s\n", msg);
}

int main() {
    printf("%d\n", yyparse());
    return 0;
}

The Flex lexer used by this parser looks like this:

/* example.lex */
%{
#include "example.parser.h"
%}
%option noyywrap

%%

[ \t]+    { /* ignore whitespace */ }
"("       { return LPAREN; }
")"       { return RPAREN; }
"+"       { return PLUS; }
"-"       { return MINUS; }
"*"       { return STAR; }
\n        { return NEWLINE; }
[0-9]+    { yylval = atoi(yytext); return NUMBER; }
.         { printf("Invalid character '%s'", yytext); }

%%

If we were writing this parser at a UNIX command line, we might generate the source files and compile the result using this sequence of commands:

bison -v -d example.y -o example.parser.c
flex -oexample.lexer.c example.lex
gcc -o example example.lexer.c example.parser.c

Now say you wanted to build the same application for Windows using Visual Studio 2005. The tools are available on Windows (see flex for Win32, bison for Win32), and you could simply run the same first two commands at the command-line and then build the resulting source files in Visual Studio. However, this simple approach sacrifices many of the advantages that Visual Studio provides for its built-in source file types: it doesn't rebuild the generated source files as needed, it doesn't allow you to jump to errors that occur during generation, and it doesn't allow you to configure build options using a nice GUI. Let's see how we can reclaim these advantages for Flex and Bison files.

Creating a simple custom build type

Our first goal is simply to be able to build Flex and Bison files. First, use the Flex and Bison setup binaries from the links above to install the tools. A bin directory will be created in the installation directory. Add this to your system path. You should be able to execute both flex and bison from a command prompt without specifying a path.

Next, we create a new C++ console application. Uncheck the option to use precompiled headers - I'll explain how to use these with Flex and Bison later. Remove the main source file created for you by the wizard. Next, right-click the project in Solution Explorer and choose Custom Build Rules. The following dialog appears:

A build rule establishes how to build a file of a particular type. A group of related build rules are stored in a build rule file, which can be saved, distributed, and reused in many projects. We'll start by creating a new build rule file for our Flex and Bison rules:

Click New Rule File.
Enter "GNU Tools" for Display Name and File Name.
Choose a suitable directory for the build rule file. If it asks you if you want to add the directory to your search path, say yes.

Now we'll create a build rule for Bison files:

Click Add Build Rule.
Enter the following values:
1. Name: Bison
2. File Extensions: *.y
3. Outputs: $(InputName).parser.c;$(InputName).parser.h
4. Command Line: bison -d [inputs] -o $(InputName).parser.c
5. Execution Description: Generating parser...
Click OK twice, then check the box labelled "GNU Tools" and click OK.
Add the example.y file above to your project. Right-click on the file and choose Compile. You should receive no errors.
Create a new file folder under the project called "Generated Files". Add the existing file example.parser.c to this folder.

If you build now, you should receive only an error complaining that yylex() is undefined. Now, go back to Custom Build Tools and click Modify Rule File on GNU Tools. Create a rule for Flex:

Click Add Build Rule.
Enter the following values:
1. Name: Flex
2. File Extensions: *.lex
3. Outputs: $(InputName).lexer.c
4. Command Line: flex -o$(InputName).lexer.c [inputs]
5. Execution Description: Generating lexer...
Click OK three times.
Add the example.lex file above to your project. Right-click on the file and choose Compile. You should receive no errors.
Add the existing file example.lexer.c to your project.

If you build now, you should receive no errors and be able to run the application successfully. Now in any project you can simply check the "GNU Tools" box, add the .lex and .y files to your project, and build. What happens if you modify the example.y and build? It runs Bison again and recompiles example.parser.c, because it was regenerated, and example.lexer.c, because it includes a header file that was regenerated. If we modify the .lex file, Flex is rerun and example.lexer.c is recompiled, but example.parser.c is not rebuilt. If you had a larger parser, you'd appreciate how much time this incremental rebuilding saves you.

Improving diagnostic support

Delete one of the "%%" marks in the .y file and build. Unsurprisingly, Bison fails. However, the Error List tells you no more than this. It'd be more helpful if you could find out what errors the tool produced. If you look at the output window, Bison did produce some errors, but if you double click on them to visit the error location, it just takes you to the top of the file. What gives?

The reason for this is that Visual Studio only recognizes one error format, that used by its own tools. Here's an example:

c:\myprojects\myproject\hello.cpp(10) : error C2065: 'i' : undeclared identifier

Bison doesn't output errors in this format, and so they aren't parsed. Flex uses yet another different format. What to do? The simplest way to deal with this is to invoke a simple script on the output of the tools as part of the build rule which parses the output and converts it to the desired format. You can write this script in any language; I wrote them in C# using the .NET Framework's regular expressions. Here's what I wrote inside the Main() function for the Bison converter tool (error checking and such omitted):

string line;
while ((line = Console.In.ReadLine()) != null)
{
    Match match = Regex.Match(line, "([^:]+):([0-9]+)\\.[^:]*: (.*)");
    if (match != null)
    {
        Console.WriteLine("{0}({1}): error BISON: {2}",
                          Path.GetFullPath(match.Groups[1].Value),
                          match.Groups[2].Value, match.Groups[3].Value);
    }
    else
    {
        Console.WriteLine(line);
    }
}

I deploy the binary, say it's called BisonErrorFilter.exe, to the same directory as bison.exe. I then change the Command Line of the Bison build rule to the following (click the arrow in the right of the field to access a multiline text box):

bison.exe -d [inputs] -o $(InputName).parser.c > bison.err 2>&1
BisonErrorFilter < bison.err

If you compile the .y file now, any errors should appear in the error list, as desired, and you can double-click them to visit their locations. I wrote a similar script for the lexer output. Be careful when doing this, though, because if you miss any errors, Visual Studio might look at the error return of the last command and interpret it as success. A better way to do this would be to wrap the tool in a script that passes its arguments to the tool, collects the tool's output and return code, converts and prints the output, and then returns its return code.

I haven't figured out how, but I believe it's possible to also create custom help entries for each error message, then have the filter tool produce the right error code for each one. This way, users can get help for each error individually by just clicking on it and pressing F1.

Properties

Properties enable you to control how the command-line tool is executed directly from the properties page for each individual file you wish to build with it. Let's start with a simple example: a handy lexer switch is -d, which prints out an informative message each time a token is recognized. We don't want it on all the time, and certainly not in release mode, but it'd be handy to be able to turn on and off as necessary.

To create a property for this, first return to the lexer build rule. Then follow these steps:

Click Add Property.
Choose Boolean for the User Property Type.
Enter the following values:
1. Name: debug
2. Display Name: Print debug traces
3. Switch: -d
4. Description: Displays an informative message each time a token is recognized.
Click OK. Then, add [debug] right after "flex" in the Command Line field.
Click OK three times.
Right-click on example.lex in Solution Explorer and choose Properties.
In the left pane, click the plus next to Flex. Click General.
You'll see your property. Click on it and its description will appear at the bottom. Set it to Yes.
Click Command Line in the left pane. You'll see that the -d flag has been added.
Click OK and build.
Run the app and type an arithmetic expression. You'll see trace messages.
View the project properties. You'll see that it now has a Flex node also. Here you can set the default settings for all files of that type in the project which don't have specific overriding settings set.

Adding more properties is just as simple. You can go through the man page for the tool and add properties for each switch, using the Category field to group them into categories. You can use the other property types for switches accepting arguments. If you want, you can create a detailed help file with additional explanation and examples for each switch. When you're done you have an impressive looking property sheet for your files reminiscent of those for built-in types:

You can also set different settings for debug and release builds. For example, for Flex, it's good to set table size to slowest and smallest for the Debug version, to speed up compilation, and to set it to the recommended full tables with equivalence classes for the Release version, which is a good tradeoff of table size and speed.

Finally, once you're done adding all the properties you like, you can take the resulting .rules files and give it to everyone on your team, or distribute it on a website, so that everyone can easily integrate the tool into Visual Studio. Perhaps eventually tools like Flex and Bison will ship with a .rules file.

Conclusion

In Visual Studio 2003 you would have had to write a plug-in to come close to achieving this level of integration with a third-party tool. Although it has limitations, I hope the problems solved by these new features help encourage you to incorporate more tools and code generation into your regular development. Now that you know how to use Flex and Bison from the Visual Studio IDE, next time I'll talk about how to use the tools themselves, going through some of the development and debugging processes that a grammar developer goes through, and show you some similar tools for other .NET languages. Thanks for reading, everyone.

Modular arithmetic and primality testing

MSDNArchive — Wed, 07 Sep 2005 20:29:00 GMT

Number theory is, roughly speaking, the study of properties of integers. Often a problem which is easy for real numbers, such as factoring or linear programming, seems to be considerably more difficult when restricted to integers (in fact, integer programming is NP-hard). Much of the focus of modern number theory is on solving these problems efficiently.

The practical importance of number theory was greatly increased by the introduction of RSA cryptography, an unprecedented system whereby information could be securely transmitted without first transmitting the key and risking its interception. Central to RSA is the ability to efficiently generate a semiprime, which is simply a product of two very large primes of roughly equal magnitude. What's useful about these numbers is that we can generate them efficiently, yet given just the semiprime, determining the two prime factors is a hard, unsolved problem. It turns out that there are a lot of primes to choose from: the probability that a randomly chosen k-bit number is prime is between 0.69/k and 0.87/k. For example, about 7-9% of 10-bit numbers are prime. Consequently, we can find a k-bit prime number in an expected O(k) random guesses. The only remaining problem is to find an efficient algorithm for testing these guesses for primality.

Brute force primality tests such as simply dividing the number by all values less than the square root are ineffective for large numbers. To get around this, we'll need an important tool called modular arithmetic. If you've written C, you may be surprised to know that you're already familiar with modular arithmetic - the "wrap around" behavior of k-bit integers is an implementation of arithmetic mod 2^k. For example, this piece of code computes the value 17:

    unsigned char c1 = 177, c2 = 96;
    unsigned char sum = c1 + c2;

In mathematical notation, we would say:

177 + 96 ≡ 17 (mod 256)

This tells us that the numbers 177 + 96 and 17 differ by a multiple of 256, or in other words have the same last 8 bits. Multiplication mod 256 is likewise similar to multiplication of unsigned chars in C:

    unsigned char c1 = 177, c2 = 23;
    unsigned char product = c1 * c2;

177 × 23 ≡ 231 (mod 256)

Modular arithmetic works exactly the same for other moduli than 256; the only difference is the number where it wraps around. For example, if you compute 10 + 10 mod 17, you get 3.

The interesting thing about modular arithmetic is that it can be shown that the values form a commutative ring. This means, among other things, that:

Addition and multiplication are commutative (you can reverse their operands without changing the result).
The associative property holds for both addition and multipliation: (177 × 23) × 45 gives the same thing as 177 × (23 × 45), even if computed with unsigned chars.
The distributive property holds: (177 + 23) × 45 is equal to (177 × 45) + (23 × 45).

If none of these surprise you, this might: if the modulus is a power of 2, as in machine arithmetic, every odd number m has a multiplicative inverse. An inverse is simply a number n such that m × n = 1. What good is this? Well, suppose you want to divide a number k by 7 on a machine with no hardware division. You know that k is divisible by 7. The inverse of 7 mod 256 is 183. Because 7 × 183 = 1, we have 7 × 183 × k = k. Divide both sides by 7, and we get 183 × k = k/7. In other words, multiplying by a number's inverse is the same as dividing by that number, provided that k is evenly divisible by the number.

Now we come back to our original problem: how can modular arithmetic help us determine primality? Well, it's not hard to show that if p is prime, then:

For all a with 0 < a < p, a^p ≡ a (mod p).

This is called Fermat's little theorem. What's useful about it is not only that it's true for all primes, but that if you find an a such that it does not hold, you have proven that the number p is composite. This can be checked efficiently because there is an efficient algorithm for computing large powers. This is simply an existence proof - it tells us nothing about what the factors are. It can be shown that for most composite numbers, if you just select several a at random and try them, one is likely to fail the test. Unfortunately, there are an infinite number of special numbers called Carmichael numbers for which the above result holds for all a, even though they are not prime.

To get around this, we design a new, slightly more complicated test which is not susceptible to this problem. I take this from the excellent book Prime Numbers: A Computational Perspective, by Richard Crandall and Carl Pomerance. Suppose p is an odd prime, and that the binary representation of p - 1 is the odd number t followed by s zeros. Then one of the following is true for every a with 0 < a < p - 1:

a^t ≡ 1 (mod p)
a^{t << i} ≡ p - 1 (mod p) for some i with 0 ≤ i < s

Above "<<" denotes the shift-left operator. It can be shown that for any odd composite number > 9, both of the above will fail for at least 3/4 of the a between 1 and p-2. We can use this to design a simple probabilistic algorithm:

Choose an a at random with 0 < a < p-1.
Check if one of the above formulas holds. If not, p is composite.
If we have iterated at least T times, claim that p is prime. Otherwise, go back to step 1.

Here's an (untested) piece of C# code which implements this simple algorithm, assuming that BigInteger is a suitable arbitrary-precision integer type (the Framework provides no such type):

    public bool IsProbablePrime(BigInteger p, int T) {
        int s = 0;
        BigInteger t = p - 1;
        while (t.IsEven()) {
            s++;
            t >>= 1;
        }
        for (int repeat = 0; repeat < T; T++) {
            BigInteger a = BigInteger.Random(1, p - 2);
            if (!PassesPrimalityTest(a, p, s, t)) {
                return false;
            }
        }
        return true;
    }

    public bool PassesPrimalityTest(BigInteger a, BigInteger p, int s, BigInteger t) {
        BigInteger b = BigInteger.ModPower(a, t, p); // b = a^t mod p
        if (b == 1 || b == p - 1) return true;
        for (int j = 1; j < s; j++) {
            b = BigInteger.ModPower(b, 2, p);
            if (b == p - 1) return true;
        }
        return false;
    }

Because each trial is independent, the chance of erroneously claiming that p is prime is 1/4^T, which is ridiculously tiny even for small T, like T = 20. We can now use this to quickly locate large prime numbers and generate semiprimes for RSA cryptography.

Unfortunately, these tests identify composite numbers but reveal nothing about their factors. This makes them useless for factoring numbers. Unlike many other hard problems, the factoring problem is believed to not be NP-complete, and so may very well have an efficient algorithm that no one has found yet. Such an algorithm would make RSA encryption useless and win you $625,000. Combinations of advanced techniques have managed to factor very large numbers, but only at an enormous expense in computing time and hardware. This is one of the most active current areas of research in number theory and computer science.

The point location problem

MSDNArchive — Fri, 02 Sep 2005 04:05:00 GMT

Computational geometry is a field that studies efficient solution of geometric problems, which are critical in mapping, manufacturing, and particularly graphics applications. If you find data structures and algorithms interesting, it's likely you'd also find computational geometry interesting, because it often combines and adapts well-known data structures and algorithms in novel ways to solve natural problems that are easy to explain. As an example of this, I'm going to discuss the point location problem.

Imagine you have a political map of Africa with a red dot somewhere on it. Someone asks you, what country is that dot in? This is a simple question for even a child to answer, but how do we make a computer answer the same question? We can represent the map as a planar decomposition, meaning a division of the plane into polygonal regions. We store the vertices making up the boundary of each region. Stop and think for a little bit about how you might solve this.

Perhaps the simplest approach is to exploit the cross product. Suppose the vertices of the region are listed in clockwise (or counterclockwise) order. Compute the cross product of the vector from vertex k to vertex k − 1 with the vector from vertex k to the point being located. The point lies in the region if and only if all of these cross products (for all k) point in the same direction. Since we visit each vertex twice and the cross product takes constant time, this requires O(n) time, where n is the total number of edges.

Another simple solution that generalizes better to higher dimensions is to start at the point and draw a line emanating from it. The first time you hit an edge, you know that edge is on the boundary of the correct region. The edge may be on the boundary of two adjacent regions, but comparing the slope of the line to the slope of the edge will tell you which one is correct. To actually compute this, we define the line parametrically with an equation like this, where (x₀, y₀) is the point being located:

(x, y) = (x₀, y₀) + t(1, 1)

Then for each edge, we compute the t value at which the line intersects that edge, if any, and choose the edge with the smallest positive t value. This algorithm is also O(n).

These algorithms are fine if you're only locating one point, but if we wish to locate many points we can do better by creating a data structure which facilitates quick point location. I'll begin by solving a much simpler problem: say your map only has rectangles in it, and these rectangles are arranged in vertical "strips" as shown to the right. To solve this problem, we can construct a balanced binary search tree. This tree will first compare the x coordinate of the point against the x coordinates of the vertical lines to determine what strip of rectangles the point falls in. It then compares the y coordinate of the point against the y coordinates of the horizontal lines to determine what rectangle it falls in. Each comparison eliminates about half of the remaining rectangles from consideration, so this takes only O(log n) time in all for n rectangles. The search tree takes O(nlog n) time to construct and takes O(n) space, since there's at most one node for each edge and there are 4n edges.

We generalize this by allowing each "strip" to contain not just rectangles but trapezoids, as shown below (both rectangles and triangles are considered to be trapezoids here). This is called a trapezoidal map. This makes the second part of the search only a little more complicated: instead of comparing the y coordinate to a fixed value, we choose an edge and use a cross product to check which side of the edge the point falls on.

Finally, take any map, and draw a vertical line through each vertex; in the image to the left the gray vertical lines were added in this manner. You'll end up with a trapezoidal map. Since we can determine what region of this map contains the point in O(log n) time, and each of these regions is part of only one region in the original map, we've solved the problem for general maps in O(log n) time.

Unfortunately, adding all these vertical lines also creates a lot of new regions, potentially as much as squaring the number of regions. Search time is still good, since log(n²) = 2 log n, but the storage required can jump from O(n) up to O(n²). This is worst-case, but even in typical practical cases it can require a lot of space, and construction time for the data structure also increases. We want to retain a quick search capability but without having quite so many regions.

To achieve this, notice that we don't actually need to draw a complete vertical line through every point to divide the map into trapezoids. All we really have to do is draw a line up and down from each vertex until we reach the edges above and below. It can be shown that the number of trapezoids is now no more than about three times the number of original regions, so the worst-case space has been cut back down to O(n). The search procedure is similar: at each node we either compare the x coordinate of the point to that of a vertical line segment, or we check on which side of an edge the point lies. Either way we eliminate about half the remaining regions. The search "tree" is no longer a tree but a directed acyclic graph, because a single trapezoid may fall on both sides of a vertical line segment. See a demonstration.

Of course I've left a lot of details out here, like how to ensure the search tree remains well-balanced during construction, what to do if some points have the same x coordinate, and so on, but I hope this gives you a taste of how clever algorithms can efficiently solve geometric problems.

The Visitor pattern and multiple dispatch

MSDNArchive — Tue, 30 Aug 2005 04:27:00 GMT

Today I'll talk about multiple dispatch, a programming language feature that increases flexibility of method calls and eliminates the need for awkward pattern constructions, at the cost of some additional complexity in the dispatch algorithm.

The Visitor pattern

Those who have read the Gang of Four's seminal Design Patterns may recall the Visitor pattern. Its purpose is to provide new operations on a class or set of related classes without actually altering the classes involved. Design Patterns provides a real-world example in the document-view paradigm, but for simplicity I'll provide a simpler example here based on some notes by Corky Cartwright.

Say you wish to represent an arithmetic expression using integers and addition such as (1 + 5) + 2. In an object-oriented language, a natural way to do this is to create an abstract base class "Expression", and then create some related subclasses:

ConstantExpression: Represents a constant number, like 2.
SumExpression: Represents the result of adding two Expression objects.

The objects form a tree, where the leaves are ConstantExpression objects. To evaluate an expression represented in such a form, it suffices to create a virtual Evaluate() method which each subclass implements appropriately:

abstract class Expression {
    public abstract int Evaluate();
}

class ConstantExpression : Expression {
    int constant;
    public override int Evaluate() {
        return constant;
    }
}

class SumExpression : Expression {
    Expression left, right;
    public override int Evaluate() {
        return left.Evaluate() + right.Evaluate();
    }
}

This may solve the problem if all you want to do is evaluate expressions, but there are all sorts of operations you can do on arithmetic expressions. Do we really want to add a virtual method for every one of them? We certainly don't want clients of the class to be forced to extend the classes just to add application-specific operations. The Visitor pattern solves this problem by replacing the Evaluate() method with an Accept() method that takes as a parameter a visitor object. The visitor object has a method for each type of expression, and Accept() invokes the appropriate one. Here's how we would evaluate expressions using the Visitor pattern:

abstract class Expression {
    public abstract object Accept(Visitor v);
}

class ConstantExpression : Expression {
    public int constant;
    public override object Accept(Visitor v) {
        return v.VisitConstant(this);
    }
}

class SumExpression : Expression {
    public Expression left, right;
    public override object Accept(Visitor v) {
        return v.VisitSum(this);
    }
}

interface Visitor {
    object VisitConstant(ConstantExpression e);
    object VisitSum(SumExpression e);
}

class EvaluateVisitor : Visitor {
    public object VisitConstant(ConstantExpression e) {
        return e.constant;
    }
    public object VisitSum(SumExpression e) {
        return (int)e.left.Accept(this) + (int)e.right.Accept(this);
    }
}

We evaluate an expression e using (int)e.Accept(new EvaluateVisitor()), and adding new operations is easy; we just create new implementors of Visitor. Clients can create application-specific algorithms without touching the original classes. As an added benefit, the visitor class can keep state internal to the algorithm; for example, this visitor counts and stores the number of constants in an expression:

class CountConstantsVisitor : Visitor {
    public int Count;
    public object VisitConstant(ConstantExpression e) {
        Count++;
        return null;
    }
    public object VisitSum(SumExpression e) {
        e.left.Accept(this);
        e.right.Accept(this);
        return null;
    }
}

But the wary designer will note that there is a maintainance problem here: if we add a new subclass of Expression, like ProductExpression, the Visitor interface and all classes implementing it have to be updated. If the Expression class rarely changes, this is fine, and may even be good: it forces the visitor classes to consider the new cases. If it's rapidly changing, however, we definitely have a problem.

Double dispatch

At this point we take a slight detour and consider another simple example. Polymorphism is the critical ability of an object-oriented method call to invoke different methods at runtime based on the type of the object, like Evaluate() in our first example. If you've ever dealt with operator overloading in C++ or C#, though, you may come up against the limits of polymorphism. If I write an operator+ for a class Complex for complex numbers, subclasses can override it in a polymorphic way - but only if the complex number is on the left. If the Complex object is on the right-hand side, there is no way to ensure the correct method is invoked polymorphically short of testing the type and invoking it yourself.

How cruel we are to the unprivileged right-hand sides! The built-in C operator + produces different types depending on both the type of its left and right argument; both a double plus a float and a float plus a double produce doubles. We want something similar for our own operators, the ability to polymorphically choose a method based on both the type of its left argument and its right argument. Then we could write code like this:

Complex operator+(Complex c, int i) { ... }
Complex operator+(ComplexSubclass c, int i) { ... }
Complex operator+(int i, Complex c) { ... }
Complex operator+(int i, ComplexSubclass c) { ... }

The right method will be chosen at runtime based on the types of both sides; the most specific (most derived) type is always chosen. This is called double dispatch, because we select the method ("dispatch" the message) based on the type of two arguments.

Unfortunately, this new flexibility also introduces some new problems. Suppose we have these two methods:

Complex operator+(ComplexSubclass s, Complex c) { ... }
Complex operator+(Complex c, ComplexSubclass s) { ... }

Now, suppose we try to add two objects of type ComplexSubclass. Neither of the above looks more appropriate than the other. There are various ways of dealing with this; for example, we can give the first argument precedence over the others, we can throw an error at runtime, or we can make the compiler require us to implement a method taking two ComplexSubclass objects. For this article, we'll sidestep the issue by saying we won't write any code that will invoke this dilemma.

Now, I'm sure you're wondering what any of this has to do with the Visitor pattern. The answer is that this feature allows us to drastically simplify our Visitor pattern example while retaining all of its benefits:

abstract class Expression { }

class ConstantExpression : Expression {
    public int constant;
}

class SumExpression : Expression {
    public Expression left, right;
}

class EvaluateVisitor {
    public int Visit(ConstantExpression e) {
        return e.constant;
    }
    public int Visit(SumExpression e) {
        return Visit(e.left) + Visit(e.right);
    }
}

Now, this isn't valid C# of course, since C# doesn't support double dispatch, but you can see the effect it would have if it did: the Visitor interface is gone, the Accept() method is gone, and the calls in the last Visit() method are simplified. Best of all, it solves our dilemma of a rapidly changing Expression class hierarchy, because we can add a "default case" method that deals with any expression types we're not familiar with:

    public int Visit(Expression e) {
        throw new Exception("Unsupported type of expression"); // or whatever
    }

One handy trick you can do in Java and .NET languages is to have a method that examines the type of its argument and then uses reflection to pass it on to the correct overload, simulating polymorphism. This is much slower than native double dispatch support, but achieves much of the same effect.

If you've studied how polymorphic methods are dispatched, you're probably most familiar with the "vptr" and "vtable" model, where each object includes a reference to a type descriptor, which in turn contains a table of function pointers giving the right function to invoke for each method for an object of that type. With double dispatch this becomes trickier. One of the simplest strategies, and an effective one in practice, is to have the initial vtable dispatch the method to a "stub" method, whose job is to perform another vtable dispatch on the second argument. Because most methods don't require double dispatch, this avoids imposing a penalty on ordinary single dispatch methods.

Multiple dispatch

But why stop at two arguments? Why not choose the method based on all the arguments? Although increasingly less useful, there are legitimate applications of methods that dispatch on 3, 4, or more arguments. This form of dispatch, called multiple dispatch, is directly supported by languages such as Dylan, Cecil, and CLOS (the Common Lisp Object System). Not only does this feature ease design of an initial system, but it also greatly increase the ability to extend it in multiple directions in an object-oriented way. Try experimenting with some of them and going through their tutorials; even in languages that don't directly support multiple dispatch, it helps to be able to recognise a situation where it would apply and simulate it.

Thanks for reading, everyone.