5 \chapter{Fine Details of Computation}
8 \section{Models and machines}
10 Traditionally, computer scientists use a~variety of computational models
11 as a~formalism in which their algorithms are stated. If we were studying
12 NP-completeness, we could safely assume that all the models are equivalent,
13 possibly up to polynomial slowdown which is negligible. In our case, the
14 differences between good and not-so-good algorithms are on a~much smaller
15 scale. In this chapter, we will replace the usual ``yardsticks'' by a~micrometer,
16 state our computation models carefully and develop a repertoire of basic
17 data structures taking advantage of the fine details of the models.
19 We would like to keep the formalism close enough to the reality of the contemporary
20 computers. This rules out Turing machines and similar sequentially addressed
21 models, but even the remaining models are subtly different from each other. For example, some of them
22 allow indexing of arrays in constant time, while on the others,
23 arrays have to be emulated with pointer structures, requiring $\Omega(\log n)$
24 time to access a~single element of an~$n$-element array. It is hard to say which
25 way is superior --- while most ``real'' computers have instructions for constant-time
26 indexing, it seems to be physically impossible to fulfil this promise regardless of
27 the size of memory. Indeed, at the level of logical gates inside the computer,
28 the depth of the actual indexing circuits is logarithmic.
30 In recent decades, most researchers in the area of combinatorial algorithms
31 have been considering two computational models: the Random Access Machine and the Pointer
32 Machine. The former is closer to the programmer's view of a~real computer,
33 the latter is slightly more restricted and ``asymptotically safe.''
34 We will follow this practice and study our algorithms in both models.
37 The \df{Random Access Machine (RAM)} is not a~single coherent model, but rather a~family
38 of closely related machines, sharing the following properties.
39 (See Cook and Reckhow \cite{cook:ram} for one of the usual formal definitions
40 and Hagerup \cite{hagerup:wordram} for a~thorough description of the differences
41 between the RAM variants.)
43 The \df{memory} of the machine is represented by an~array of \df{memory cells}
44 addressed by non-negative integers, each of them containing a~single non-negative integer.
45 The \df{program} is a~finite sequence of \df{instructions} of two basic kinds: calculation
46 instructions and control instructions.
48 \df{Calculation instructions} have two source arguments and one destination
49 argument, each \df{argument} being either an~immediate constant (not available
50 as destination), a~directly addressed memory cell (specified by its number)
51 or an~indirectly addressed memory cell (its address is stored in a~directly
52 addressed memory cell).
54 \df{Control instructions} include branches (to a~specific instruction in
55 the program), conditional branches (e.g., jump if two arguments specified as
56 in the calculation instructions are equal) and an~instruction to halt the program.
58 At the beginning of the computation, the memory contains the input data
59 in specified cells and arbitrary values in all other cells.
60 Then the program is executed one instruction at a~time. When it halts,
61 specified memory cells are interpreted as the program's output.
64 In the description of the RAM family, we have omitted several details
65 on~purpose, because different members of the family define them differently.
66 These are: the size of the available integers, the time complexity of a~single
67 instruction, the space complexity assigned to a~single memory cell and the set
68 of operations available in calculation instructions.
70 If we impose no limits on the magnitude of the numbers and we assume that
71 arithmetic and logical operations work on them in constant time, we get
72 a~very powerful parallel computer --- we can emulate an~exponential number
73 of parallel processors using arithmetics and suddenly almost everything can be
74 computed in constant time, modulo encoding and decoding of input and output.
75 Such models are unrealistic and there are two basic possibilities how to
79 \:Keep unbounded numbers, but increase costs of instructions: each instruction
80 consumes time proportional to the number of bits of the numbers it processes,
81 including memory addresses. Similarly, space usage is measured in bits,
82 counting not only the values, but also the addresses of the respective memory
84 \:Place a~limit on the size of the numbers ---define the \df{word size~$W$,}
85 the number of bits available in the memory cells--- and keep the cost of
86 instructions and memory cells constant. The word size must not be constant,
87 since we can address only~$2^W$ cells of memory. If the input of the algorithm
88 is stored in~$N$ cells, we need~$W\ge\log N$ just to be able to read the input.
89 On the other hand, we are interested in polynomial-time algorithms only, so $\Theta(\log N)$-bit
90 numbers should be sufficient. In practice, we pick~$W$ to be the larger of
91 $\Theta(\log N)$ and the size of integers used in the algorithm's input and output.
92 We will call an integer which fits in a~single memory cell a~\df{machine word.}
95 Both restrictions easily avoid the problems of unbounded parallelism. The first
96 choice is theoretically cleaner and Cook et al.~show nice correspondences to the
97 standard complexity classes, but the calculations of time and space complexity tend
98 to be somewhat tedious. What more, when compared with the RAM with restricted
99 word size, the complexities are usually exactly $\Theta(W)$ times higher.
100 This does not hold in general (consider a~program which uses many small numbers
101 and $\O(1)$ large ones), but it is true for the algorithms we are interested in.
102 Therefore we will always assume that the operations have unit cost and we make
103 sure that all numbers are limited by the available word size.
106 As for the choice of RAM operations, the following three instruction sets are often used:
109 \:\df{Word-RAM} --- allows the ``C-language operators'', i.e., addition,
110 subtraction, multiplication, division, remainder, bitwise $\band$, $\bor$, exclusive
111 $\bor$ ($\bxor$) and negation ($\bnot$), and bitwise shifts ($\shl$ and~$\shr$).
112 \:\df{${\rm AC}^0$-RAM} --- allows all operations from the class ${\rm AC}^0$, i.e.,
113 those computable by constant-depth polynomial-size boolean circuits with unlimited
114 fan-in and fan-out. This includes all operations of the Word-RAM except for multiplication,
115 division and remainders, and also many other operations like computing the Hamming
116 weight (number of bits set in a~given number).
117 \:Both restrictions at once.
120 Thorup discusses the usual techniques employed by RAM algorithms in~\cite{thorup:aczero}
121 and he shows that they work on both Word-RAM and ${\rm AC}^0$-RAM, but the combination
122 of the two restrictions is too weak. On the other hand, the intersection of~${\rm AC}^0$
123 with the instruction set of modern processors is already strong enough (e.g., when we
124 add some floating-point operations and multimedia instructions available on the Intel's
125 Pentium~4~\cite{intel:pentium}).
127 We will therefore use the Word-RAM instruction set, mentioning differences from the
128 ${\rm AC}^0$-RAM where necessary.
131 When speaking of the \df{RAM,} we implicitly mean the version with numbers limited
132 by a~specified word size of $W$~bits, unit cost of operations and memory cells and the instruction
133 set of the Word-RAM. This corresponds to the usage in recent algorithmic literature,
134 although the authors rarely mention the details. In some cases, a~non-uniform variant
135 of the Word-RAM is considered as well (e.g., in~\cite{hagerup:dd}):
137 \defn\id{nonuniform}%
138 A~Word-RAM is called \df{weakly non-uniform,} if it is equipped with $\O(1)$-time
139 access to a~constant number of word-sized constants, which depend only on the word
140 size. These are called \df{native constants} and they are available in fixed memory
141 cells when the program starts. (By analogy with the high-level programming languages,
142 these constants can be thought of as computed at ``compile time.'')
145 The \df{Pointer Machine (PM)} also does not have any well established definition. The
146 various kinds of pointer machines are mapped by Ben-Amram in~\cite{benamram:pm},
147 but unlike the RAM's they turn out to be equivalent up to constant slowdown.
148 Our definition will be closely related to the \em{linking automaton} proposed
149 by Knuth in~\cite{knuth:fundalg}, we will only adapt it to use RAM-like
150 instructions instead of an~opaque control unit.
152 The PM works with two different types of data: \df{symbols} from a~finite alphabet
153 and \df{pointers}. The memory of the machine consists of a~fixed amount of \df{registers}
154 (some of them capable of storing a~single symbol, each of the others holds a~single pointer)
155 and an~arbitrary amount of \df{cells}. The structure of all cells is the same: each of them
156 again contains a~fixed number of fields for symbols and pointers. Registers can be addressed
157 directly, the cells only via pointers --- by using a~pointer stored either in a~register,
158 or in a~cell pointed to by a~register (longer chains of pointers cannot be followed in
161 We can therefore view the whole memory as a~directed graph, whose vertices
162 correspond to the cells (the registers are stored in a~single special cell).
163 The outgoing edges of each vertex correspond to pointer fields of the cells and they are
164 labelled with distinct labels drawn from a~finite set. In addition to that,
165 each vertex contains a~fixed amount of symbols. The program can directly access
166 vertices within distance~2 from the register vertex.
168 The program is a~finite sequence of instructions of the following kinds:
171 \:\df{symbol instructions,} which read a~pair of symbols, apply an~arbitrary
172 function on them and write the result to a~symbol register or field;
173 \:\df{pointer instructions} for assignment of pointers to pointer registers/fields
174 and for creation of new memory cells (a~pointer to the new cell is assigned
176 \:\df{control instructions} --- similarly to the RAM; conditional jumps can decide
177 on~arbitrary unary relations on symbols and compare pointers for equality.
180 Time and space complexity are defined in the straightforward way: all instructions
181 have unit cost and so do all memory cells.
183 Both input and output of the machine are passed in the form of a~linked structure
184 pointed to by a~designated register. For example, we can pass graphs back and forth
185 without having to encode them as strings of numbers or symbols. This is important,
186 because with the finite alphabet of the~PM, symbolic representations of graphs
187 generally require super-linear space and therefore also time.\foot{%
188 The usual representation of edges as pairs of vertex labels uses $\Theta(m\log n)$ bits
189 and as a~simple counting argument shows, this is asymptotically optimal for general
190 sparse graphs. On the other hand, specific families of sparse graphs can be stored
191 more efficiently, e.g., by a~remarkable result of Tur\'an~\cite{turan:succinct},
192 planar graphs can be encoded in~$\O(n)$ bits. Encoding of dense graphs is of
193 course trivial as the adjacency matrix has only~$\Theta(n^2)$ bits.}
196 Compared to the RAM, the PM lacks two important capabilities: indexing of arrays
197 and arithmetic instructions. We can emulate both with poly-logarithmic slowdown,
198 but it will turn out that they are rarely needed in graph algorithms. We are
199 also going to prove that the RAM is strictly stronger, so we will prefer to
200 formulate our algorithms in the PM model and use RAM only when necessary.
203 Every program for the Word-RAM with word size~$W$ can be translated to a~PM program
204 computing the same with $\O(W^2)$ slowdown (given a~suitable encoding of inputs and
205 outputs, of course). If the RAM program does not use multiplication, division
206 and remainder operations, $\O(W)$~slowdown is sufficient.
209 Represent the memory of the RAM by a~balanced binary search tree or by a~radix
210 trie of depth~$\O(W)$. Values are encoded as~linked lists of symbols pointed
211 to by the nodes of the tree. Both direct and indirect accesses to the memory
212 can therefore be done in~$\O(W)$ time. Use standard algorithms for arithmetic
213 on big numbers: $\O(W)$ per operation except for multiplication, division and
214 remainders which take $\O(W^2)$.\foot{We could use more efficient arithmetic
215 algorithms, but the quadratic bound is good enough for our purposes.}
218 \FIXME{Add references, especially to the unbounded parallelism remark.}
221 Every program for the PM running in polynomial time can be translated to a~program
222 computing the same on the Word-RAM with only $\O(1)$ slowdown.
225 Encode each cell of the PM's memory to $\O(1)$ integers. Store the encoded cells to
226 the memory of the RAM sequentially and use memory addresses as pointers. As the symbols
227 are finite and there is only a~polynomial number of cells allocated during execution
228 of the program, $\O(\log N)$-bit integers suffice ($N$~is the size of the program's input).
232 There are also \df{randomized} versions of both machines. These are equipped
233 with an~additional instruction for generating a~single random bit. The standard
234 techniques of design and analysis of randomized algorithms apply (see for
235 example Motwani and Raghavan~\cite{motwani:randalg}).
237 \FIXME{Consult sources. Does it make more sense to generate random words at once on the RAM?}
240 There is one more interesting machine: the \df{Immutable Pointer Machine} (see
241 the description of LISP machines in \cite{benamram:pm}). It differs from the
242 ordinary PM by the inability to modify existing memory cells. Only the contents
243 of the registers are allowed to change. All cell modifications thus have to
244 be performed by creating a~copy of the particular cell with some fields changed.
245 This in turn requires the pointers to the cell to be updated, possibly triggering
246 a~cascade of further cell copies. For example, when a~node of a~binary search tree is
247 updated, all nodes on the path from that node to the root have to be copied.
249 One of the advantages of this model is that the states of the machine are
250 persistent --- it is possible to return to a~previously visited state by recalling
251 the $\O(1)$ values of the registers (everything else could not have changed
252 since that time) and ``fork'' the computations. This corresponds to the semantics
253 of pure functional languages, e.g., Haskell~\cite{jones:haskell}.
255 Unless we are willing to accept a~logarithmic penalty in execution time and space
256 (in fact, our emulation of the Word-RAM on the PM can be easily made immutable),
257 the design of efficient algorithms for the immutable PM requires very different
258 techniques. Therefore, we will concentrate on the imperative models instead
259 and refer the interested reader to the thorough treatment of purely functional
260 data structures in the Okasaki's monograph~\cite{okasaki:funcds}.
262 %--------------------------------------------------------------------------------
264 \section{Bucket sorting and contractions}\id{bucketsort}%
266 The Contractive Bor\o{u}vka's algorithm (\ref{contbor}) needed to contract a~given
267 set of edges in the current graph and flatten it afterwards, all this in time $\O(m)$.
268 We have spared the technical details for this section and they will be useful
269 in further algorithms, too.
271 As already suggested, the contractions can be performed by building an~auxiliary
272 graph and finding its connected components. Thus we will take care of the flattening
276 On the RAM, it is sufficient to sort the edges lexicographically (each edge viewed
277 as an~ordered pair of vertex identifiers with the smaller of the identifiers placed
278 first). We can do that by a two-pass bucket sort with~$n$ buckets corresponding
279 to the vertex identifiers.
281 However, there is a~catch in this. Suppose that we use the standard representation
282 of graphs by adjacency lists whose heads are stored in an array indexed by vertex
283 identifiers. When we contract and flatten the graph, the number of vertices decreases,
284 but if we inherit the original vertex identifiers, the arrays will still have the
285 same size. Hence we spend a~super-linear amount of time on scanning the increasingly
286 sparse arrays, most of the time skipping unused entries.
288 To avoid this, we have to renumber the vertices after each contraction to component
289 identifiers from the auxiliary graph and we create a~new vertex array. This way,
290 the representation of the graph will be kept linear with respect to the size of the
294 The pointer representation of graphs does not suffer from sparsity as the vertices
295 are always identified by pointers to per-vertex structures. Each such structure
296 then contains all attributes associated with the vertex, including the head of its
297 adjacency list. However, we have to find a~way how to perform bucket sorting
300 We will keep a~list of the per-vertex structures which defines the order of~vertices.
301 Each such structure will contain a~pointer to the head of the corresponding bucket,
302 again stored as a~list. Putting an~edge to a~bucket can be done in constant time then,
303 scanning all~$n$ buckets takes $\O(n+m)$ time.
305 \FIXME{Add an example of locally determined orders, e.g., tree isomorphism?}
307 %--------------------------------------------------------------------------------
309 \section{Data structures on the RAM}
311 There is a~lot of data structures designed specifically for the RAM, taking
312 advantage of both indexing and arithmetics. In many cases, they surpass the known
313 lower bounds for the same problem on the~PM and they often achieve constant time
314 per operation, at least when either the magnitude of the values or the size of
315 the data structure are suitably bounded.
317 A~classical result of this type are the trees of van Emde Boas~\cite{boas:vebt},
318 which represent a~subset of the integers $\{0,\ldots,U-1\}$, allowing insertion,
319 deletion and order operations (minimum, maximum, successor etc.) in time $\O(\log\log U)$,
320 regardless of the size of the subset. If we replace the heap used in the Jarn\'\i{}k's
321 algorithm (\ref{jarnik}) by this structure, we immediately get an~algorithm
322 for finding the MST in integer-weighted graphs in time $\O(m\log\log w_{max})$,
323 where $w_{max}$ is the maximum weight. We will show later that it is even
324 possible to achieve linear time complexity for arbitrary integer weights.
326 A~real breakthrough has been made by Fredman and Willard, who introduced
327 the Fusion trees~\cite{fw:fusion} which again perform membership and predecessor
328 operation on a~set of $n$~integers, but this time with complexity $\O(\log_W n)$
329 per operation on a~Word-RAM with $W$-bit words. This of course assumes that
330 each element of the set fits in a~single word. As $W$ must at least~$\log n$,
331 the operations take $\O(\log n/\log\log n)$ and we are able to sort $n$~integers
332 in time~$o(n\log n)$. This was a~beginning of a~long sequence of faster and
333 faster sorting algorithms, culminating with the work by Thorup and Han.
334 They have improved the time complexity of integer sorting to $\O(n\log\log n)$ deterministically~\cite{han:detsort}
335 and expected $\O(n\sqrt{\log\log n})$ for randomized algorithms~\cite{hanthor:randsort},
336 both in linear space.
338 Despite the recent progress, the corner-stone of most RAM data structures
339 is still the representation of data structures by integers introduced by Fredman
340 and Willard. It will also form a~basis for the rest of this chapter.
342 \FIXME{Add more history.}
344 %--------------------------------------------------------------------------------
346 \section{Bits and vectors}\id{bitsect}
348 In this rather technical section, we will show how the RAM can be used as a~vector
349 computer to operate in parallel on multiple elements, as long as these elements
350 fit in a~single machine word. At the first sight this might seem useless, because we
351 cannot require large word sizes, but surprisingly often the elements are small
352 enough relative to the size of the algorithm's input and thus also relative to
353 the minimum possible word size. Also, as the following lemma shows, we can
354 easily emulate slightly longer words:
356 \lemman{Multiple-precision calculations}
357 Given a~RAM with $W$-bit words, we can emulate all calculation and control
358 instructions of a~RAM with word size $kW$ in time depending only on the~$k$.
359 (This is usually called \df{multiple-precision arithmetics.})
362 We split each word of the ``big'' machine to $W'$-bit blocks, where $W'=W/2$, and store these
363 blocks in $2k$ consecutive memory cells. Addition, subtraction, comparison and
364 bitwise logical operations can be performed block-by-block. Shifts by a~multiple
365 of~$W'$ are trivial, otherwise we can combine each block of the result from
366 shifted versions of two original blocks.
367 To multiply two numbers, we can use the elementary school algorithm using the $W'$-bit
368 blocks as digits in base $2^{W'}$ --- the product of any two blocks fits
371 Division is harder, but Newton-Raphson iteration (see~\cite{ito:newrap})
372 converges to the quotient in a~constant number of iterations, each of them
373 involving $\O(1)$ multiple-precision additions and multiplications. A~good
374 starting approximation can be obtained by dividing the two most-significant
375 (non-zero) blocks of both numbers.
377 Another approach to division is using the improved elementary school algorithm as described
378 by Knuth in~\cite{knuth:seminalg}. It uses $\O(k^2)$ steps, but the steps involve
379 calculation of the most significant bit set in a~word. We will show below that it
380 can be done in constant time, but we have to be careful to avoid division instructions.
383 \notan{Bit strings}\id{bitnota}%
384 We will work with binary representations of natural numbers by strings over the
385 alphabet $\{\0,\1\}$: we will use $\(x)$ for the number~$x$ written in binary,
386 $\(x)_b$ for the same padded to exactly $b$ bits by adding leading zeroes,
387 and $x[k]$ for the value of the $k$-th bit of~$x$ (with a~numbering of bits such that $2^k[k]=1$).
388 The usual conventions for operations on strings will be utilized: When $s$
389 and~$t$ are strings, we write $st$ for their concatenation and
390 $s^k$ for the string~$s$ repeated $k$~times.
391 When the meaning is clear from the context,
392 we will use $x$ and $\(x)$ interchangeably to avoid outbreak of symbols.
395 The \df{bitwise encoding} of a~vector ${\bf x}=(x_0,\ldots,x_{d-1})$ of~$b$-bit numbers
396 is an~integer~$x$ such that $\(x)=\(x_{d-1})_b\0\(x_{d-2})_b\0\ldots\0\(x_0)_b$, i.e.,
397 $x = \sum_i 2^{(b+1)i}\cdot x_i$. (We have interspersed the elements with \df{separator bits.})
399 \notan{Vectors}\id{vecnota}%
400 We will use boldface letters for vectors and the same letters in normal type
401 for their encodings. The elements of a~vector~${\bf x}$ will be written as
402 $x_0,\ldots,x_{d-1}$.
405 If we want to fit the whole vector in a~single word, the parameters $b$ and~$d$ must satisty
406 the condition $(b+1)d\le W$.
407 By using multiple-precision arithmetics, we can encode all vectors satisfying $bd=\O(W)$.
408 We will now describe how to translate simple vector manipulations to sequences of $\O(1)$ RAM operations
409 on their codes. As we are interested in asymptotic complexity only, we prefer clarity
410 of the algorithms over saving instructions. Among other things, we freely use calculations
411 on words of size $\O(bd)$, assuming that the Multiple-precision lemma comes to save us
415 First of all, let us observe that we can use $\band$ and $\bor$ with suitable constants
416 to write zeroes or ones to an~arbitrary set of bit positions at once. These operations
417 are usually called \df{bit masking}. Also, any element of a~vector can be extracted or
418 replaced by a~different value in $\O(1)$ time by masking and shifts.
421 \def\setslot#1{\setbox0=#1\slotwd=\wd0}
422 \def\slot#1{\hbox to \slotwd{\hfil #1\hfil}}
423 \def\[#1]{\slot{$#1$}}
424 \def\9{\rack{\0}{\hss$\cdot$\hss}}
428 \halign{\hskip 0.15\hsize\hfil $ ##$&\hbox to 0.6\hsize{${}##$ \hss}\cr
433 \algn{Operations on vectors with $d$~elements of $b$~bits each}\id{vecops}
437 \:$\<Replicate>(\alpha)$ --- creates a~vector $(\alpha,\ldots,\alpha)$:
439 \alik{\<Replicate>(\alpha)=\alpha\cdot(\0^b\1)^d. \cr}
441 \:$\<Sum>(x)$ --- calculates the sum of the elements of~${\bf x}$, assuming that
442 the result fits in $b$~bits:
444 \alik{\<Sum>(x) = x \bmod \1^{b+1}. \cr}
446 This is correct because when we calculate modulo~$\1^{b+1}$, the number $2^{b+1}=\1\0^{b+1}$
447 is congruent to~1 and thus $x = \sum_i 2^{(b+1)i}\cdot x_i \equiv \sum_i 1^i\cdot x_i \equiv \sum_i x_i$.
448 As the result should fit in $b$~bits, the modulo makes no difference.
450 If we want to avoid division, we can use double-precision multiplication instead:
452 \setslot{\hbox{~$\0x_{d-1}$}}
455 \def\dd{\slot{$\cdots$}}
456 \def\vd{\slot{$\vdots$}}
457 \def\rule{\noalign{\medskip\nointerlineskip}$\hrulefill$\cr\noalign{\nointerlineskip\medskip}}
460 \[\0x_{d-1}] \dd \[\0x_2] \[\0x_1] \[\0x_0] \cr
461 *~~ \z \dd \z\z\z \cr
463 \[x_{d-1}] \dd \[x_2] \[x_1] \[x_0] \cr
464 \[x_{d-1}] \[x_{d-2}] \dd \[x_1] \[x_0] \. \cr
465 \[x_{d-1}] \[x_{d-2}] \[x_{d-3}] \dd \[x_0] \. \. \cr
466 \vd\vd\vd\vd\.\.\.\cr
467 \[x_{d-1}] \dd \[x_2]\[x_1]\[x_0] \. \. \. \. \cr
469 \[r_{d-1}] \dd \[r_2] \[r_1] \[s_d] \dd \[s_3] \[s_2] \[s_1] \cr
472 This way, we even get the vector of all partial sums:
473 $s_k=\sum_{i=0}^{k-1}x_i$, $r_k=\sum_{i=k}^{d-1}x_i$.
475 \:$\<Cmp>(x,y)$ --- element-wise comparison of~vectors ${\bf x}$ and~${\bf y}$,
476 i.e., a~vector ${\bf z}$ such that $z_i=1$ if $x_i<y_i$ and $z_i=0$ otherwise.
478 We replace the separator zeroes in~$x$ by ones and subtract~$y$. These ones
479 change back to zeroes exactly at the positions where $x_i<y_i$ and they stop
480 carries from propagating, so the fields do not interact with each other:
482 \setslot{\vbox{\hbox{~$x_{d-1}$}\hbox{~$y_{d-1}$}}}
483 \def\9{\rack{\0}{\hss ?\hss}}
485 \1 \[x_{d-1}] \1 \[x_{d-2}] \[\cdots] \1 \[x_1] \1 \[x_0] \cr
486 -~ \0 \[y_{d-1}] \0 \[y_{d-2}] \[\cdots] \0 \[y_1] \0 \[y_0] \cr
488 \9 \[\ldots] \9 \[\ldots] \[\cdots] \9 \[\ldots] \9 \[\ldots] \cr
491 It only remains to shift the separator bits to the right positions, negate them
492 and mask out all other bits.
494 \:$\<Rank>(x,\alpha)$ --- returns the number of elements of~${\bf x}$ which are less than~$\alpha$,
495 assuming that the result fits in~$b$ bits:
498 \<Rank>(x,\alpha) = \<Sum>(\<Cmp>(x,\<Replicate>(\alpha))). \cr
501 \:$\<Insert>(x,\alpha)$ --- inserts~$\alpha$ into a~sorted vector $\bf x$:
503 We calculate $k = \<Rank>(x,\alpha)$ first, then insert~$\alpha$ as the $k$-th
504 field of~$\bf x$ using masking operations and shifts.
506 \:$\<Unpack>(\alpha)$ --- creates a~vector whose elements are the bits of~$\(\alpha)_d$.
507 In other words, inserts blocks~$\0^b$ between the bits of~$\alpha$. Assuming that $b\ge d$,
508 we can do it as follows:
511 \:$x\=\<Replicate>(\alpha)$.
512 \:$y\=(2^{b-1},2^{b-2},\ldots,2^0)$. \cmt{bitwise encoding of this vector}
514 \:Return $\<Cmp>(z,y)$.
517 Let us observe that $z_i$ is either zero or equal to~$y_i$ depending on the value
518 of the $i$-th bit of the number~$\alpha$. Comparing it with~$y_i$ normalizes it
519 to either zero or one.
521 \:$\<Unpack>_\pi(\alpha)$ --- like \<Unpack>, but changes the order of the
522 bits according to a~fixed permutation~$\pi$: The $i$-th element of the
523 resulting vector is equal to~$\alpha[\pi(i)]$.
525 Implemented as above, but with a~mask $y=(2^{\pi(b-1)},\ldots,2^{\pi(0)})$.
527 \:$\<Pack>(x)$ --- the inverse of \<Unpack>: given a~vector of zeroes and ones,
528 it produces a~number whose bits are the elements of the vector (in other words,
529 it crosses out the $\0^b$ blocks).
531 We interpret the~$x$ as an~encoding of a~vector with elements one bit shorter
532 and we sum these elements. For example, when $n=4$ and~$b=4$:
534 \setslot{\hbox{$x_3$}}
536 \def\|{\hskip1pt\vrule height 10pt depth 4pt\hskip1pt}
537 \def\.{\hphantom{\|}}
540 \|\z\.\z\.\z\.\z\.\[x_3]\|\z\.\z\.\z\.\z\.\[x_2]\|\z\.\z\.\z\.\z\[x_1]\|\z\.\z\.\z\.\z\.\[x_0]\|\cr
541 \|\z\.\z\.\z\.\z\|\[x_3]\.\z\.\z\.\z\|\z\.\[x_2]\.\z\.\z\|\z\.\z\[x_1]\.\z\|\z\.\z\.\z\.\[x_0]\|\cr
544 However, this ``reformatting'' does not produce a~correct encoding of a~vector,
545 because the separator zeroes are missing. For this reason, the implementation
546 of~\<Sum> using modulo does not work correctly (it produces $\0^b$ instead of $\1^b$).
547 We therefore use the technique based on multiplication instead, which does not need
548 the separators. (Alternatively, we can observe that $\1^b$ is the only case
549 affected, so we can handle it separately.)
554 We can use the aforementioned tricks to perform interesting operations on individual
555 numbers in constant time, too. Let us assume for a~while that we are
556 operating on $b$-bit numbers and the word size is at least~$b^2$.
557 This enables us to make use of intermediate vectors with $b$~elements
560 \algn{Integer operations in quadratic workspace}\id{lsbmsb}
564 \:$\<Weight>(\alpha)$ --- computes the Hamming weight of~$\alpha$, i.e., the number of ones in~$\(\alpha)$.
566 We perform \<Unpack> and then \<Sum>.
568 \:$\<Permute>_\pi(\alpha)$ --- shuffles the bits of~$\alpha$ according
569 to a~fixed permutation~$\pi$.
571 We perform $\<Unpack>_\pi$ and \<Pack> back.
573 \:$\<LSB>(\alpha)$ --- finds the least significant bit of~$\alpha$,
574 i.e., the smallest~$i$ such that $\alpha[i]=1$.
576 By a~combination of subtraction with $\bxor$, we create a~number
577 which contains ones exactly at the position of $\<LSB>(\alpha)$ and below:
580 \alpha&= \9\9\9\9\9\1\0\0\0\0\cr
581 \alpha-1&= \9\9\9\9\9\0\1\1\1\1\cr
582 \alpha\bxor(\alpha-1)&= \0\9\9\9\0\1\1\1\1\1\cr
585 Then we calculate the \<Weight> of the result and subtract~1.
587 \:$\<MSB>(\alpha)$ --- finds the most significant bit of~$\alpha$ (the position
588 of the highest bit set).
590 Reverse the bits of the number~$\alpha$ first by calling \<Permute>, then apply \<LSB>
591 and subtract the result from~$b-1$.
596 As noted by Brodnik~\cite{brodnik:lsb} and others, the space requirements of
597 the \<LSB> operation can be reduced to linear. We split the input to $\sqrt{b}$
598 blocks of $\sqrt{b}$ bits each. Then we determine which blocks are non-zero and
599 identify the lowest such block (this is a~\<LSB> of a~number whose bits
600 correspond to the blocks). Finally we calculate the \<LSB> of this block. In
601 both calls to \<LSB,> we have a $\sqrt{b}$-bit number in a~$b$-bit word, so we
602 can use the previous algorithm. The same trick of course works for finding the
605 The following algorithm shows the details.
607 \algn{LSB in linear workspace}
610 \algin A~$w$-bit number~$\alpha$.
611 \:$b\=\lceil\sqrt{w}\,\rceil$. \cmt{size of a~block}
612 \:$\ell\=b$. \cmt{the number of blocks is the same}
613 \:$x\=(\alpha \band (\0\1^b)^\ell) \bor (\alpha \band (\1\0^b)^\ell)$.
615 \cmt{encoding of a~vector~${\bf x}$ such that $x_i\ne 0$ iff the $i$-th block is non-zero}%
616 \foot{Why is this so complicated? It is tempting to take $\alpha$ itself as a~code of this vector,
617 but we unfortunately need the separator bits between elements, so we create them and
618 relocate the bits we have overwritten.}
619 \:$y\=\<Cmp>(0,x)$. \cmt{$y_i=1$ if the $i$-th block is non-zero, otherwise $y_0=0$}
620 \:$\beta\=\<Pack>(y)$. \cmt{each block compressed to a~single bit}
621 \:$p\=\<LSB>(\beta)$. \cmt{the index of the lowest non-zero block}
622 \:$\gamma\=(\alpha \shr bp) \band \1^b$. \cmt{the contents of that block}
623 \:$q\=\<LSB>(\gamma)$. \cmt{the lowest bit set there}
624 \algout $\<LSB>(\alpha) = bp+q$.
628 We have used a~plenty of constants which depend on the format of the vectors.
629 Either we can write non-uniform programs (see \ref{nonuniform}) and use native constants,
630 or we can observe that all such constants can be easily manufactured. For example,
631 $(\0^b\1)^d = \1^{(b+1)d} / \1^{b+1} = (2^{(b+1)d}-1)/(2^{b+1}-1)$. The only exceptions
632 are the~$w$ and~$b$ in the LSB algorithm \ref{lsb}, which we are unable to produce