Deep Learning and Solving NP-Complete Problems

This recent article about AlphaGeometry being able to solve Olympiad geometry proof problems brings to my attention that there have been active research in recent years to combine deep learning with theorem proving. The novelty of AlphaGeometry is that it specializes in constructing geometry proofs, using synthetic training data guided by the theorem prover, which previous approaches have had difficulty with.

As a formal method like type checking of programming languages, theorem proving tries to overcome the ambiguity of mathematical proofs written in a natural language by formalizing proof terms using an algebraic language with rigorous semantics (e.g. Lean), so that proofs written in the language can be verified algorithmically for correctness. In the past, proofs had to be tediously written by humans, as computers had great difficulty constructing the proofs. Computers can do exhaustive proof search, but the search space blows up exponentially on the number of terms, which renders the exhaustive search intractable. However, once the proof is constructed, verification is quick. Problems that are exponentially intractable to solve but easy to verify are called NP-complete.

Traditional proof searching limits the search space to simple deductions, which can offload some tedium from manually writing the proofs. In recent years, machine learning is used to prune or prioritize the search space, making it possible to write longer proofs and solve larger problems. It is now possible for computers to construct complete mathematical proofs for a selected number of high school level math. The idea is a natural extension to what AlphaGo does to prune the game's search space, which I have written about before.

This is quite promising for a programming language that embeds theorem proving in its type system to express more elaborate safety properties, such as Applied Type System (ATS) designed by my advisor, which I have written about before. Programs written in ATS enjoy a high degree of safety: think of it as a super-charged Rust, even though ATS predates it. Now that generative AI can help people write Rust, a logical next step would be to create a generative AI for ATS.

I would be interested to see how machine learning can be applied to solve other NP-complete problems in general, such as bin-packing, job shop scheduling, or the traveling salesman problem. However, machine learning will face fierce competition here, as these problems are traditionally solved using approximation algorithms, which comes up with a close-enough solution using fast heuristics. Also, only a few of them have real-world applications. Admittedly, none are as eye catching as "AI is able to solve math problems" giving the impression of achieving general intelligence; in reality, the general intelligence is built into the theorem proving engine when it is used in conjunction with machine learning.

However, I think proof construction is a very prolific direction that machine learning is heading. Traditionally, machine learning models are a black box: even though it can "do stuff," it could not rationalize how to do it, let alone teaching a human or sharing the same insight with a different machine learning model. With proof construction, at least the proofs are human readable and have a universal meaning, which means we can actually learn something from the new insights gained by machine learning. It is still not general intelligence, and the symbolic insight to be gained might very well be a mere probabilistic coincidence, but it does allow machine learning to begin contributing to human knowledge.

Embedding a Binary Search Tree in an Array

Recently I started researching what would be the most efficient data structure to represent an IP address lookup table. An IP address table is useful for making routing decisions or access control, which checks whether a single IP address is a member of a range of netblocks. A netblock is represented by a base address and a prefix length, e.g. 192.168.42.0/24 represents a range of IP addresses from 192.168.42.0 through 192.168.42.255. Hash table is not suitable for this purpose because hashes are computed on exact matches. With the given example, this would require entering all 256 addresses in the range into the hash table, which is a significant blowup.

Hardware routers use Ternary Content Addressable Memory which can tell whether a specific value matches anything in the memory, after applying a bitmask. A software implementation might use a prefix tree (trie), which theoretically deduplicates common prefixes among the entries, resulting in a smaller table. However, any data structure employing pointers has the added overhead of the pointers (a machine word) and cache locality problems, which are not considered in theoretical papers.

The most compact data structure for the vast majority of lookup tables is going to be just an array of IP netblocks. For an unsorted array, lookup takes \( O(n) \) time which will be slow for a big table. If we presort the array and use binary search, the lookup takes \( O(\lg n) \) time. But binary search is a pathological case for cache locality because lookups tend to jump back and forth until it settles in to a target. Rather than a binary search which partitions candidates by two parts, it was shown that ternary search (three partition) improves performance, while quaternary search (four partition) improves little over binary search.

It had occurred to me that cache locality for binary search could be improved if we had simply rearranged the presorted array to become an embedding of binary search tree (BST), such that for node at index \( i \), its left and right children are at \( 2i + 1 \) and \( 2i + 2 \) respectively. The scheme is inspired by how binary heap is embedded in an array, but we are embedding a complete BST. The rationale is that the search always occurs left to right in the array, with the most accessed items located at the beginning of the array, which gets the most cache hits.

Here is an example showing how a presorted array with the number sequence \( 0 \dots 9 \) could be embedded.

Note that this embedding works for any presorted array. The values for \(i\) are the indices of the BST embedding, and the values for \(x\) are the indices in the presorted array. A presorted array can always be represented by a complete BST, so the array is dense. However, we can't embed a prefix tree this way because a prefix tree is sparse.

The implementation strategy is to presort the array, create a new array, and copy the items from the presorted array to the new array in the order of the BST embedding. The preprocessing only needs to be done once, and it results in the most compact representation of the lookup table. The assumption is that once the table is preprocessed, it won't need to be changed again.

I implemented the embedded BST in Go and compared its lookup performance with the builtin sort.Search() which performs binary search on a presorted array. The runs are repeated for increasing array sizes in the powers of 2 as well as stride width (both in terms of the number of machine words, not in bytes). The items to be searched are incremented by the stride width for each iteration; small stride approximates sequential lookup, whereas large stride approximates random lookup. Each of the data point is the average number after running some 2,000,000 to 100,000,000 iterations (the number of iterations varies due to the idiosyncrasies of the Go benchmark framework, not my fault; it tries to run each case for about 5 seconds).

The folly of this is that sequential lookup (smaller stride widths) may not cover a significant portion of a very large array size \(2^n\) where \(n \ge 20\), but the experiment is fair because (1) both algorithms are subject under the same condition, and (2) the array size above \(2^{20}\) already exceeds L2 cache size on the machine running the experiment, which is sufficient to demonstrate performance difference due to cache locality. The performance numbers are collected on a late 2013 model of MacBook Pro 13" whose Core i5 "Haswell" I5-4258U CPU runs at 2.4 GHz and has 32KB L1 cache and 3MB L2 cache. This machine has 8GB DDR3 ram. The program is compiled in 64-bits.


[popup figure in a new window]

You can use the drop-down menus to highlight the series for a specific algorithm (binary search, embedded BST) to see how they perform at various stride widths, or to highlight specific stride widths to compare how the two algorithms perform.

Here are some observations.

• At small array sizes \(2^n\) where \(n \le 7\), all lookups fit in the L1 cache. The two algorithms have identical performance at all stride widths.
• From array sizes \(2^7\) to \(2^{12}\), the L1 misses start to happen. Array size \(2^{12}\) occupies 32KiB, the size of L1 cache.
• Array sizes \(2^{12}\) through \(2^{18}\) still fit within the L2 cache, after which the performance starts skyrocketing at the largest stride width (which simulates random access). Even so, embedded BST held up much better than binary search.
• At small strides, embedded BST has a more predictable performance than binary search across all array sizes.
• Embedded BST performs at least as well as or better than binary search across all array sizes and strides.

Admittedly, embedded BST is optimized for lookup. After the initial preprocessing, we don't try to insert or remove entries from the table. If a modifiable table is desired, we can relax the density of the array so that the BST only needs to be approximately balanced. This allows us to use any classical balanced BST algorithms like red-black tree or AVL tree. It may be worthwhile to investigate the cache performance of array embedded versions of these algorithms as future work.

Some IP address lookup use cases require inserting and deleting individual addresses, such as for the purpose of intrusion detection. Since they only need to match individual addresses rather than netblocks, a hash table would suffice.

Otherwise, embedded BST is a great data structure for immutable IP address range lookup tables that is initialized once and optimized for lookup. It has the most compact representation and has much improved cache locality over binary search.

Fermat Numbers and Bit Mask

There is a little known application of Fermat numbers for generating bit masks, that I discovered back in 2010 while implementing "counting the number of one bits" using C++ template meta-programming. In the algorithm, bit masks such as 0x55555555, 0x33333333, 0x0f0f0f0f, etc. are used to mask out bits in specific locations. These bit masks can be generated by dividing the maximum integer 0xffffffff (or simply ~0u) against Fermat numbers 3, 5, 17, 257, 65537, etc.

Hex	Binary	Expression
0x5555...	0101010101010101...	~0u / 3
0x3333...	0011001100110011...	~0u / 5
0x0f0f...	0000111100001111...	~0u / 17
0x00ff...	0000000011111111...	~0u / 257
and so on...

It's not a coincidence that integers with all one bits set can be divided cleanly by a Fermat number, and this works with 8-bit, 16-bit, 32-bit, 64-bit, and beyond. To understand why, a \(k\)-bit all-one integer is \(2^k - 1\). In particular, \(k = 2^n\) in our case. We can factor \(2^{2^n} - 1\) using the formula \(a^2 - b^2 = (a + b)(a - b)\) we learned in grade school, here using the special case \(b = 1\), so \(a^2 - 1 = (a + 1)(a - 1)\).
\begin{align}
2^{2^n} - 1 & = (2^{2^{n - 1}} + 1)(2^{2^{n - 1}} - 1) \\
 & = (2^{2^{n - 1}} + 1)(2^{2^{n - 2}} + 1)(2^{2^{n - 2}} - 1) \\
 & = (2^{2^{n - 1}} + 1)(2^{2^{n - 2}} + 1)(2^{2^{n - 3}} + 1)(2^{2^{n - 3}} - 1) \\
 & = \dots
\end{align}
Conveniently, the last term \(2^{2^0} - 1\) is just \(1\), so we can write the product in terms of \(F_i = 2^{2^i} + 1\) for the \(i\)-th Fermat number,
\[ 2^{2^n} - 1 = \prod_{i=0\dots n-1}{2^{2^i} + 1} = \prod_{i=0\dots n-1}{F_i} \]
For example,

\(k\)	\(2^k - 1\)	Expression
2	3	\(3 \)
4	15	\(3 \times 5 \)
8	255	\(3 \times 5 \times 17 \)
16	65535	\(3 \times 5 \times 17 \times 257\)
32	4294967295	\(3 \times 5 \times 17 \times 257 \times 65537\)
and so on...

So how are Fermat numbers related to bit mask construction? It turns out that if you write 3, 5, 17, 257, 65537 in binary, you get 11, 101, 1 0001, 1 0000 0001, and 1 0000 0000 0000 0000 0000 0000 0000 0001 respectively. If you multiply 101 (5) and 1 0001 (17), you get alternating patterns of 0101 0101 (85). If you multiply this by 11 (3), you get an 8-bit all-one integer 1111 1111 (255). Conversely, you can take a pattern "out of circulation" by dividing the all-one integer by a Fermat number which is one of its factors.

Appendix: Counting the number of one bits using C++ template meta-programming

We start with the straightforward bit vector algorithm:

x = ((x >> 1) & 0x55555555) + (x & 0x55555555);
x = ((x >> 2) & 0x33333333) + (x & 0x33333333);
x = ((x >> 4) & 0x0f0f0f0f) + (x & 0x0f0f0f0f);
// etc.

which can be represented by the following for-loop:

for (int n = 1; n < sizeof(Tp) * 8; n *= 2)
  Tp mask = ~0 / ((1 << n) + 1));
  x = ((x >> n) & mask) + (x & mask);
}

In the template meta-programming version, the template __unroll_shift_mask_add essentially unrolls into the for-loop above, but in a manner that is agnostic to the actual integer type.

template<unsigned int n, typename Tp>
struct __unroll_shift_mask_add {
  static Tp compute(Tp x) {
    Tp y = __unroll_shift_mask_add<n / 2, Tp>::compute(x);
    Tp mask = ~Tp(0) / ((Tp(1) << n) + 1);
    return ((y >> n) & mask) + (y & mask);
  }
};

template<typename Tp>
struct __unroll_shift_mask_add<0u, Tp> {
  static Tp compute(Tp x) { return x; }
};

template<typename Tp>
Tp one_bits(Tp x) {
  return __unroll_shift_mask_add<sizeof(Tp) * 8 / 2, Tp>::compute(x);
}

How to Play Go Like a Computer Scientist

The board game of Go, or 圍棋 (weiqi), started receiving a lot of attention when AlphaGo beat the European champion Fan Hui in all of the five games. Now people are worried when AlphaGo is on its way to beat the world champion Lee Sedol, winning 3 out of 5 games so far. AlphaGo plays the game like no human as seen the game played before, inciting the awe of sadness and beauty. It is touted as an accomplishment of machine intelligence because traditionally the game of Go is very difficult for a machine to play well, whereas games like Chess is well-understood.

Board games are fascinating because a few number of rules allow an exponential number of possible game moves, some with winning outcomes, and some with losing. Traditionally, a naive game solver for any game tries all the possible positions of each alternate side's game play some number of steps ahead and chooses the first move that leads to the best gameplay after that number of steps, known as the min-max algorithm. Pruning heuristics is applied to reduce the search space so the computer does not waste time searching down "lost cause" paths, known as alpha-beta pruning. Even with pruning, the search space is still exponential to the number of positions on the board.

The keyword here is exponential. The machine can spend more time to compute, but for the amount of time spent \( t \), it can only search by \( \log{t} \) more steps. We can attempt to parallelize exponential search by enlisting the help of some number of machines, but it will still only allow us to advance the search depth by \( \log{n} \), i.e. the added intelligence is logarithmic to the number of machines.

Chess is well-understood. The number of game play combinations are vast, but we have found ways to prune the search space so that it is tractable with computer clusters that are practical to build. Computers can beat human playing Go on a smaller \( 9 \times 9 \) board with an upper bound of \( (9 \times 9)! \) or \( 10^{121} \) possible game plays, but on a professional championship level, the board size \( 19 \times 19 \) sports an exponential explosion of an upper bound of \( (19 \times 19)! \) or approximately \( 10^{768} \) possible game plays.

Even though Google might undertake the most ambitious of data center builds, they couldn't possibly shut down all the kitten videos on YouTube so AlphaGo could play against Lee Sedol with super-intelligence. Even so, they would achieve only logarithmic increase in intelligence. Building machines to brute force an exponential algorithm is a sure way to go bankrupt (NSA is an exception because the government prints the money).

I feel that a disclaimer is in order before I write any further. I don't have insider knowledge about how AlphaGo works. Everything I'm about to speculate comes from common sense and public sources.

AlphaGo algorithm reduces the search space using Monte-Carlo tree search, which means it randomly prunes some moves to reduce the branching factor, which allows the search to go deeper. The pruning is guided by deep neural network so it decreases the probability of pruning winning moves. Deep neural network works by building a hierarchy of signal processing nodes that amplifies certain signals after training.

When using DNN to recognize pictures, we can visualize how it works by looking at the images synthesized iteratively in order to boost certain signals.

In computer science, a pervasive strategy to design a super-efficient algorithm is by divide and conquer: basically, a big problem is solved by breaking down the problem into smaller sub-problems, until the problem becomes trivially solvable. Here are some examples:

• Sorting is taught to every first year computer science undergraduate. Given \( n \) items to sort, a naive sorting algorithm like bubble sort will take \( O(n^2) \) steps, whereas better ones using divide and conquer like quick sort or merge sort will take \( O(n \lg{n}) \) steps. If we are to sort 1000 items, this makes the difference between 1 million steps verses just 10000 steps.
• First year graduate student would likely learn about discrete Fourier transform. A naive DFT takes \( O(n^2) \) steps, whereas fast Fourier transform using divide and conquer takes \( O(n \lg{n}) \) steps.
• A naive matrix multiplication takes \( O(n^3) \) steps for an \( n \times n \) matrix, but Strassen's algorithm using divide and conquer takes \( O(n^{\lg{7}}) \) or approximately \( O(n^{2.8}) \) steps.
• The n-th Fibonacci number can be computed in \( O(\lg{n}) \) time using a divide and conquer method. The straightforward iterative method would take \( O(n) \) time, and a naive recursive method would take \( O(2^n) \) time. See Fibonacci Number Program for some examples in code.

For other algorithms, divide and conquer may not reduce the number of steps, but it still improves computation performance due to better cache locality. Divide and conquer is the basic idea of cache oblivious algorithms.

In other words, the real innovation of AlphaGo is not because machines achieved singularity, but because deep neural network approximated the divide and conquer method which happens to allow Go to be played efficiently.

The outcome of Go is scored by determining whether the pieces can survive locally, and the territorial occupation by the surviving pieces. This makes the winning strategy two-fold: on one hand, a player needs to ensure local survival of the pieces, but he has to take the territory into account. If a local battle is losing, it could be better for a player to strengthen his position on a different part of the board. This division of local and global thinking lends itself very well to using the divide and conquer method.

Here are some ideas how to divide. On the first level, we divide the board into subregions of \( 10 \times 10 \) each (figure left). We can also add overlapping subregions of \( 9 \times 9 \) or \( 9 \times 10 \) (figure right). This results in the first level division into \( 3 \times 3 \) subregions.

Left: Level 1

Right: Level 1 overlaps

Each subregion can be divided further into \( 5 \times 5 \) or \( 5 \times 4 \) grids with overlaps. In the following figure, only the division of the first subregion is shown, but it should be understood that the same division will take place for all first level subregions.

Left: Level 2

Right: Level 2 overlaps

A sketch of the divide and conquer algorithm would play Go by first choosing the region to play by valuation of the amount of territory it can likely occupy. Then within each region, it could brute force the positions that maximizes the survival of its own pieces.

Here is a very rough guesstimate of the number of steps needed to solve the game. These moves can be grouped into streaks of local battles taking place in \( 3 \times 3 \) subregions, and each streak is further divided to settle the \( 3 \times 3 \) sub-subregions which are about \( 5 \times 5 \) by grid. This reduces the amount of search space to be just:
\[ (5 \times 5)! \cdot (3 \times 3)! \cdot (3 \times 3)! \approx 10^{36} \]
This is a lot more computationally tractable. Of course, the caveat of my calculations is that I have never implemented a divide and conquer player for Go, and I don't know how well it works in practice. But even if we have to add a few levels of \( 3 \times 3 \) division, multiplying the search space by \( 9! \) still grows a lot slower than the exponential explosion of \( (9 \times 9)! \).
\[ (5 \times 5)! \cdot (9!)^9 \approx 10^{75} \]
Using this divide and conquer method, we can envision a game of Go played on a much bigger board, maybe like the 3-dimensional chess depicted in Star Trek.

The difference between the world champion Go player such as Lee Sedol and a modest Go player who happens to be a computer scientist, is that a computer scientist understands the power of divide and conquer. There is an immense beauty of divide and conquer that reaches far beyond game play. Nobody is worried when a computer doing quick sort can sort numbers faster than a human can with no sorting strategy, so why should we worry that a computer using divide and conquer can play a board game better than humans with only an ad-hoc strategy?

No, the computer is not about to achieve singularity. It is just a tool that helped us understand the game better.

Recursively Extensible Open Addressing Hash Table

There are two common ways to implement a hash table: open addressing and separate chaining. Both use a hash function to map keys to indices within an array for fast lookup. They differ in the conflict resolution strategy when multiple keys map to the same index. Upon conflict, open addressing searches for the next available index in the array through a probing policy such as linear probing (increment index by the number of probing attempts), quadratic probing (add index by the number of probing attempts squared), or by enumerating through different hash functions. Separate chaining avoids conflict by maintaining keys of the same index as a linked list.

Open Addressing

Separate Chaining

Both tables can be subject to O(n) lookup for a naive implementation when there are many collisions. A better one would dynamically resize the array depending on the load factor, but resizing is an expensive operation. For separate chaining, performance degradation can also be mitigated by using a balanced binary search tree instead of a linked list. Since the worst case access time is bound to O(lg n), it can retain some efficiency without resizing.

So far, this is text book material prior art.

In terms of cache efficiency, separate chaining has the disadvantage that list nodes may be scattered, so accessing the list can cause many cache misses. Open addressing is more cache friendly because the key and values are stored in a contiguous array. The question is: can we do away without resizing like separate chaining, but enjoy the cache efficiency of open addressing?

Yes we can. Here I introduce the recursively extensible open addressing hash table.

A recursively extensible open addressing hash table keeps the key/value in the table, but also reserves a pointer to an extension subtable per entry. To resolve conflict, it would first probe within the same table for an available slot but limited to k tries with a small k. If the search fails, then it would descend to the extension subtable at the original hashed index, creating the subtable only if necessary. It would use an orthogonal hash function to compute the index within this subtable, repeat the probing, and descend if necessary.

An orthogonal hash function is linearly independent (in the linear algebra sense) to the hash functions used in parent tables. An easy way to compute the orthogonal hash is to consider a universal hash function producing a 32-bit unsigned integer output. The first level hash table has 256 entries, and the index is based on the value in the 8 least significant bits. The second level tables also have 256 entries, with index drawn from the next 8 bits, and so on for the third and fourth levels.

As for the probing policy, we could use any policy like any open addressing hash table would, but linear probing has the best cache locality, and it is sufficient to achieve high load factor without compromising lookup performance. That is because the subtables all enjoy O(1) access as well.

Ideas for improving remote free for StreamFlow allocator

StreamFlow is a memory allocator that uses a simple idea to deal with producer-consumer threads, where an object is allocated from the producer thread, sent to the consumer, and subsequently freed. The object is "remote freed" back to the original heap using an atomic list insertion involving just a compare and swap. The remote free list can be vacated by the producer again using an atomic exchange. There is no need for complex lock-free wait-free algorithm.

But the remote-free lists themselves can become a performance bottleneck, as each remote free causes exclusive cache line to be migrated to the CPU performing the free. A pathological case is when any of the \(N\) threads could free to any other threads, causing cache line ping-pong and memory bus contention.

An idea is to cache an object to be remote-freed in a local free list first. For object sizes smaller than a cache line, reusing this local free list for allocation could cause false sharing, so this local free list is only for later remote free. Once this list reaches a certain length (say 1024), the objects in the list are sorted by the heap into batches, and then each batch is remote-freed in a single compare and swap.

Here is how the sorting would work. The sorter would allocate an array of 1024 pointers on stack, copy the object addresses to this array, and run in-place quicksort. Copying to the array might incur the most cache miss, but subsequent sorting would benefit by \(O(lgn)\) of cache hits. This idea came from What's the fastest algorithm for sorting a linked list? on StackExchange.

For object sizes greater than a cache line, false sharing occurs rarely, so the local free list could be used to satisfy allocation request first in the fashion of tcmalloc. It would be cleared using batched remote free if the list overflows.

Zero-Latency Computer Audio Effect Notes

When it comes to audio processing, it is believed that the human ear can detect at least 13ms of delay. In a computer audio setting, a large delay can be attributed to the underlying operating system. For example, it could be that the OS or driver does not handle interrupts (making input data available) fast enough, not scheduling threads to process the data soon enough, and making too many copies of audio data which delays the input as well as output. On Windows, the latency could be as high as 100ms. On Linux with preemptible and low-latency kernel, latency lower than 1ms is achievable. (These numbers are pulled from memory, and if any reader cares to point me to some latency measurement results—more recent the better—that would be very welcome).

But beyond the operating system, computer audio latency is also limited by the sampling rate as well. For example, when recording stereo audio at 44100Hz with 16-bit samples, filling a 1024 byte buffer takes 6ms. A buffer size for sub-millisecond latency would need to be 128 bytes, which is about the size of a memory cache line. This translates to about 32 samples. A very small number for audio effects processing. For example, computing Fourier Transform on a block size this small is not very interesting. This means that in order to achieve zero-latency computer audio effect, the algorithm has to be designed as a stream of samples.

Here I'm just noting two algorithms that are streamable: low-pass and high-pass filters. Both algorithms have a “discrete-time realization” (taken from Wikipedia) as follows:

• Low-pass filter: \[ y_i = \alpha x_i + (1 - \alpha) y_{i-1} \qquad \text{where} \qquad \alpha \triangleq \frac{\Delta_T}{RC + \Delta_T} \]
• High-pass filter: \[ y_i = \alpha y_{i-1} + \alpha (x_{i} - x_{i-1}) \qquad \text{where} \qquad \alpha \triangleq \frac{RC}{RC + \Delta_T} \]

In both cases, the cutoff frequency is \( f_c = \frac{1}{2\pi RC} \), and \( \Delta_T \) is the duration of time between samples (the reciprocal of the sampling rate).

The Tower Hotel Problem

Dynamic storage allocation problem, both as an online problem and an offline problem, is widely studied in literature, such as this article by Jordan Gergov. Here I present a similar problem which I call "The Tower Hotel Problem."

The Tower Hotel is peculiar. It has an infinite number of floors, but each floor has the same fixed number of rooms. The proprietor wanted to keep the number of floors open to the public minimal in order to save costs on maintenance; a closed floor incurs no cost. A floor is closed if none of the rooms are occupied. But as long as there is at least one occupied room, the floor must be open to public. Once a guest arrives and settles in a room, the proprietor must not ask the guest to move to another room.

In a zealous attempt to minimize cost, the proprietor requires all guests to make reservations before the year begins. The guest would indicate upon reservation the arrival and departure dates. Once the year begins, no more reservation can be made. The proprietor then optimizes the placement of guests on floors. This eccentric policy may be one reason the hotel is not very well known.

A number of years passed, and the proprietor passed away. His son knows nothing about hotel management, so the son interviewed three potential managers. In order to make the hotel more popular, the son required that the managers must not assume the reservations are made in advance. In fact, he went a step further and said that a guest could arrive at any time and leave at any time.

Each manager had a different strategy. The first manager would place the guest at random on any of the opened floors with vacancy (and open up a new floor if none of the existing floors are available). The second manager would place the guest on the floor with the most number of occupied rooms (but still has vacancy). The third manager would place the guest on the lowest numbered floor with vacancy. The third manager has an assistant who, due to her desire to be promoted, secretly proposes to the proprietor's son a slightly different strategy than her manager: she would locate the lowest numbered floors with vacancy, fill it entirely first, then find again the lowest numbered floors with vacancy, which may be lower than the previously chosen floor because some guests below may have checked out.

How does the proprietor's son determine which manager has the most cost-effective strategy? How do these managers and the third manager's assistant compare with what the proprietor could do, when the proprietor had full knowledge of all the reservations in the whole year?

---

This problem describes segregate fits dynamic storage allocation with virtual memory where all objects are of the same size, and a memory page could fit a fixed number of these objects. A virtual memory page has to be backed by a physical memory page if there is at least one object in the page. Once all objects in the page are freed, the physical memory page backing can be removed. This presents a different sort of fragmentation problem where certain long-lived objects can pin a memory page even though the page is mostly unused.

Evaluating RC4 as pseudo-random number generator

Update (September 3, 2012): I made some mistakes in my own implementation of RC4. An analysis of the mistake is documented in Mistake in my RC4 implementation. After fixing the mistake, I do not see the skewed distribution shown here. Please be advised that none of the blog posts in this blog are peer reviewed, and mistakes happen.

On the quest to look for a pseudo-random number generator that is easy to implement, I gave up on Mersenne Twister but took Wikipedia's advice that a stream cipher can be used for generating pseudo-random sequences. But it also mentions that RC4 is not a particularly uniform stream cipher which makes statistical attacks possible. Because of the simplicity of RC4, I decided to implement it and see really how bad it is. Since RC4 outputs a stream of bytes, I took the histogram of the byte value. The initialization key is the integer value of a stack-allocated object pointer, repeated for 64 bytes, and it's run on a 64-bit machine. I'm not sure if the fact that the key is highly regular has anything to do with the following non-uniform result that I observed.

byte	count	byte	count	byte	count	byte	count	byte	count	byte	count	byte	count	byte	count
0	5499938	32	50376	64	81740	96	253697	128	50004	160	144433	192	147856	224	147031
1	134854	33	50103	65	50386	97	115718	129	50024	161	49930	193	149593	225	50260
2	50352	34	50179	66	954060	98	50362	130	50185	162	50553	194	134910	226	78249
3	115332	35	50006	67	50091	99	50170	131	50543	163	50289	195	50516	227	49853
4	50214	36	50065	68	50256	100	50052	132	50288	164	49950	196	50532	228	50303
5	50089	37	50589	69	50173	101	50118	133	143783	165	49753	197	50182	229	50135
6	207241	38	49946	70	50416	102	50395	134	50447	166	50270	198	50145	230	134634
7	68820	39	50147	71	49985	103	49774	135	50252	167	50502	199	49368	231	49987
8	50212	40	103906	72	49950	104	50501	136	49803	168	50038	200	50032	232	50486
9	695154	41	49905	73	50137	105	50389	137	50085	169	50619	201	66006	233	49849
10	50439	42	49865	74	50472	106	50284	138	49508	170	50045	202	50336	234	50083
11	49682	43	948555	75	50141	107	50230	139	50294	171	50719	203	50300	235	50267
12	50378	44	49965	76	140698	108	50119	140	49772	172	50399	204	49963	236	50117
13	50466	45	50171	77	49814	109	50012	141	50385	173	50031	205	50072	237	49945
14	49591	46	50290	78	50219	110	90366	142	50177	174	50126	206	50585	238	50157
15	49776	47	50340	79	49916	111	49974	143	50224	175	50513	207	50350	239	50234
16	287052	48	205931	80	151232	112	131863	144	50079	176	50026	208	50226	240	131317
17	50022	49	50076	81	50134	113	49801	145	50487	177	50688	209	81920	241	50161
18	118991	50	50125	82	50108	114	50370	146	66545	178	78182	210	49998	242	50092
19	50209	51	50706	83	50292	115	50297	147	50236	179	50411	211	50458	243	50180
20	50170	52	50025	84	221044	116	116139	148	247793	180	215784	212	50185	244	50124
21	50371	53	49717	85	50259	117	50274	149	50327	181	758974	213	49808	245	90855
22	50365	54	49979	86	50302	118	50523	150	50228	182	50272	214	50685	246	50228
23	50467	55	50145	87	49797	119	105677	151	50108	183	50390	215	113514	247	50316
24	100861	56	78827	88	113581	120	138358	152	57005	184	53395	216	147492	248	50546
25	50158	57	49803	89	50264	121	50209	153	50095	185	50334	217	50396	249	50078
26	50483	58	50279	90	50176	122	49674	154	50292	186	50529	218	49915	250	50111
27	50211	59	49997	91	50337	123	50396	155	50453	187	49836	219	50170	251	50524
28	49859	60	50507	92	50313	124	50280	156	50059	188	50229	220	248053	252	49981
29	50580	61	50869	93	50417	125	49877	157	50355	189	50313	221	49937	253	50533
30	50142	62	50358	94	50292	126	50240	158	49821	190	50204	222	50139	254	49975
31	50265	63	50165	95	50352	127	244313	159	50410	191	50529	223	202938	255	246009

The result is pretty staggering non-uniform and heavily biased to 0. I quickly eyeballed the figures and labeled the counter with colors. Red means a strong bias. Yellow is a slight bias. Green is the ballpark figure for unbiased result, but they are still above a typical count in this table. The coloring is quite subjective, but I'm only using it to highlight the non-uniform bias.

Now, this observation also indicates that RC4 is a very poor stream cipher. The way stream cipher works is essentially that a stream of pseudo-random bytes are generated and XOR-ed byte-wise with the plaintext in order to obtain ciphertext. The prevalence of 0 values means that most of the plaintext will be identical to the ciphertext because 0 is an XOR identity. The fact that this bad encryption is possible with some key choices should worry a lot of people who are still using RC4.

By a stroke of luck, I tried a slightly modified RC4 that actually gives me a pretty uniform distribution. The difference (from the pseudo-code listing of the RC4 wikipedia article) is highlighted in grey below.

i := 0
j := 0
while GeneratingOutput:
    i := (i + 1) mod 256
    j := (j + S[i]) mod 256
    swap values of S[i] and S[j]
    K := S[(S[i] + S[j]) mod 256]
    output K ^ j  // was just K before
endwhile

This modification gives me the following histogram.

byte	count	byte	count	byte	count	byte	count	byte	count	byte	count	byte	count	byte	count
0	99723	32	100388	64	99739	96	99866	128	100026	160	100026	192	99831	224	99581
1	99941	33	99908	65	100363	97	100328	129	99741	161	100338	193	99792	225	99776
2	99968	34	100010	66	99589	98	100231	130	100649	162	100084	194	100218	226	99839
3	99896	35	99924	67	100465	99	100002	131	100093	163	99815	195	100478	227	100137
4	100439	36	99677	68	100004	100	99960	132	100165	164	99756	196	100064	228	100114
5	100081	37	99362	69	99664	101	99872	133	99887	165	100059	197	99523	229	100590
6	100570	38	99924	70	99648	102	99945	134	99459	166	99854	198	100170	230	100069
7	99812	39	99706	71	100066	103	99598	135	99532	167	99499	199	99956	231	99892
8	99969	40	100484	72	100074	104	99506	136	99831	168	100291	200	100164	232	99562
9	99840	41	99756	73	99128	105	100178	137	99646	169	100348	201	100121	233	99412
10	99793	42	100133	74	100019	106	100488	138	99882	170	99661	202	99864	234	99878
11	100208	43	100059	75	99717	107	100141	139	100374	171	100116	203	99643	235	100394
12	100157	44	99772	76	99846	108	100085	140	100180	172	99632	204	99928	236	100123
13	100098	45	100232	77	100017	109	100654	141	99955	173	99993	205	99858	237	100188
14	100017	46	100135	78	99566	110	99104	142	99640	174	99256	206	100243	238	99918
15	100259	47	100403	79	99459	111	99887	143	100255	175	99818	207	99612	239	99840
16	100139	48	99888	80	100828	112	99999	144	99820	176	100378	208	99663	240	99622
17	100093	49	100173	81	99633	113	100480	145	99804	177	99826	209	100564	241	99916
18	100257	50	99987	82	100502	114	100253	146	100554	178	99998	210	100183	242	100319
19	99835	51	100206	83	99776	115	99674	147	100726	179	99909	211	100548	243	99803
20	100010	52	99879	84	100188	116	99657	148	100255	180	99687	212	99937	244	100122
21	100424	53	100257	85	100197	117	100375	149	99866	181	99741	213	100669	245	99972
22	100062	54	99376	86	100123	118	100396	150	99623	182	100523	214	100195	246	100034
23	99729	55	100404	87	100620	119	100275	151	99594	183	99879	215	99921	247	100092
24	100076	56	99792	88	100383	120	99712	152	100641	184	99514	216	100226	248	100600
25	100107	57	100214	89	100079	121	99682	153	99723	185	99764	217	100120	249	99574
26	99785	58	100284	90	99835	122	100090	154	99805	186	99752	218	99850	250	100211
27	99839	59	99459	91	100332	123	100185	155	99557	187	100176	219	100167	251	100213
28	99700	60	99795	92	99966	124	100599	156	99709	188	100422	220	100028	252	100368
29	99776	61	99663	93	99596	125	99479	157	100266	189	100427	221	100172	253	99842
30	100319	62	99734	94	99962	126	99874	158	100068	190	100238	222	100058	254	99921
31	99537	63	100025	95	100043	127	100382	159	99910	191	99824	223	100697	255	99923

This is quite amazingly uniform, at least appears good enough for my purpose. It remains to be seen whether two or more consecutive bytes of output will still retain the same uniformity.

I think it is quite possible that this slight modification could be immune to the existing statistical attacks on RC4.

A test on gcc's ability to optimize

I have an implementation of some singly linked data structure which I wish to adapt for owner reference. An owner reference is like a "parent" pointer in the case of a binary search tree, or the "previous" pointer in the case of a doubly linked list, i.e. a back pointer. However, the owner reference does not point to a node.

To understand what the owner reference is, I've written about tracking object ownership before, where we use a linear pointer to keep track of who is the owner of that object. Each object would have exactly one pointer that points to the object, and whoever has that pointer is the owner. The owner is responsible for disposing the object.

The principle problem with linear pointer is that it cannot be used to implement doubly linked data structure by definition, since doubly linked data structure implies that a node could have two incoming pointers. In a situation like that, only one pointer can be the owner. The other pointer must be a loan.

You can implement a doubly linked data structure with a linear pointer and a loan pointer, but we have to realize that the chief reason to have doubly linked data structure is so we can easily remove a node, knowing only its address, in O(1) time without needing to look up where it is in the data structure.

The owner reference in a node is a device to tell the address of the linear pointer that is the node's owner. With owner reference, it is possible to borrow an object and then steal ownership of it from the previous owner.

Using owner reference is preferred in case of binary tree node that requires in-place removal. If we use a parent node pointer, we would need additional logic to figure out whether the node we want to remove is the left or the right child of its parent. With owner reference, we don't care; the parent could even be the root pointer. The owner reference leads us right to the correct pointer that we need to modify.

However, the owner reference is additional book-keeping that needs to be updated whenever we use another routine to change the structure, such as list reversal. In the case of binary search tree, I would need to quickly augment a splay tree I implemented before with owner references. Obviously, owner reference changes whenever a linear pointer move occurs, and in C++ I could fortify the assignment operator of my linear pointer to take care of that for me. This would allow the same code that previously dealt only with singly linked data structure to become owner reference aware. However, the code could be doing extra work.

To illustrate, suppose we have a move function that supports owner reference update when doing assignment.

class node {
 public:
  node *next() const throw() { return next_; }
  node *& next() throw() { return next_; }

  node **owner() const throw() { return owner_; }
  node **& owner() throw() { return owner_; }

 private:
  node *next_;
  node **owner_;
};

void move(node *& to, node *& from) {
  to = from;
  to->owner() = &to;
  from = 0;
}

The move assignment updates the pointers linearly and keeps the owner reference updated. Once you move ownership from one pointer to another, the original pointer would lose the pointer to the node, and that's why it's filled with a NULL pointer, or 0. The owner would be the new "to" pointer.

Now consider the swap function.

void swap_naive(node *& a, node *& b) {
  node *tmp;
  move(tmp, a);
  move(a, b);
  move(b, tmp);
}

Notice that we can achieve the same effect with much less work.

void swap_fast(node *& a, node *& b) {
  node *tmp = a;
  a = b;
  b = tmp;
  a->owner() = &a;
  b->owner() = &b;
}

The question I ask in the subject of this post is, if we write swap_naive, is GCC smart enough to emit code in the fashion of swap_fast? If it can, then I could merrily overload the assignment operator of an augmented linear_ptr<> class to update owner references. However, this means I will not be able to write swap_fast() directly. If not, that means I'll need to manually optimize away the redundant owner reference updates.

I compiled opt.cc with g++ -S -O3 -fomit-frame-pointer with different gcc versions, and here is the result for the number of essential instructions (not counting boilerplate context saving and return) for either function.

Compiler (output)	swap_naive()	swap_fast()
gcc-4.4.1 -march=i686	11 (9)	9 (8)
gcc-4.4.1 x86_64	7	7 (6)
apple gcc-4.2 -arch i686	15	10
apple gcc-4.2 -arch x86_64	11	8
apple gcc-4.2 -arch ppc	12	8
apple gcc-4.2 -arch ppc64	12	8

The reason I decided to do x86_64 and ppc/ppc64 is because these architectures have more register count than i386, and therefore is able to retain more data in the register rather than needing to spill data to the stack frame.

It is definitely surprising how gcc-4.4.1 on x64_64 can achieve the same number of instructions for both swap_naive() and swap_fast(), although you can see from the assembly listing that the naive version does assigns zero to the next_ pointer. And closer inspection in the assembly of swap_fast() indicates that the movq (%rdi), %rdx instruction is redundant because we had movq %rdx, (%rdi) before. The compiler could have saved us one more instruction, but the optimization has a bug.

Regarding optimization bug, a closer inspection on gcc-4.4.1 on i686 assembly code reveals that the subl  $16, %esp and addl $16, %esp instructions are redundant. This means, with proper optimization, we could bring instruction count for swap_naive() down to 9 as well. Again, in swap_fast(), the movl (%eax), %ebx instruction is redundant, which brings the instruction count down to 8.

The compiler optimizer is indeed impressive, but not without flaws.

After this study, I determined that it is better to not blanket overload the assignment operator. I would modify singly linked routines to use move(), which is identical to simple assignment in singly linked data structures. However, for augmented linked structure, it will be possible to manually optimize some of these routines further by not using move(). For example, I could leverage a singly linked list sort() to sort an augmented list with owner reference, but I only relink the owner reference after sorting is done. This would not be possible had assignment operator been overloaded.

Finding n Consecutive Bits of Ones in a Bit Vector in O(lg n) Time

We want to find out whether a bit vector X of length m has n consecutive bits of ones, as well as all the positions where we can find it.

For example, the following bit vector:

(msb) 11111111 01111111 00111111 00011111 (lsb)

Has 5 consecutive ones starting at bit 0, 6 consecutive ones at bit 8, 7 consecutive ones at bit 16, and 8 consecutive ones at bit 24. In addition, we count 4 consecutive ones at bit 1, 3 consecutive ones at bit 2, and so on. The overlaps are also considered. Therefore, if we were to find all positions with 6 consecutive ones, we could generate a positional bit mask like this:

(msb) 00000111 00000011 00000001 00000000 (lsb)

The result indicates that we can find 6 consecutive ones in the original bit vector if we start at bits 8, 16, 17, 24, 25, 26.

Of course, if we can only access one bit at a time, the best running time we can do is O(mn). Most computers provide bitwise operators to compute logical AND and OR on the matching bits of two bit vectors simultaneously, as well as bit shift operators. This, as we show, will dramatically reduce our running time.

Our first attempt is to simplify the problem of finding two consecutive ones, that is, bit i should be set if both bit i and bit i + 1 are ones. This is easy. We shift X to the right by one bit, and AND it against the original unshifted X. In C, this is written as X & (X >> 1). We can generalize the first attempt to find n consecutive ones by doing X & (X >> 1) & (X >> 2) & ... & (X >> (n - 1)). This naive algorithm takes O(n) time.

The part that gets really interesting is when we observe that 4 consecutive ones can be detected by 2 consecutive ones followed by another 2 consecutive ones. In the original example,

(msb) 11111111 01111111 00111111 00011111 (lsb) [n = 1] Y0 = X
(msb) 01111111 00111111 00011111 00001111 (lsb) [n = 2] Y1 = Y0 & (Y0 >> 1)
(msb) 00011111 00001111 00000111 00000011 (lsb) [n = 4] Y2 = Y1 & (Y1 >> 2)
(msb) 00000001 00000000 00000000 00000000 (lsb) [n = 8] Y3 = Y3 & (Y2 >> 4)

The reason why this works is that in computing Y2, we use the fact that Y1 has gaps of 2 consecutive zeroes to shorten our span of ones, and in Y3 uses the fact that Y2 has gaps of 4 consecutive zeroes, and so on. This is another type of SWAR algorithm.

We will keep doing this until we reach a power of 2 that is <= n. Then the last step we just apply a shift to truncate off the remainder of the bits.

In the original example, suppose we want to find n = 7 consecutive ones.

(msb) 11111111 01111111 00111111 00011111 (lsb) [n = 1] Y0 = X
(msb) 01111111 00111111 00011111 00001111 (lsb) [n = 2] Y1 = Y0 & (Y0 >> 1)
(msb) 00011111 00001111 00000111 00000011 (lsb) [n = 4] Y2 = Y1 & (Y1 >> 2)
(msb) 00000011 00000001 00000000 00000000 (lsb) [n = 7] R = Y2 & (Y2 >> 3)

The last step R, instead of shifting by 4, we only need to shift by 3.

The final algorithm is this:

typedef unsigned int bitvec_t;

bitvec_t consecutive_mask(bitvec_t x, unsigned int n) {
  unsigned int stop = prev_pow2(n);
  for (unsigned int i = 1; i < stop; i <<= 1)
    x &= (x >> i);
  if (stop < n)
    x &= (x >> (n - stop));
  return x;
}

The for-loop runs in O(lg n) time. We still need to show that prev_pow2(), which computes the power of two that is <= n, also runs efficiently. Again, this can be done using a SWAR algorithm.

unsigned int prev_pow2(unsigned int n) {
  return (fold_bits(n) >> 1) + 1;
}

where fold_bits() is the SWAR algorithm defined as follows:

unsigned int fold_bits(unsigned int x) {
  x |= x >> 1;
  x |= x >> 2;   /* last for 4-bit */
  x |= x >> 4;   /* last for 8-bit */
  x |= x >> 8;   /* last for 16-bit */
  x |= x >> 16;  /* last for 32-bit */
  x |= x >> 32;  /* last for 64-bit */
  return x;
}

The remark "last for k-bit" indicates that if x only has k bits, then that line is the last operation for the fold_bits algorithm. Since our n only has 5 or 6 bits, we only need to do up to "last for 8-bit." It turns out that this part is also O(lg n).

Here we presented the problem of finding the positions of all consecutive bits of ones in a bit vector, and an algorithm to do it in O(lg n) + O(lg n) = O(lg n) time. An application for this algorithm is dynamic storage allocation. Suppose I have each bit in the bit vector to indicate the availability of a storage cell, then this algorithm allows me to find places where I can reserve n consecutive storage cells in O(lg n) time.

Are fast integer min and max operators portable?

Sometimes when writing high performance code, hackers have come to use hard coded assembly to compute integer min(x, y) and max(x, y) in a way described in The Aggregate MAGIC. The formula for integer min is x + (((y - x) >> (WORDBITS - 1)) & (y - x)). Of course, the common expression (y - x) can be computed just once. The constant (WORDBITS - 1) can be computed at compile time. This expression requires 4 instructions. This code is reasonably portable, as most computer architectures nowadays are two's complement, and they support addition, subtraction, and sign-preserving bit shifts.

Compared with the usual min(x, y) defined as:

int min(int x, int y) {
  int result;
  if (x < y)
    result = x;
  else
    result = y;
  return result
}

which is essentially what the expression (x < y)? x : y gets compiled to at the assembly level, the idea is that the "magic" version does not use branching, so it will work best with processors with pipelined instruction execution. The branching would cause the pipeline to be flushed. A naive processor with no branch prediction and speculative execution will suffer a heavy penalty.

On modern CPU, a lot of transistor die is dedicated to branch prediction and speculative execution. One chief reason is that it alleviates the burden on the programmer, so the programmer can continue to write simple code. It also alleviates the burden on compiler writer, so that the compiler does not need the heavy machinery to recognize pattern in code like (x < y)? x : y and substitute it for the branch-free version. However, the transistors doing branch prediction and speculative execution draw a lot of electric power. They also produce a lot of heat, requiring more power in a data center for heat management. It is becoming unfashionable in a "green conscious" economy.

Here comes the ARM processor that powers mobile phones and embedded devices. The reason it is used for low-power applications is because of its energy efficiency. Words are that ARM will enter the data center market in 2011 allowing data centers to reduce its energy use. In terms of technical merits as a server, ARM began to support atomic compare and swap instruction, a key ingredient in coordinating multi-threaded programs, in ARMv6 instruction set. In terms of hardware implementation, ARM11 MPCORE and ARM Cortex-A9 are multi-core processors.

One design choice that leads to the low-power requirement of ARM processor is its instruction set. Instead of having branch prediction and speculative execution, each instruction is conditionally executed. It allows a branch-free rendition of the min and max operator essentially like this:

int min(int x, int y) {
  bool cond = (x <= y);
  int result;
  cond && (result = x);
  !cond && (result = y);
  return result;
}

Hence, the expression (x < y)? x : y can be compiled to 3 instructions.

This presents an interesting dilemma to the programmer.

• If he were to program the "magic" version, his code will run faster on all pipelined architectures, but it may run slightly slower on ARM.
• If he were to use the standard (x < y)? x : y expression, the compiler may or may not optimize this code. He may not get consistent performance when using a different compiler.

Furthermore, it presents an interesting challenge to compiler writers. When writing a multi-target compiler, the source code is often compiled to a platform-neutral intermediate code by a front-end, and then from the intermediate code to the target machine code by the back-end. The intermediate code may have scrambled the source code enough so that it is hard for the back-end to recognize simple patterns like (x < y)? x : y and perform pin-hole optimization.

Finally, the programmer who wants the performance guarantee would probably just write in-line assembly for each target machine architecture. This is just another fine example where, if you want to maximize the ability of hardware, you still need to write code for that particular architecture.

Ordered Intervals

Although we have efficient data structures to express sets, such as balanced binary search tree and sparse hash table, some sets are most efficiently expressed as an interval. An example is the character set, which is typically used for regular expression parsing, e.g. the set [0-9A-Fa-f] denotes characters for hexadecimal number digits. Historically, ASCII only has 127 characters, so a direct-index table would suffice; the table could perhaps be bitmapped, where each bit represents one code point, for a total of 16 bytes in the table. The state transition table for finite state automata representing a regular expression will only use a couple of hundred bytes. However, supporting regular expression parsing for Unicode is more challenging. Unicode currently allocates 20-bits per character (Unicode 5.0), and a typical set could easily contain several dozen code points, so a direct-mapped table certainly will not work well (one bitmapped row in the state transition table would be 128KB). A plain old list of character interval would work much better.

For the purpose of discussion, let's forget about character sets.

Assume we have a set expressed by a set of intervals, so we determine if a value is a member of the set by examine if the value lies inside any of these intervals. Set union is straightforwardly defined by union of the set of intervals. We'll focus on set membership for now and ignore set intersection. If intervals are arranged in a list, searching still takes linear time. In order to allow binary searching in logarithm time, we need a way to order these intervals.

Ordering is tricky if intervals overlap. Fortunately, overlapping intervals can be trivially combined, e.g. [5, 10] ∪ [7, 12] = [5, 12]. However, there is still a corner case. How can we tell if the interval [5, 9] ∪ [10, 14] equals to [5, 14] or not? It depends. We can combine them if the intervals range over natural numbers; not if the intervals range over real numbers. However, we shouldn't make any assumption about the continuity of elements. In order to avoid corner case like this, we make one side of the interval an open interval. It is customary in computer science to use intervals like [a, b), where a is inclusive, and b is exclusive, because for-loops are usually written this way.

We'll have a rule saying [a, b) ∪ [c, d) = [a, d) if b ≥ c. In particular, [a, b) ∪ [b, c) = [a, c).

The rest of interval tree implementation can be seen in Introduction of Algorithms, pages 311-315.

Splay trees with duplicate items

I've implemented a splay tree, and I'd like it to store duplicate items, which are key-value mappings, as follows: if a key to be inserted already exists, then the new binding shadows the old one. If a key is to be removed, then the old binding is revealed. In other words, the inserting and removing the same key reveals the value in LIFO, or stack, order.
The most obvious solution is to have each node be a pointer to a singly linked list, with the cons/uncons operation trivially obeying stack order. But we could embed the singly linked list in a simple binary search tree if we relaxed the criteria a bit, for example making the left child <= parent as opposed to left child < parent. This effectively makes duplicate nodes a chain of singly linked list down the left path. The idea is that we always reach the parent first, which is the most recently inserted item. And if the parent is removed, it reveals the second next item of the same key.
The complication with Splay tree is that it could rotate the tree, which breaks the chain of duplicate nodes. Sometimes it is as bad as splitting the chain to left and right children to a node != the key of the chain, invalidating the binary search tree invariant. I figured out that, in the zig-zig and zig-zag cases, if the grand-parent node has the same key as the parent node, then it is possible to come up with special case rotation that preserve the chain. However, the terminal case actually needs to descend down the whole chain of duplicate notes. At this moment I gave up.
Also, another problem with embedding singly linked list in a tree is that it will increase the depth of the tree, hurting the overall running time. In this case, it seems that the simple approach, making each node a pointer to a singly linked list, is better. It only took me 5 minutes to change the implementation this way and it worked.
Update (9/15): Wikipedia claims that identical keys in a splay tree are preserved in its order like stable sort. It suffices to find the leftmost or rightmost node of the identical key. I haven't tried this.
Update (9/16): Here is the Wikipedia edit that added the aforementioned claim about stable sort. However, I don't think retrieving duplicate item is as simple as finding the leftmost or rightmost node. Even if insertion is stable, the duplicate nodes may not be continuous (e.g. the balanced binary search tree for 1, 2, 3, 3, 3, 4, 5, where the 3's are separated by nodes 2 and 4). Looking for the "largest" duplicate complicates the search by a great deal because you will have to probe either the minimal of the right child, or the maximal of the left child, and may have to back-track. Curiously, Splaysort did not mention how they deal with duplicates.

Manifesto of a Very Hard To Find Memory Leak (To Be)

I'm implementing a hash table that will become the basis of a free list memory allocator. The hash table uses singly linked list for buckets. In the list implementation, there is a RemoveAll function written like this:

void RemoveAll(Predicate& p, Consume& c)
{
Node **result;
for (result = &first_; *result != NULL; result = &(*result)->next_) {
if (p(**result)) {
Node *node = *result;
*result = node->next_;
node->next_ = NULL;
c(node);

if (*result == NULL)
break;
}
}
}

Admittedly, the code is copied from a similar function called Remove(), which removes only the first node that satisfies the predicate. However, when it is modified to remove all nodes, an unintended side effect is, if a node is to be removed, it always skips the next node. When I wrote a unit test, I inserted nodes for 0 through 16, and removed all the even nodes, so this problem was not detected. However, when I used this function with a predicate that always returns true, I noticed that the list is not cleared properly.
Now, this function is used in the hash table to clear all the buckets. Each item in the bucket would be a free list for a certain size of objects. When a free list memory allocator uses this hash table, it will see memory leak on a certain size of objects, and the size would appear to be at random.
Fortunately, when a list is destructed, I have it do nothing but to assert that the list is empty. This helped detecting the problem early. Unit test for the hash table to clear all entries helped bringing this problem to the surface.

Hash table using a variety of things for bucket

A hash table stores key-value mapping of data in an array of buckets. The bucket number for a particular key-value pair is decided by a hash of the key modulo the number of buckets. The idea is that if the average number of key-value pairs in buckets is low overall, then we can achieve constant lookup time. One way to organize the bucket is by using a singly linked list. This has the side effect that, if the buckets are overloaded, then lookup time approaches the time of the singly linked list, which is O(n).

An obvious improvement is to make the bucket a binary search tree, perhaps a splay tree (so that item being looked up more often is quicker to find). This cuts down the worst case time of overloading to O(logn), while achieving average case lookup time of O(1).

We can nest hash tables, so a bucket can be implemented by another hash table. In this case, if the load of buckets are small, then we can let two or more buckets initially share a hash table. As buckets grow in size, they split. The size of nested hash tables must be relatively prime, or the inner tables would have high collision rate and won't be effective. The hashing scheme can be organized in the fashion of radix sort where the outer table use one hash function on the most significant bits, and the inner table use another hash function on the least significant bits.