```pascal
let pop q = match !q with
| None ->
    invalid_arg "pop"
```

When `q` is not empty, it is convenient to treat the case in which `q` contains only one element separately. Here, the final element points to itself. It thus suffices to empty the queue, giving it the value `None` and returning the only element that it contained.

```pascal
| Some last when last.next == last ->
    q := None;
    last.elt
```

In the general case, when `q` contains at least two elements, we begin by determining the first element of the queue, `first`. We then remove it from the linked list by making the final element point to the second one. The value to be returned is contained in the cell `first`.

```pascal
| Some last ->
    let first = last.next in
    last.next <- first.next;
    first.elt
```

The complete code is given in program 58 (see following page).

## Complexity

It is clear that the operations `create`, `is_empty`, `push`, and `pop` all execute in constant time. A queue containing `N` elements occupies three words per element in memory (being a record with two fields, `elt` and `next`), plus two words for the reference.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Program 58 — Imperative queues using linked lists.

```pascal
type 'a cell = { elt : 'a; mutable next : 'a cell }

type 'a t = 'a cell option ref

let create () =
    ref None

let is_empty q =
    !q = None

let push x q = match !q with
| None ->
    let rec c = { elt = x; next = c } in
        q := Some c
| Some last ->
    let c = { elt = x; next = last.next } in
        last.next <- c;
        q := Some c

let pop q = match !q with
| None ->
    invalid_arg "pop"
| Some last when last.next == last ->
    q := None;
    last.elt
| Some last ->
    let first = last.next in
        last.next <- first.next;
        first.elt
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
<figure><img src="image_96.png" /></figure>

**For Further Information**

The module `Queue` of the OCaml standard library also uses a singly linked list, but in a slightly different way. First, it uses a type with two constructors, `Nil` and `Cons`, instead of using the combination of a record type and an option type. Second, it maintains two pointers, to the first and last elements, rather than using a circular list.

## 6.2 Persistent Queues

In this section, we present a data structure of persistent queues. The signature of such a structure is given in program 59. We see, in particular, that the operations `push` and `pop` return a modified version of the data structure.

A naive idea would consist in representing a queue by a list. This solution is unfortunately very inefficient, because either the insertion or the extraction operation would have to be performed at the end of the list. To implement these two operations more efficiently, a simple idea consists in using not one but two lists. The first list contains the elements at the front of the queue, in order. The second list contains the elements at the back of the queue, in reverse order. In this way, both insertion and extraction are done at the head of a list. For example, the queue containing the elements 1, 2, 3, 4, 5, inserted in this order, may be represented by the two lists `[1; 2]` and `[5; 4; 3]`. However, the queue could equally be represented by the two lists `[1]` and `[5;4;3;2]`, or even by `[1;2;3;4;5]` and `[]`, which are all equivalent representations of the same queue.

We therefore introduce the type of persistent queues simply as a synonym for a pair of lists.

```hl
type 'a t = 'a list * 'a list
```

The first list represents the front of the queue. The second represents the back of the queue. The definitions of the empty queue `empty` and the function `is_empty` are clear.

```hl
let empty = [], []
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 59 — Minimal signature for persistent queues.

```python
module type PersistentQueue = sig

type 'a t
val empty : 'a t
val is_empty : 'a t -> bool
val push : 'a -> 'a t -> 'a t
val pop : 'a t -> 'a * 'a t
end
```

```python
let is_empty = function
| [], [] -> true
| _ -> false
```

To insert an element x, it suffices to add it to the head of the second list.

```python
let push x (o, i) =
  (o, x :: i)
```

Extraction is a more delicate operation. There are three possible cases. First, when the queue is empty, we raise an exception.

```python
let pop = function
| [], [] ->
    invalid_arg "pop"
```

Second, when the first list, which represents the front of the queue, contains at least one element, it suffices to extract this element.

```python
| x :: o, i ->
  x, (o, i)
```

Third, when all the elements are in the back of the queue, it suffices to reverse the list i and proceed as in the previous case, where the element to be extracted was at the head of the list. The tail of the list thus obtained becomes the

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
new front of the queue, and the back of the new queue no longer contains any
elements.

| [], i ->
 match List.rev i with
 | x :: o -> x, (o, [])

Since the list `i` is not empty, the list `List.rev i` is also not empty. Rather
than write a non-exhaustive pattern matching, which would trigger a compiler
warning, we indicate explicitly that the case of an empty list cannot occur.

| [] -> assert false

The complete code is given in program 60.

Complexity

It is clear that the operations `is_empty` and `push` execute in constant time. The
operation `pop` can, however, have cost $O(N)$ for a queue containing $N$ elements.
This is the case when all the elements are at the back of the queue so that the
list has to be inverted. Nevertheless, if we consider a sequence of $M$ successive
`push` and `pop` operations, starting from an empty queue, then the total cost
cannot be greater than $O(M)$. Indeed, each element may be involved in, at
most, one list reversal. The `pop` operation therefore has amortized cost $O(1)$.
The space complexity is the same as in case of imperative queues, with three
words per element of the queue.

6.3 Imperative Priority Queues

In this section and the next, we will consider queues in which the elements are
associated with priorities. In such queues, called *priority queues*, the elements
are removed in the order fixed by their priority rather than in the order in which
they were entered. We present here an *imperative* data structure of priority
queues, whose signature is given in program 61 (see the following page).

In this interface, the notion of minimality coincides with the notion of great-
est priority. Unlike with ordinary queues, here we distinguish the access and

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 60 — Persistent queues represented by pairs of lists.

```javascript
type 'a t = 'a list * 'a list
```

```javascript
let empty =
    ([], [])
```

```javascript
let is_empty = function
  | [], [] -> true
  | _ -> false
```

```javascript
let push x (o, i) =
    (o, x :: i)
```

```javascript
let pop = function
  | [], [] ->
    invalid_arg "pop"
  | x :: o, i ->
    x, (o, i)
  | [], i ->
    match List.rev i with
    | x :: o -> x, (o, [])
    | [] -> assert false
```
removal of the first element using two separate operations. This is for the sake of efficiency, as will be explained later. Accordingly, the function `get_min` returns the element of the queue with the highest priority, and the function `remove_min` removes it.

To implement an efficient priority queue, it is necessary to use a more complex data structure than in case of ordinary queues. One solution consists in organizing the elements as a *heap*. A heap is a binary tree in which the element stored at each node has a higher priority than the two elements situated immediately below it. The element with the highest priority is therefore found at the root. Thus, a heap containing the elements {3, 7, 9, 12, 21}, ordered by magnitude, may have the following shape:

<figure><img src="image_97.png" /></figure>

We note that there are other heaps containing the same elements. However, for the sake of efficiency, it is preferable to choose a heap of minimal height. Heaps represented by complete binary trees (in which all levels are filled except, perhaps, the last one) have minimal height. It turns out that a complete binary tree can be easily represented using an array. The idea consists in numbering the nodes of the tree from top to bottom and left to right, beginning from 0. The result of this numbering of the heap above yields the following labeling:

<figure><img src="image_98.png" /></figure>

This numbering permits the representation of the heap in an array of five elements as follows:

<table><thead><tr><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>3</td><td>7</td><td>12</td><td>21</td><td>9</td></tr></tbody></table>

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 61 — Minimal signature for imperative priority queues.

```pascal
module type ImperativePriorityQueue = sig
type t
  type elt
  val create : unit -> t
  val is_empty : t -> bool
  val add : elt -> t -> unit
  val get_min : t -> elt
  val remove_min : t -> unit
end
```

In general, the root of the tree occupies the cell at index 0, and the roots of the two subtrees of the node stored at cell *i* are stored in cells 2*i* + 1 and 2*i* + 2, respectively. Conversely, the parent of node *i* is stored at $\lfloor (i - 1)/2 \rfloor$.

One problem remains: We do not know the size of the priority queue before-hand and therefore cannot fix in advance the maximum size of the array. An elegant solution consists in using resizeable arrays. Such arrays were presented in section 4.1 *Resizeable Arrays*. It then suffices to write a functor parametrized by a structure of resizeable arrays. A signature for such arrays, *ResizeableArray*, is given in program 17, at the start of chapter 4. As far as the elements are concerned, we require a type equipped with a total order and a default value. The corresponding signature, *OrderedWithDummy*, is given in program 62.

The functor of priority queues then takes the following form:

```pascal
module Make(X: OrderedWithDummy)(A: ResizeableArray)
  : ImperativePriorityQueue with type elt = X.t =
    struct
```

Module X is then that of elements, and module A that of resizeable arrays.

The types elt and t are respectively synonyms for X.t and elt A.t.

```pascal
type elt = X.t
type t = elt A.t
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
**Program 62 — Ordered elements with a default value.**

```python
module type OrderedWithDummy = sig
type t
val compare: t -> t -> int
val dummy : t
end
```

In fact, a heap is nothing but a resizeable array. The function `create` constructs an empty array with `A.make`, and the function `is_empty` tests whether the array is empty.

```python
let create () = A.make 0 X.dummy
let is_empty h = A.length h = 0
```

The function `get_min` returns the root of the heap, if it exists, and raises an exception otherwise.

```python
let get_min h =
    if A.length h = 0 then invalid_arg "get_min";
    A.get h 0
```

The code is given in program 63.

## Inserting an Element

Inserting an element `x` in a heap `h` consists in adding `x` into the tree that represents `h` while preserving the heap structure. One approach consists in adding `x` at the bottom right and then moving it up the tree as long as is necessary. In our implementation, this boils down to extending the array by one cell, inserting the value `x` there, and then moving `x` up until it reaches the correct position. For this, we use the following algorithm: if `x` is smaller than its parent, we exchange the two values and repeat.

For example, let us consider again the following heap:

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
<figure><img src="image_99.png" /></figure>

The insertion of the element 1 in this heap is performed in three steps:

<figure><img src="image_100.png" /></figure>

We begin by writing a function `move_up` that inserts an element `x` in the heap `h`, starting at position `i`. This function assumes that the tree of root `i`, obtained by placing `x` at `i`, is a heap. The function `move_up` first considers the case in which `i` is 0, that is, when we have reached the root. It then suffices to insert `x` at position `i`.

```hkl
let rec move_up h x i =
    if i = 0 then
        A.set h i x
```

If, however, we are at an internal node, we calculate the index `fi` of the parent of `i` and the value `y` stored in this node. If `y` is greater than `x`, we need to move `x` up while moving `y` down to position `i`. We then call `move_up` recursively beginning from `fi`.

```hkl
else
    let fi = (i - 1) / 2 in
    let y = A.get h fi in
    if X comparisons to y > 0 then begin
        A.set h i y;
        move_up h x fi
    end
```

If, however, `y` is less than or equal to `x`, then `x` is at its final position, and it suffices to store it there.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
else A.set h i x

The function `add` proceeds in two steps. It increases the size of the array by one unit, adding a cell at the end of the array. Then it calls the function `move_up`, beginning from this cell.

let add x h =
let n = A.length h in
A_resize h (n + 1);
move_up h x n

The code is given in program 63 (see next page).

## Removing the Smallest Element

Removing the smallest element of a heap is more delicate. The reason for this is that we remove the root of the tree and must, therefore, find an element to replace it. The idea consists in choosing the element at the bottom right of the heap, that is, the element occupying the last cell of the array, placing it at the root, and then moving it down the heap until it reaches its place. This is roughly what we did when inserting an element, where we moved the new element up.

Let us assume, for example, that we wish to remove the smallest element of the following heap:

<figure><img src="image_101.png" /></figure>

We replace the root, 1, by the element at the bottom right, 8. Then we make 8 move down until it reaches its place. For this, we compare 8 with the roots `a` and `b` of the two subtrees. If `a` and `b` are both greater than 8, the descent is complete. Otherwise, we exchange 8 with the smaller of the two nodes `a` and `b`, and continue the descent. In the example, 8 is successively exchanged with 4 and 5:

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Program 63 — An imperative heap data structure (1/2).

```pascal
module type OrderedWithDummy = sig
type t
val compare: t -> t -> int
tval dummy: t
end

module Make(X: OrderedWithDummy)(A: ResizeableArray)
  : ImperativePriorityQueue with type elt = X.t =
    struct
      type elt = X.t
      type t = elt A.t

    let create () =
      A.make 0 X.dummy

    let is_empty h =
      A.length h = 0

    let get_min h =
      if A.length h = 0 then invalid_arg "get_min";
      A.get h 0

    let rec move_up h x i =
      if i = 0 then A.set h i x else
        let fi = (i - 1) / 2 in
        let y = A.get h fi in
        if X compare y x > 0 then begin
          A.set h i y;
          move_up h x fi
        end else
          A.set h i x

    let add x h =
      let n = A.length in Aresize h (n + 1) move_up h x n
```
<figure><img src="image_102.png" /></figure>

We begin by writing a function that compares two nodes, given by their indices `l` and `r`, and returns the smaller of the two.

```python
let min h l r =
    if X compare (A.get h r) (A.get h l) < 0 then r else l
```

Next, we write a function `smallest_node` that determines if a value `x`, situated at node `i`, needs to move down the heap or not. The value `x` is not that of the node `i` because we do not know yet if this will be its final position. However, the function `smallest_node` acts as if the node `i` had the value `x`. The function `smallest_node` accordingly takes as parameters the heap `h`, the value `x`, and the index `i`. We begin by calculating the left child `l` of node `i`.

```python
let smallest_node h x i =
    let l = 2 * i + 1 in
```

If the node `i` has no children, we return `i`.

```python
let n = A.length h in
if l >= n then
    i
```

If, however, the node `l` is part of the heap, we determine the index `j` of the smallest child of `i`. For this, it is convenient to determine if the right child `r = l+1` exists. If it does, we compare the nodes `l` and `r` with the function `min`. Otherwise, the index `j` is equal to `l`.

```python
else
    let r = l + 1 in
    let j = if r < n then min h l r else l in
```

Finally, we compare the value of the node `j` with that of `x`.

```python
if X compare (A.get h j) x < 0 then j else i
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Moving down the heap is implemented by a recursive function `move_down` that takes as parameters the heap `h`, the value `x` that is moving down, and the index `i` of the candidate node. We begin by calling the function `smallest_node` to calculate the index `j` that should receive the value `x`.

```pascal
let rec move_down h x i =
    let j = smallest_node h x i in
```

If `j = i`, the descent is complete and it suffices to store the value `x` at index `i`.

```pascal
if j = i then
    A.set h i x
```

Otherwise, we move the value of the node `j` up to the node `i` and then continue moving down from `j`.

```pascal
else begin
    A.set h i (A.get h j);
    move_down h x j
end
```

The function for removing the smallest element of a heap `h` is then the following:

```pascal
let remove_min h =
    let n = A.length h - 1 in
    if n < 0 then invalid_arg "remove_min";
    let x = A.get h n in
    A_resize h n;
    if n > 0 then move_down h x 0
```

We consider three cases: If the heap is empty, we raise an exception. Next, if it contains only one element, it is not necessary to call the function `move_down`. Finally, if the heap contains more than one element, it is worthwhile to resize the array with `A_resize`. The complete code is given in programs 63 (see page 281) and 64.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 64 — An imperative heap data structure (2/2).

```pascal
let min h 1 r =
    if X comparisons (A, get h r) (A, get h 1) < 0 then r else 1

let smallest_node h x i =
    let l = 2 * i + 1 in
    let n = A.length h in
    if l >= n then i else
        let r = l + 1 in
        let j = if r < n then min h l r else l in
        if X comparisons (A, get h j) x < 0 then j else i

let rec move_down h x i =
    let j = smallest_node h x i in
    if j = i then A.set h i x
    else begin A.set h i (A.get h j); move_down h x j end

let remove_min h =
    let n = A.length h - 1 in
    if n < 0 then invalid_argv "remove_min";
    let x = A.get h n in
    A_resize h n;
    if n > 0 then move_down h x 0
end
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Complexity

It is clear that the functions `create`, `is_empty`, and `get_min` have constant cost. In case of the functions `add` and `remove_min`, their cost is bounded by the height of the heap. Indeed, in one case, we move from the leaves towards the root, and in the other, we move from the root towards the leaves. Since the heap is a complete binary tree, its height is logarithmic in the number of elements. For a priority queue containing $N$ elements, the functions `add` and `remove_min` therefore have cost $O(\log N)$. A simpler way of looking at it consists in noting that the function `move_up` divides its argument $i$ by two at each recursive call. Therefore, it cannot be called more than $O(\log N)$ times. Similarly, the function `move_down` multiplies its argument by two at each recursive call and, therefore, cannot be called more than $O(\log N)$ times. Of course, all of this applies only in an amortized manner by virtue of the calls to `resize`, which are amortized $O(1)$, as we showed in chapter 4.

The space complexity depends on the implementation of the resizeable array. If we use the implementation presented at the start of chapter 4, where the size of the array was doubled when it had to be increased, the array contains at most two times more elements than required. A priority queue containing $N$ elements therefore occupies at most $2N$ memory words.

6.4 Persistent Priority Queues

In this section, we present a data structure of *persistent* priority queues, whose signature is given in program 65.

As in the previous section, we use a heap data structure, except that here it is directly represented by a binary tree. Unlike in case of AVL trees, no extra information is stored in the nodes for the purposes of balancing. We will see later that these heaps nevertheless have amortized logarithmic complexity. We say these trees are *self-balancing*.

The functor of persistent priority queues therefore takes the following form:

module Make(X: Ordered) :
PersistentPriorityQueue with type elt = X.t =

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 65 — Minimal signature for persistent priority queues.

```python
module type PersistentPriorityQueue = sig
type t
    type elt
    val empty : t
    val is_empty : t -> bool
    val add : elt -> t -> t
    val get_min : t -> elt
    val remove_min : t -> t
end
```

```python
struct
    type elt = X.t
```

The argument `X`, of type `Ordered`, defines the type of elements. These must be equipped with a comparison function corresponding to the priority order (see program 33 of chapter 5 for the signature `Ordered`). The type `elt` of the elements of the queue is thus a synonym for `X.t`. The type `t` is the type of binary trees.

```python
type t = Empty | Node of t * elt * t
```

The definitions of the empty heap and the function `is_empty` are clear.

```python
let empty = Empty
let is_empty h =
    h = Empty
```

The structure of the heap gives immediate access to the smallest element.

```python
let get_min = function
    | Empty -> invalid_arg "get_min"
    | Node (., x, _ ) -> x
```

All the subtlety of this structure lies in the function `merge` that merges two heaps. Assuming we have already written this function, it is easy to write the

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
functions `add` and `remove_min`. Indeed, the addition of an element `x` to a heap `h` consists in merging the heap containing the single element `x` with `h`.

```h
let add x h =
    merge (Node (Empty, x, Empty)) h
```

Similarly, the removal of the smallest element of a heap consists in merging the two child heaps of the root.

```h
let remove_min = function
  | Empty -> invalid_arg "remove_min"
  | Node (a,_, b) -> merge a b
```

All we have to do now is write the function `merge` that merges two heaps `ha` and `hb`. If one of the two heaps is empty, how to do so is clear.

```h
let rec merge ha hb = match ha, hb with
  | Empty, h | h, Empty ->
    h
```

If, however, none of the heaps is empty, we construct the resulting heap as follows. Its root is the smallest of the two roots of `ha` and `hb`. Assume that the root of `ha` is smaller.

```h
| Node (la, xa, ra), Node (lb, xb, rb) ->
  if X comparisons xa xb <= 0 then
```

We must now determine the two subtrees of `Node` `(`, `xa`, `_`)`. There are several possibilities that result from recursively calling `merge` on two of the three trees `la`, `ra`, and `hb`, and from choosing the result as the left or right subtree. Among these possibilities, we choose the one that performs a rotation of the subtrees from right to left, so as to ensure self-balancing. Thus, `ra` takes the position of `la`, and `la` is merged with `hb`.

```h
Node (ra, xa, merge la hb)
```

The other situation, in which the root of `hb` is smaller, is symmetrical.

```h
else Node (rb, xb, merge lb ha)
```

The complete code of persistent priority queues is given in program 66.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Program 66 — Persistent priority queues.

```python
module Make(X: Ordered):
    PersistentPriorityQueue with type elt = X.t =
struct
    type elt = X.t
    type t = Empty | Node of t * elt * t

    let empty =
        Empty

    let is_empty h =
        h = Empty

    let get_min = function
        | Empty -> invalid_arg "get_min"
        | Node (., x,_) -> x

    let rec merge ha hb = match ha, hb with
        | Empty, h | h, Empty ->
            h
        | Node (la, xa, ra), Node (lb, xb, rb) ->
            if X comparisons xa xb <= 0 then
                Node (ra, xa, merge la hb)
            else
                Node (rb, xb, merge lb ha)

    let add x h =
        merge (Node (Empty, x, Empty)) h

    let remove_min = function
        | Empty -> invalid_arg "remove_min"
        | Node (a,_, b) -> merge a b
end
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Complexity

The functions `empty`, `is_empty`, and `get_min` clearly have constant cost. The analysis is subtler in case of the function `merge`, and hence `add` and `remove_min`. In the worst case, `merge` can have a cost as large as the total number of elements. However, for a sequence of `add` and `remove_min` operations, executed successively on a heap, this worst case cannot occur each time. We can show that the amortized cost of each operation is, in fact, $O(\log N)$, where N is the number of elements in the heap. For an analysis of this complexity, we refer the reader to the article by Sleator and Tarjan introducing self-balancing heaps [22] or the chapter on heaps by Okasaki in *The Fun of Programming* [10].

6.5 Exercises

Linked Lists

6.1 Write a function iter : ('a -> unit) -> 'a t -> unit that traverses the elements of a queue in the order in which they were inserted.

6.2 Write a function transfer : 'a t -> 'a t -> unit that takes two queues q1 and q2 as arguments, moves all the elements of q1 into q2, and empties the queue q1. Make sure that this operation executes in constant time.

6.3 Modify the queue structure of program 58 to keep track of the total num-
ber of elements throughout. You can make type t into a record with a mutable
field containing the number of elements of the queue.

Pairs of Lists

6.4 Consider the case of a queue q, in which the list representing the front is empty. If the operation pop is performed multiple times, then the reversal of the list i, representing the back of the queue, will also be performed multiple times. To remedy this, we can memoize the reversal of i in the original queue q by means of a side effect. Accordingly, we adopt an equivalent representation of the queue q that is still persistent. Modify the type 'a t as follows:

— openly licensed via CC BY SA 4.0 —
```type
    type 'a t = { mutable o : 'a list; mutable i : 'a list }
```

Reimplement the operations on these queues. The mutable nature of the fields of this type does not allow defining a polymorphic value empty. It suffices to replace it by a function `create : unit -> 'a t`.

**6.5** Note that the representation of queues as pairs of lists is symmetrical. We can therefore perform insertion and removal operations from the two ends of the queue with equal facility. We call the resulting structure a double-ended queue or `dequeue`. We can name the four operations as follows:

```hkl
val push_front : 'a -> 'a t -> 'a t
val push_back : 'a -> 'a t -> 'a t
val pop_front : 'a t -> 'a * 'a t
val pop_back : 'a t -> 'a * 'a t
```

The functions `push_back` and `pop_front` correspond to the functions `push` and `pop` already written. Write the functions `push_front` and `pop_back`.

**6.6** The double-ended queues of the previous exercise may prove inefficient. Indeed, alternating removals from the two ends of the queue can end up requiring many list reversals. To remedy this, we can try to balance the lengths of the two lists, so as to ensure a minimum proportion of elements at each end. For this, we let $c \ge 2$ be a constant and impose the following invariant on the lengths $l_i$ and $l_o$ of the two lists:

$$l_i \le c \times l_o + 1 \quad \text{and} \quad l_o \le c \times l_i + 1$$

Modify the type `'a t` as follows, to keep track of the length of each list.

```type
    type 'a t = { li : int; i : 'a list; lo : int; o : 'a list }
```

Reimplement the operations for double-ended queues. For more details concerning this structure, and notably for its analysis, you can consult the book by Okasaki [19, sec. 8.4].

## Heap Structures

**6.7** You can use the heap structure to easily implement an efficient sorting algorithm, called `heapsort`. The idea is as follows: We insert all the elements

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
to be sorted in a heap, and then extract them successively with the functions `get_min` and `remove_min`. Add a function `sort`: `X.t array -> unit` to the functor `Make` to sort an array using this algorithm. If the complexity of the heap operations is logarithmic (which is the case when you use resizeable arrays of constant amortized complexity), then the complexity of this sorting algorithm is $O(n \log n)$, that is, optimal. Heapsort will be described in detail in chapter 12.

**6.8** Several different elements of a heap can have the same priority. This is, for example, the case when the elements are pairs for which the function `compare` takes into account only one of the two components. Show that the order of insertion in a priority queue is preserved, that is, that the elements of equal priority are extracted from the queue in the same order in which they were inserted.

**6.9** We wish to add an operation `remove` to remove an element from a priority queue. Several problems arise. On the one hand, the same element may occur several times in the queue. On the other, to implement the removal of an element $x$ efficiently, it is necessary to determine at which index it appears without traversing the entire array, which is not possible with the current structure.

We propose to resolve these two problems as follows. When an element is added in the priority queue, we return a “pointer” to the user, that may be used subsequently to remove this element. The signature is accordingly modified as follows:

```hkl
type pointer
val add: elt -> t -> pointer
```

Internally, the type `pointer` contains the index where the corresponding element of the array is stored. This index is updated when the element is moved within the array. In practice, it is simplest to also include the value of the element along with the index:

```hkl
type pointer = { elt: elt; mutable index: int }
```

The resizeable array then contains values of type `pointer` rather than `elt`. Modify the priority queue structure following this idea and provide a function `remove: pointer -> t -> unit`. Hint: To remove the element situated at index $i$, replace it with the final element of the resizeable array. It suffices then

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
to re-establish the heap property, by using either `move_up` or `move_down`. Show that the complexity of `remove` is in $O(\log n)$.

**6.10** In the context of the above exercise, a value of type `pointer` may reference an element that is no longer in the queue. A call to `remove` could therefore corrupt the structure. Remedy this problem by invalidating all pointers corresponding to an element that is no longer in the queue, so that the function `remove` is able to fail cleanly.

**6.11** It can be interesting to have priority queues that provide an operation `change_priority` to modify the priority of an existing element. Reusing the idea of exercise 6.9, provide the following new operation:

```hkl
val change_priority: pointer -> elt -> t -> unit
```

This operation modifies the value of an element already present in the queue. Show that the complexity of `change_priority` is in $O(\log n)$.

**6.12** In certain cases, we may wish to distinguish the priority of an element from its value. Modify the implementation of priority queues to give them the following signature:

```python
module Make(P: Ordered)(V: WithDummy):
    sig
        type t
        type priority = P.t
        type value = V.t
        val create: unit -> t
        val is_empty: t -> bool
        val add: priority -> value -> t -> unit
        val get_min: t -> priority * value
        val remove_min: t -> unit
    end
```

The type `P.t` is that of priorities and `V.t`, that of values. As the signatures indicate, the first is equipped with a function `compare` and the second with a default value.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Self-Balancing Heaps

6.13 Redo exercise 6.7 (heapsort) with self-balancing heaps. What is the complexity of this sorting algorithm?

6.14 Using a reference, wrap the code of program 66 in a module with the imperative signature of program 57.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by/4.0/) —
— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by/4.0/) —
Graphs

Graphs are a fundamental data structure in computer science. A graph consists of a set of *vertices* linked by *edges*. We are accustomed to visualizing graphs as illustrated in figure 7.1.

<figure><img src="image_103.png" /><figcaption>Figure 7.1: Example of a graph.</figcaption></figure>

Formally, a graph is a set $V$ of vertices and a set $E$ of edges, which are *unordered pairs* of vertices. If $\{x, y\} \in E$, we say that the vertices $x$ and $y$ are adjacent. As the adjacency relation is symmetric, we speak of an *undirected graph*.

We can also define the notion of *directed graph* by taking $E$ to be a set of *ordered pairs* of vertices rather than unordered pairs. We then speak of *arcs* rather than edges. If $(x, y) \in E$, we say that $y$ is a successor of $x$ and write

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by/4.0/) —
$x \to y$. An example of a directed graph is given in figure 7.2.

<figure><img src="image_104.png" /><figcaption>Figure 7.2: Example of a directed graph.</figcaption></figure>

An arc from a vertex to itself, as in this example, is called a *cycle*. The *indegree* (resp. *outdegree*) of a vertex is the number of arcs that point towards it (resp. that point away from it). A *path* from a vertex $u$ to a vertex $v$ is a sequence $x_0, \dots, x_n$ of vertices such as $x_0 = u$, $x_n = v$ and $x_i \to x_{i+1}$ for $0 \le i < n$. Such a path is of length $n$; it contains $n$ arcs.

Vertices, like arcs, may carry information. We then speak of a *labeled graph*. An example of a directed, labeled graph is given in figure 7.3.

<figure><img src="image_105.png" /><figcaption>Figure 7.3: Example of a directed, labeled graph.</figcaption></figure>

It is important to note that the label of a vertex is not the same as the vertex itself. In particular, two vertices may bear the same label. Formally, a labeled graph therefore comprises two additional functions giving, respectively, the label of a vertex in $V$ and the label of an arc in $E$.

In what follows, we will use the term edge for both directed and undirected graphs.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by/4.0/) —
Basic Operations on Graphs

Before defining data structures to represent graphs, it is necessary to consider
the operations that these structures should provide. We do this independently
of any particular representation, with a view to constructing generic algorithms.
We may distinguish the following basic operations:

* construction operations, such as the creation of a graph, the addition or removal of a vertex or an edge;

* access operations, such as testing whether a vertex or an edge belongs to a graph, accessing the degree of a vertex, its label, the number of vertices or edges;

* basic traversal operations, such as the traversal of all vertices, of all successors of a given vertex, or of all edges.

As we will see shortly, there is no single structure of graphs that provides all
these operations in an optimal manner. At the end of this chapter, we com-
pare different graph data structures according to the complexity of the different
operations above.

**7.1 Adjacency Matrix**

We consider in this section the case in which vertices are represented by integers
and, more precisely, by consecutive integers 0, ..., *N* − 1. In other words, we
have *V* = {0, ..., *N* − 1}.

Program 67 (see following page) gives a minimal signature for such graphs.
The type *t* is that of graphs. It is deliberately made abstract as we are going to
present two different implementations of this signature. The creation function,
*create*, takes the value *N* as argument, and the function *nbVertex* gives it
back. Crucially, we assume here that the graph contains *N* vertices from the
outset. Consequently, we do not provide operations to add or remove vertices, or
to test whether they belong to the graph. In case of edges, however, we do have
these operations, namely, *addEdge*, *removeEdge*, and *memEdge*. Note that

— openly licensed via CC BY SA 4.0 —
Program 67 — Minimal signature for graphs with integer vertices.

```pascal
type vertex = int
type t

val create : int -> t
val nbVertex : t -> int

val memEdge : t -> vertex -> vertex -> bool
val addEdge : t -> vertex -> vertex -> unit
val removeEdge : t -> vertex -> vertex -> unit

val iterSucc : (vertex -> unit) -> t -> vertex -> unit
val iterEdge : (vertex -> vertex -> unit) -> t -> unit
```

the functions to add and remove edges correspond to an imperative structure. The signature provides only the functions to traverse the successors of a vertex (*iter*) and all edges (*iterEdge*). Traversing all vertices can be easily achieved with a **for** loop from 0 to *N* − 1. Finally, note that this signature works equally well for both directed and undirected graphs.

The signature is deliberately minimal. In particular, graphs are not labeled. Some of the exercises of this section propose extensions.

A straightforward representation of such a graph is a matrix *M* of size *N* × *N*, where each element *M*<sub>*i*,*j*</sub> indicates the presence of an edge between the vertices *i* and *j*. As the graphs are assumed to be unlabeled, it suffices to use a matrix of booleans:

```pascal
type vertex = int
type t = bool array array
```

The definition of the function that creates a graph without any edges is clear: We construct a square matrix, initialized with `false`.

```pascal
let create n = Array.make_matrix n n false
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 68 — Graphs via adjacency matrices.

```pascal
type vertex = int
type t = bool array array

let create n = Array.make_matrix n n false
let nbVertex = Array.length

let mem_edge g v1 v2 = g.(v1).(v2)
let add_edge g v1 v2 = g.(v1).(v2) <- true
let remove_edge g v1 v2 = g.(v1).(v2) <- false

let iterapi f g v =
  Array.iteri (fun w b -> if b then f w) g.(v)
let iterEdge f g =
  for v = 0 to nbVertex g - 1 do iterapi (f v) g v done
```

```pascal
let nbVertex = Array.length
```

In what follows, we assume that all graphs are directed. Exercise 7.4 indicates the modifications required for undirected graphs. Adding, removing, or testing for the presence of an edge can all be done in constant time:

```pascal
let mem_edge g v1 v2 = g.(v1).(v2)
let add_edge g v1 v2 = g.(v1).(v2) <- true
let remove_edge g v1 v2 = g.(v1).(v2) <- false
```

To apply a function `f`: vertex -> unit to all successors of a vertex v of a graph g, it suffices to traverse the row `g.(v)` of the matrix, for example, with `Array.iteri`. We test whether an edge is present before applying the function `f` to the corresponding vertex.

```pascal
let iterapi f g v =
  Array.iteri (fun w b -> if b then f w) g.(v)
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
To apply a function `f`: `vertex -> vertex -> unit` to all edges of a graph `g`, we begin by traversing the vertices of `g` with a `for` loop. Then, for each vertex `v`, we traverse its outgoing edges using `iter`. The function passed to `iterapi` is obtained by partially applying `f` to `v`.

```pascal
let iter_edge f g =
  for v = 0 to nb vertex g - 1 do iterapi (f v) g v done
```

The complete code for graphs via adjacency matrices is given in program 68. The structure of the adjacency matrix may be adapted to the case of labeled graphs (see exercise 7.3) and/or undirected graphs (see exercise 7.4).

## Complexity

The table in figure 7.4 summarizes the complexity of the main operations on adjacency matrices. In general, the cost is expressed as a function of the number *N* of vertices and the number *E* of edges. In the particular case of adjacency matrices, however, the cost depends only on *N*.

<table><caption>Figure 7.4: Cost of operations on adjacency matrices.</caption><thead><tr><td>operation</td><td>cost</td></tr></thead><tbody><tr><td>memEdge</td><td>O(1)</td></tr><tr><td>addEdge</td><td>O(1)</td></tr><tr><td>iterSucc</td><td>O(N)</td></tr><tr><td>iterVertex</td><td>O(N)</td></tr><tr><td>iterEdge</td><td>O(N<sup>2</sup>)</td></tr></tbody></table>

Of course, we assume that the functions passed as arguments of `iterapi`, `iterVertex`, and `iterEdge` have constant cost. From this table, we deduce that adjacency matrices are well adapted to *dense* graphs, that is, graphs in which the number of edges *E* is of the order of *N*². Indeed, the outdegree of each vertex is then of the order of *N*, so that `iterapi` is optimal. Similarly, `iterEdge` must traverse a number of edges of the order of *N*² and is therefore also optimal.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by/4.0/) —
As for the cost in space, an adjacency matrix clearly occupies quadratic
space, namely ($N + 1)^2$ words in case of a boolean matrix (see section 3.2
Runtime Model). We can reduce this cost by a constant factor by using bit
vectors (see section 4.2 Bit Vectors). Exercise 7.5 proposes such an optimization.
It is important to note that this optimization is only possible if the edges are
not labeled.

## 7.2 Adjacency Lists

In case of sparse graphs, an alternative to adjacency matrices consists in using an
array containing the list of the successors of each vertex. We speak of adjacency
lists. The type of graphs follows directly from this definition.

```c
type vertex = int
type t = vertex list array
```

An empty graph is represented by an array containing only empty lists.

```javascript
let create n = Array.make n []
let nbVertex = Array.length
```

To test whether an edge is present between vertices v1 and v2 in a graph g, we
traverse the adjacency list g.(v1) in search of v2. We do this here with the
function List.mem.

```javascript
let memEdge g v1 v2 =
    List.mem v2 g.(v1)
```

We add an edge between v1 and v2 by adding v2 to the list g.(v1). We only
do this if the edge does not already exist, to avoid unnecessary duplicates in the
adjacency list.

```javascript
let addEdge g v1 v2 =
    if not (memEdge g v1 v2) then g.(v1) <- v2 :: g.(v1)
```

We remove the edge v1 → v2 by removing the possible instance of v2 in the list
g.(v1). We can do this easily using List.filter.

— openly licensed via CC BY SA 4.0 —
```python
let remove_edge g v1 v2 =
    g.(v1) <- List.filter ((<>) v2) g.(v1)
```

Note, nevertheless, that this use of `List.filter` has two defects: On the one hand, it does not take advantage of the absence of duplicates in the list to stop early. On the other, it unnecessarily reconstructs the entire list `g.(v1)` in the case where `v2` does not appear. Exercise 7.8 proposes a remedy for these issues. Similarly, the fact that the vertices are integers could be used to optimize the three preceding operations by keeping the adjacency lists sorted (see exercise 7.9).

The traversal of the successors of a vertex `v` of a graph `g` is simpler than in case of adjacency matrices. It is trivially implemented by applying `List.iter` to the adjacency list `g.(v)`.

```python
let iter f g v = List.iter f g.(v)
```

The code of `iter_edge`, by contrast, is the same as that for adjacency matrices.

```python
let iter_edge f g =
    for v = 0 to nbVertex g - 1 do iterapi (f v) g v done
```

The complete code for graphs via adjacency lists is given in program 69.

## Complexity

The table in figure 7.5 summarizes the complexity of the main operations on adjacency lists. Here, $\delta$ denotes the outdegree of a vertex.

Note that, unlike in case of adjacency matrices, the operations `mem_edge` and `add_edge` are no longer in $O(1)$ but in $O(\delta)$. By contrast, `iterapi` and `iter_edge` now have optimal complexity.

As for the cost in space, adjacency lists have optimal complexity $O(N + E)$. To be precise, the cost is $N + 3E + 1$ words (see section 3.2 Runtime Model).

## 7.3 Adjacency Dictionaries

One clear drawback of the two preceding data structures is the condition that vertices must be represented by integers and, moreover, by consecutive integers.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Program 69 — Graphs via adjacency lists.

```pascal
type vertex = int
type t = vertex list array

let create n = Array.make n []
let nbVertex = Array.length

let mem_edge g v1 v2 =
  List.mem v2 g.(v1)

let add_edge g v1 v2 =
  if not (mem_edge g v1 v2) then g.(v1) <- v2 :: g.(v1)

let remove_edge g v1 v2 =
  g.(v1) <- List.filter ((<> ) v2) g.(v1)

let iter f g v =
  List.iter f g.(v)

let iter_edge f g =
  for v = 0 to nbVertex g - 1 do iter (f v) g v done
```
<table><caption>Figure 7.5: Cost of operations on adjacency lists.</caption><thead><tr><td>operation</td><td>cost</td></tr></thead><tbody><tr><td>mem_edge</td><td>O(δ)</td></tr><tr><td>add_edge</td><td>O(δ)</td></tr><tr><td>iter</td><td>O(δ)</td></tr><tr><td>iter vertex</td><td>O(N)</td></tr><tr><td>iter edge</td><td>O(E)</td></tr></tbody></table>

This representation is adapted for graphs whose sets of vertices do not change. If, on the other hand, we wish to add or remove vertices dynamically, we have to modify the preceding structures. We can, for example, use resizeable arrays to permit the addition of vertices. To remove vertices, we could imagine marking them so as to exclude them from the graph, but this would have an impact on the complexity of the functions `iterapi`, `iterapi`, and `iterapi`.

Ideally, we would like a graph structure whose vertices are not necessarily integers and can be added or removed dynamically. The structure would also combine the respective advantages of adjacency matrices and adjacency lists in time and in space. Such a structure exists: It suffices to combine the structure of a dictionary with that of a set. A graph is thus nothing but a dictionary in which each vertex is associated with the set of its successors. In other words, we retain the idea of adjacency lists, but replace the array with a dictionary, and lists by sets. We will call this structure an *adjacency dictionary*. In this way, vertices are not limited to integers, and may be easily added and removed. As for efficiency, it suffices to use a hash table for the dictionary to obtain optimal complexity. We may also use a hash table to represent the set of successors of a vertex.

We can therefore parametrize our adjacency dictionary structure by a module `H` of hash tables (such as the one in section 5.3 *Hash Tables* or the module `Hashtbl` of the OCaml standard library).

module Graph(H: HashTable) = struct

For clarity, we introduce a module `V` for the dictionary, as a synonym of the module `H`. The type of the vertices is then that of the keys of the module `V`.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
```python
module V = H
type vertex = V.key
```

Similarly, we introduce a module `E` for the sets of adjacent vertices, also a synonym of the module `H`. The idea here is to represent a set of vertices by a hash table, by associating each element of the set with the value `()`.

```python
module E = H
```

The type of graphs is therefore the following:

```python
type t = (unit E.t) V.t
```

If the edges were labeled, it would suffice to associate each successor `w` of `v` with the label of the edge $v \to w$ rather than with the value `()` (see exercise 7.10).

Creating a graph becomes merely a matter of constructing an empty dictionary. The number of vertices is obtained simply as the number of keys of the dictionary.

```python
let create () = V.create ()
let nbVertex g = V.length g
```

We check if a vertex `v` belongs to the graph `g` with `V.mem`, that is, by testing if it is a key in the dictionary. We add `v` by creating a binding associating `v` with a new hash table.

```python
let memVertex g v = V.mem g v
let addVertex g v = V.add g v (E.create())
```

We note that the function `addVertex` assumes that `v` is not yet present in the graph `g`. Otherwise, we would have to first verify the presence of `v` using `memVertex`. The removal of a vertex `v` is more delicate: Not only must its binding be removed from the dictionary, it is also necessary to remove every instance of `v` from the adjacency sets of other vertices (here written as `s`).

```python
let removeVertex g v =
    V.remove g v;
    V.iter (fun _ s -> E.remove s v) g
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
For greater efficiency, we should access only the predecessors of v. Exercise 7.14 proposes such an improvement.

To add, remove, or test for the presence of an edge between two vertices v1 and v2, we begin by retrieving the adjacency set of v1 using `V.find g v1`. We then use respectively the functions `E.replace`, `E.remove`, and `E.mem`.

```hkl
let mem_edge g v1 v2 = E.mem (V.find g v1) v2
```

```hkl
let add_edge g v1 v2 = E.replace (V.find g v1) v2 ()
```

```hkl
let remove_edge g v1 v2 = E.remove (V.find g v1) v2
```

These three functions assume that v1 is already present in the graph, that is, that `V.find g v1` does not fail. As for the function `add_edge`, we use `E.replace` rather than `E.add`, to avoid duplicates in the adjacency set.

To traverse the set of vertices of a graph g, it suffices to traverse the set of keys of the dictionary using `V.iter`. Similarly, we traverse the successors of a vertex v by traversing its adjacency set using `E.iter`.

```hkl
let iterVertex f g = V.iter (fun v _ -> f v) g
```

```hkl
let iterSucc f g v = E.iter (fun w _ -> f w) (V.find g v)
```

Here too, we assume that `iterSucc` is called on an existing vertex v. To traverse all edges, it suffices to compose `iterVertex` and `iterSucc` as we did for adjacency matrices and adjacency lists. However, it would be unnecessarily costly to access the adjacency set of v using `V.find g v`. Instead, it suffices to compose directly `V.iter` and `E.iter`.

```hkl
let iterEdge f g = V.iter (fun v s -> E.iter (f v) s) g
```

The complete code of this graph structure is given in program 70.

For undirected, labeled graphs or graphs with multiple edges between two vertices, the reader is referred to exercises 7.10–7.12.

## Complexity

If we assume that dictionaries and sets are implemented by hash tables, we obtain the performances given in the table of figure 7.6.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 70 — Graphs via adjacency dictionaries.

```python
module Graph(H: HashTable) = struct
    module V = H
    type vertex = V.key

module E = H

type t = (unit E.t) V.t

let create () = V.create ()
let nbVertex g = V.length g

let memVertex g v = V.mem g v
let addVertex g v = V.add g v (E.create())
let removeVertex g v =
    V.remove g v;
    V.iter (fun _ s -> E.remove s v) g

let memEdge g v1 v2 = E.mem (V.find g v1) v2
let addEdge g v1 v2 = E.replace (V.find g v1) v2 ()
let removeEdge g v1 v2 = E.remove (V.find g v1) v2

let iterVertex f g = V.iter (fun v _ -> f v) g
let iterSucc f g v = E.iter (fun w _ -> f w) (V.find g v)
let iterEdge f g = V.iter (fun v s -> E.iter (f v) s) g
end
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
<table><caption>Figure 7.6: Cost of operations on adjacency dictionaries.</caption><thead><tr><th>operation</th><th>cost</th><th>operation</th><th>cost</th><th>operation</th><th>cost</th></tr></thead><tbody><tr><td>mem vertex</td><td>O(1)</td><td>mem edge</td><td>O(1)</td><td>iter succ</td><td>O(δ)</td></tr><tr><td>add vertex</td><td>O(1)</td><td>add edge</td><td>O(1)</td><td>iter vertex</td><td>O(N)</td></tr><tr><td>remove vertex</td><td>O(N)</td><td>remove edge</td><td>O(1)</td><td>iter edge</td><td>O(E)</td></tr></tbody></table>

We observe that the complexities are optimal, except in case of `remove vertex` (see exercise 7.14). As for the cost in space, it is necessary to make an assumption regarding the sizes of the hash tables. If we assume that the number of buckets is at least twice the number of bindings in each hash table then the cost in space of a hash table containing $K$ inputs is $6K + 4$ words. Hence, for a graph containing $N$ vertices and $E$ edges, the total cost is:

$$ (6N + 4) + \sum_{v} (6\delta_{v} + 4) $$

This is $10N + 6E + 4$ words in total.

## Persistent Graphs

Graphs implemented using adjacency dictionaries as presented above are imperative due to the imperative nature of hash tables. We naturally obtain persistent graphs by substituting hash tables with a persistent data structure. For example, we may use any dictionary having the signature `Map.S` of OCaml (see exercise 7.15), in particular all those discussed in chapter 5.

## 7.4 Performance Comparison

We recapitulate here the performance of the different graph structures. The time complexity of the main operations are summarized in the table of figure 7.7. As we have already highlighted above, adjacency dictionaries offer the best possible complexity for each operation.

As for the cost in space, recall the cost in number of words for each graph structure:

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
<table><caption>Figure 7.7: Comparison of different graph data structures.</caption><thead><tr><th>operation</th><th>adjacency matrices</th><th>adjacency lists</th><th>adjacency dictionaries</th></tr></thead><tbody><tr><td>mem vertex</td><td>—</td><td>—</td><td>O(1)</td></tr><tr><td>add vertex</td><td>—</td><td>—</td><td>O(1)</td></tr><tr><td>remove vertex</td><td>—</td><td>—</td><td>O(N)</td></tr><tr><td>mem edge</td><td>O(1)</td><td>O(δ)</td><td>O(1)</td></tr><tr><td>add edge</td><td>O(1)</td><td>O(δ)</td><td>O(1)</td></tr><tr><td>remove edge</td><td>O(1)</td><td>O(δ)</td><td>O(1)</td></tr><tr><td>iter succ</td><td>O(N)</td><td>O(δ)</td><td>O(δ)</td></tr><tr><td>iter vertex</td><td>O(N)</td><td>O(N)</td><td>O(N)</td></tr><tr><td>iter edge</td><td>O(N<sup>2</sup>)</td><td>O(E)</td><td>O(E)</td></tr></tbody></table>

* adjacency matrix: $(N + 1)^2$ words;

* adjacency lists: $N + 3E + 1$ words;

* adjacency dictionary: $10N + 6E + 4$ words.

It is clear that adjacency lists are more economical in terms of memory than adjacency dictionaries. The comparison with adjacency matrices is more difficult. Indeed, it depends on the density of the graph, that is, the ratio between the number of vertices $N$ and the number of edges $E$. For a very dense graph, in which $E$ is close to $N^2$, the adjacency matrix is more economical. In the case of a sparse graph, for example where $E \sim N$, adjacency lists are more economical. In general, lists are preferable to matrices when $E < N^2/3$ and, similarly, dictionaries are preferable to matrices when $E < N^2/6$. Lists are always preferable to dictionaries in terms of space.

Of course, in choosing a structure for graphs, the complexity in space of the different operations must be weighed against their complexity in time.

## 7.5 Exercices

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
## Adjacency Matrices

**7.1** The transposition of a directed graph $G$, denoted by $G^R$, is the graph having the same vertices as $G$ and which contains an edge $x \to y$ if and only if there exists an edge $y \to x$ in $G$. Write a function `reverse: graph -> graph` that calculates the transposition of a given graph.

**7.2** Add an operation `nb_edge: t -> int` that gives the number of edges in constant time. Hint: Keep track of the number of edges in the structure of the graph by updating its value in `add_edge` and `remove_edge`.

**7.3** Modify adjacency matrices for graphs in which edges are labeled by a given type `label`.

**7.4** The easiest way to represent undirected graphs is by preserving the same structure as for directed graphs, but by maintaining the invariant that, for each edge $a \to b$, there is also an edge $b \to a$. Modify the operations `add_edge` and `removeEdge` of adjacency matrices and lists accordingly. Also modify `iterEdge` so that it traverses each edge only once.

**7.5** Modify the adjacency matrix structure to use bit vectors (see section 4.2 *Bit Vectors*) rather than boolean arrays. What is the gain in terms of space?

**7.6** The transitive closure of a graph $G$ is a graph $T$ having the same vertices as $G$, and an edge between $i$ and $j$ when there is a *path* between $i$ and $j$ in $G$. When the graph is represented by a boolean matrix, the transitive closure may be calculated in $O(N^3)$ time using the Floyd-Warshall algorithm. The pseudo-code is as follows:

$$ 
\begin{array}{l}
T \leftarrow G \\
\quad \text{for } k \text{ from } 0 \text{ to } N-1 \\
\quad \quad \text{for } i \text{ from } 0 \text{ to } N-1 \\
\quad \quad \quad \text{for } j \text{ from } 0 \text{ to } N-1 \\
\quad \quad \quad \quad T_{i,j} \leftarrow T_{i,j} \text{ or } (T_{i,k} \text{ and } T_{k,j})
\end{array}
$$

The idea consists in determining, for increasing values of $k$, if there exists a path between $i$ and $j$ that uses only intermediate vertices smaller than $k$. The last line of the pseudo-code considers the two cases of a path between $i$ and $j$, passing via $k$, or not. Write a function `transitive_closure: t -> t` implementing this algorithm.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
## Adjacency Lists

7.7 Redo exercises 7.2 to 7.4 for adjacency lists.

7.8 Modify the function `remove_edge` on adjacency lists so that (1) it does not modify the adjacency list if the edge to be removed does not exist, and (2) it makes use of the fact that adjacency lists do not contain duplicates to stop once the edge is found.

7.9 Modify the operations on adjacency lists to maintain the invariant that adjacency lists are sorted. Make use of this property to optimize some of the operations.

## Adjacency Dictionaries

7.10 Modify the code of program 70 for graphs whose edges are labeled by values of any type. More precisely, give graphs the polymorphic type `type 'a t = ('a E.t) V.t`, where 'a is the type of labels.

7.11 Modify the adjacency dictionary structure for undirected graphs drawing on exercise 7.4.

7.12 Modify the adjacency dictionary structure for graphs with *multiple edges*, that is, where two vertices may be linked by more than one edge. For unlabeled multi-edges, the dictionary must associate each vertex with the *multi-set* of its successors. For labeled multi-edges, we may distinguish two cases, depending on whether the edges between $v_1$ and $v_2$ can bear the same label or not.

7.13 Modify the functor of program 70 so that, instead of taking a module `H` of hash tables as argument, it takes:

* either a module `V` introducing a type `t` equipped with hash and comparison functions;

* or two modules, `M` and `S`, that respectively implement dictionaries and sets on the same type of vertices.

7.14 Modify the adjacency dictionary structure to keep track of the set of predecessors of a vertex in an efficient manner. In other words, provide an operation

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
```hkl
val iter_pred: (vertex -> unit) -> t -> vertex -> unit
```

whose cost is proportionate to the indegree. Use this function to reimplement
the function `removeVertex`.

**7.15** Modify the functor of program 70 to obtain persistent graphs by replac-
ing the argument `H` by a module `M` of signature `Map.S`.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
# Disjoint Sets

This chapter presents an imperative data structure that solves the problem of partitioning a finite set into disjoint subsets, called “classes.” We wish to determine whether two elements belong to the same class and merge two classes into one. It is these two operations that give the structure its name: *union-find*.

## 8.1 The Basic Idea

Without loss of generality, we may assume that the set to be partitioned consists of the $n$ integers $\{0, 1, \dots, n-1\}$. Program 71 (see page 315) gives the signature of such a structure. The operation `create n` constructs a new partition of $\{0, 1, \dots, n-1\}$ in which each element forms a class by itself. The operation `find` determines the class of an element by returning one of the integers as a distinguished representative. In particular, we can determine if two elements are in the same class by comparing the result given by `find` for each one. Finally, the operation `union` merges two classes of the partition, modifying the data structure in place.

The basic idea is to link the elements of a class between themselves. In each class, these links form a graph in which all paths lead to the representative,

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
<figure><img src="image_106.png" /><figcaption>Figure 8.1: A partition of {0, 1, ..., 7} into two sets.</figcaption></figure>

which is the only element linked to itself. Figure 8.1 illustrates a situation in which the set {0, 1, ..., 7} is partitioned into two classes, whose representatives are respectively 3 and 4.

It is possible to represent such a structure using individually allocated nodes (see exercise 8.5). Nevertheless it is simpler, and often more efficient, to use an array that links each integer to another integer of the same class. These links lead to the class representative, which is associated with itself in the array. Thus, the partition of figure 8.1 is represented by the following array:

<table><tr><td>0</td><td>1</td><td>...</td><td>7</td></tr><tr><td>7</td><td>5</td><td>6</td><td>3</td><td>4</td><td>3</td><td>4</td><td>3</td></tr></table>

The **find** operation simply follows the links until the representative is found. The **union** operation begins by finding the representatives of both elements and then links them together. We propose two improvements in the interest of efficiency. The first involves performing *path compression* during the search carried out by **find**: All the elements found along the path traversed to reach the representative are linked directly to it. The second consists in keeping track of the *rank* of each representative, that is, the maximum length of a path in its class. This information is stored in a second array and is used by the function **union** to choose the representative for the union of the classes.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
<caption>Program 71 — Signature of the union-find structure.</caption>

```pascal
module type UnionFind = sig
type t
val create : int -> t
val find : t -> int -> int
val union : t -> int -> int -> unit
end
```

## 8.2 Implementation

We now describe the code of the union-find structure. The type `t` is a record containing two arrays: `rank`, containing the rank of each class; and `link`, containing the links between elements.

```pascal
type t = {
    rank: int array;
    link: int array;
}
```

The information contained in `rank` is only meaningful for those elements `i` that are representatives, that is, for which `link.(i) = i`. Initially, each element forms a class in itself (that is, it is its own representative), and the rank of every class is 0.

```pascal
let create n =
    { rank = Array.make n 0;
      link = Array.init n (fun i -> i) }
```

The function `find` computes the representative of an element `i`. It is naturally written as a recursive function. We begin by calculating the element `p` linked to `i` in the array `t.link`. If it is `i` itself, `i` is the class representative, and we are done.

```pascal
let rec find t i =
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 72 — Union-find structure.

```pascal
type t = {
    rank: int array;
    link: int array;
}

let create n =
    { rank = Array.make n 0;
      link = Array.init n (fun i -> i) }

let rec find t i =
    let p = t.link.(i) in
    if p = i then
        i
    else begin
        let r = find t p in
        t.link.(i) <- r;
        r
    end

let union t i j =
    let ri = find t i in
    let rj = find t j in
    if ri <> rj then begin
        if trank.(ri) < trank.(rj) then
            t.link.(ri) <- rj
        else begin
            t.link.(rj) <- ri;
            if trank.(ri) = trank.(rj) then
                t(rank.(ri) <- t(rank.(ri) + 1)
        end
    end
```
```pascal
let p = t.link.(i) in
if p = i then
  i
```

Otherwise, we recursively compute the representative `r` of `p` with `find t p`. However, before returning `r`, we perform path compression, that is, we link `i` to `r` directly.

```pascal
else begin
  let r = find t p in
t.link.(i) <- r;
  r
end
```

This way, the next time we call `find` on `i`, we will immediately find `r`.

The operation `union` merges the classes of two elements, `i` and `j`. We begin by calculating their representatives, `ri` and `rj`, respectively. If they are equal, there is nothing more to be done.

```pascal
let union t i j =
  let ri = find t i in
  let rj = find t j in
  if ri <> rj then begin
```

Otherwise, we compare the ranks of the two classes. If that of `ri` is strictly less than that of `rj`, we make `rj` the representative of the union.

```pascal
  if trank.(ri) < trank.(rj) then
    t.link.(ri) <- rj
```

The rank does not have to be updated for this new class. Indeed, only the paths in the previous class of `ri` have had their lengths increased by one unit, and this new length does not exceed the rank of `rj`. When the rank of `rj` is smaller, we proceed symmetrically.

```pascal
else begin
  t.link.(rj) <- ri;
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
When the two classes have the same rank, we choose one of the two representatives arbitrarily as the representative of the union (here, $ri$). The rank information must then be updated, because in this case the length of the longest path in the class may increase by one unit.

```pascal
if trank.(ri) = trank.(rj) then
trank.(ri) <- trank.(ri) + 1
end
end
```

Importantly the function **union** uses the function **find**, and therefore performs path compression even in the case when i and j belong to the same class. The complete code of the *union-find* structure is given in program 72 (see page 316).

## Complexity

We can show that, thanks to path compression and the rank associated with each class, a sequence of $m$ **find** and **union** operations, performed on a structure containing $n$ elements, is executed in $O(m\alpha(n, m))$ total time, where $\alpha$ is a function that grows extremely slowly. In fact, it grows so slowly that we may consider it to be constant for all practical purposes, given the constraints that memory and time impose on the values of $n$ and $m$. We may therefore assume that each operation takes amortized constant time. This analysis is complex and lies beyond the scope of this book. A detailed explanation may be found in *Introduction to Algorithms* [7, chap. 22].

<figure><img src="image_107.png" /></figure>

**For Further Information**

The data structure presented in this chapter is attributed to McIlroy and Morris [3], and its complexity has been analyzed by Tarjan [23].
## 8.3 Exercises

8.1 The signature of the *union-find* structure is constraining because it requires prior knowledge of the number of elements. For greater flexibility, we may wish to add new elements dynamically, and without any restriction on the number of elements. The signature of such a structure could be the following:

```h
type t
val create: unit -> t
val add: t -> int
val find: t -> int -> int
val union: t -> int -> int -> unit
```

The function `create` creates an empty structure. The function `add` adds a new element and returns it. (Elements are always integers). The signatures of `find` and `union` remain unchanged. Implement such a data structure using resizeable arrays (see section 4.1 *Resizeable Arrays*).

8.2 Extend the *union-find* structure with an operation `num_classes` that returns the total number of classes.

```h
val num_classes: t -> int
```

Ensure that `num_classes` executes in constant time by keeping track of this value as an additional field of the record.

8.3 Extend the *union-find* structure with an operation `iter_classes` that permits the traversal of the set of all class representatives.

```h
val iter_classes: (int -> unit) -> t -> unit
```

Ensure that `iter_classes` executes in time proportionate to the total number of classes.

8.4 If the elements are not consecutive integers, we can replace the two arrays `rank` and `link` by a pair of hash tables. Reimplement the operations `create`, `find`, and `union` using this idea.

8.5 Instead of using arrays, another way of implementing the *union-find* structure consists in representing each class directly as an acyclic graph whose nodes

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
are records containing the values `rank` and `link`. Each node is of the following type:

```hkl
type elt = { mutable rank: int; mutable link: elt }
```

If we wish to keep track of extra information for each element, we may use the following type instead:

```hkl
type 'a elt =
    { mutable rank: int; mutable link: 'a elt; data: 'a }
```

It is no longer necessary to keep track of any global information in the `union-find` structure, because each node contains all the necessary information. The signature of the `union-find` structure is then modified as follows:

```hkl
type 'a elt
val create_node: 'a -> 'a elt
val find: 'a elt -> 'a elt
val union: 'a elt -> 'a elt -> unit
```

The function `create_node` constructs a class consisting of a single element, that is, a record whose field `link` points to itself. Write the operations `make`, `find`, and `union`. Note: The comparison of values of type `'a elt` must be done using physical equality.

8.6 Replace the arrays `rank` and `link` with persistent arrays (see section 4.4 *Persistent Arrays*), and write a persistent `union-find` structure with the following signature:

```hkl
type t
val create: int -> t
val find: t -> int -> int
val union: t -> int -> int -> t
```

Note: To perform path compression, it is convenient to modify the contents of the field `link` using side effects. We accordingly define the following type:

```hkl
type t = { rank: int A.t; mutable link: int A.t }
```

Here, the type `A.t` is the type of persistent arrays, and the mutability of the field `link` can be used to record the effect of path compression before returning the

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
result of the function `find`. Although the data structure is modified by a side effect, the modification cannot be observed: the representative of each element remains the same.

**8.7** The *union-find* structure can be used to efficiently construct a perfect maze, that is, one in which there is one and only one path between any two cells. Here is an example of such a maze:

<figure><img src="image_108.png" /></figure>

You may proceed as follows: Create a *union-find* structure whose elements are the different cells. The idea is that two cells are in the same class if and only if they are linked by a path. Initially, all the cells of the maze are separated from each other by walls. Next, consider all pairs of adjacent cells (vertically and horizontally) in random order. For each pair $(c_1, c_2)$, compare the classes of the cells $c_1$ and $c_2$. If they are identical, there is nothing more to be done. Otherwise, delete the wall that separates $c_1$ and $c_2$, and merge the two classes using *union*. Write code that constructs a maze using this method.

**Hint:** The easiest way to traverse all pairs of adjacent cells in random order is by constructing an array containing all pairs, and then mixing them randomly using the *Knuth shuffle* (exercise 2.11, page 134).

Justify that at the end of the construction every cell is linked to every other cell by a single path.

**8.8** The aim of this exercise is to color the *connected* components of a black-and-white image using different colors. A connected component is formed by pixels that have the same color as their neighbors. Two pixels are neighbors if they have a common border.

For example, in the following image of 6 × 6 pixels, there are two white connected components, and four black ones.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
<figure><img src="image_109.png" /></figure>

A black-and-white image is represented by a matrix of booleans. To compute its connected components, we use a *union-find* structure to group neighboring pixels having the same color in a single class. Write a function that takes as argument a matrix of booleans and returns a *union-find* structure containing its equivalence classes.

To assign a different color to each connected component, use a hash table that associates a different color with each representative. Write a function that displays the image with its connected components colored. Assume that you are given a function `new_color` of type `unit -> color` and a function `draw pixel` of type `int -> int -> color -> unit`.
9

Zippers

In this chapter, we present a general technique to traverse data structures. This technique does not involve a specific traversal order and allows for the possibility of local modifications (inserting, removing, etc.). It goes by the name of *zipper*. We present it here in the context of lists and trees, but it can be adapted to many other data structures.

9.1 Zippers on Lists

Suppose we wish to “navigate” within a list of type 'a list, that is, move from element to element, occasionally performing modifications. The image that comes to mind here is that of a cursor in a text editor, and the action of the keys of the keyboard, used to move, insert, delete characters, etc.

Consider a list containing seven elements, with the cursor placed after the
third, as illustrated in figure 9.1.

We represent this situation with two lists, one containing the elements situated to the left of the cursor, and another containing the elements situated to the right. It is convenient to represent the first of these lists in reverse order, so that the element situated immediately to the left of the cursor is easily

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by/4.0/) —
<figure><img src="image_110.png" /><figcaption>Figure 9.1: A cursor on a list.</figcaption></figure>

accessible. Figure 9.2 illustrates this representation on the preceding example.

<figure><img src="image_111.png" /><figcaption>Figure 9.2: A zipper on a list.</figcaption></figure>

This structure consisting of two lists is called a *zipper*. We define it using the following type:

```h
type 'a zipper = { left: 'a list; right: 'a list; }
```

To construct such a zipper from a list `l`, with the cursor placed all the way to the left, we initialize the left list to `[]` and the right list to `l`.

```h
let of_list l =
    { left = []; right = l }
```

Let us consider next the navigation and modification operations for the zipper.

## Navigation Operations

The first navigation operation consists in moving the cursor to the right. This is only possible if the right list contains at least one element. In this case, we move the first element of this list to the head of the left list.

```h
let move_right z = match z.right with
| [] -> invalid_arg "move_right"
| x :: r -> { left = x :: z.left; right = r }
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
The application of the function `move_right` to the example of figure 9.2 returns the zipper illustrated in figure 9.3.

<figure><img src="image_112.png" /><figcaption>Figure 9.3: The zipper after a move to the right.</figcaption></figure>

We similarly write a function `move_left` to move the cursor to the left (see program 73, following page).

To convert a zipper `z` into a list, it suffices to concatenate the two lists `z.left` and `z.right`, reversing the former. This is precisely what the function `List.rev_append` does.

```hkl
let to_list z =
    List.rev_append z.left z.right
```

Equivalently, we could have moved the cursor all the way to the left and then retrieved the right list.

## Modification Operations

To insert an element at the cursor position, it suffices to place that element at the head of one of the two lists. If we insert it at the head of the left list, the cursor will be to the right of the new element (as in case of a text editor).

```hkl
let insert z x =
    { z with left = x :: z.left }
```

If we wish instead to leave the cursor to the left of the inserted element, it suffices to replace `left` by `right` in the code above.

Another modification operation consists in removing an element. For example, we may remove the element to the left of the cursor, if such an element exists. (This is what the `backspace` key does in a text editor.) This operation is equivalent to removing the head of the `left` list.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 73 — Structure of a zipper on a list.

```h
type 'a zipper = { left: 'a list; right: 'a list; }

let of_list l =
    { left = []; right = l }

let move_right z = match z.right with
| [] -> invalid_arg "move_right"
| x :: r -> { left = x :: z.left; right = r }

let move_left z = match z.left with
| [] -> invalid_arg "move_left"
| x :: l -> { left = l; right = x :: z.right }

let to_list z =
    List.rev_append z.left z.right

let insert z x =
    { z with left = x :: z.left }

let delete_left z = match z.left with
| [] -> invalid_arg "delete_left"
| _ :: l -> { z with left = l }

let delete_right z = match z.right with
| [] -> invalid_arg "delete_right"
| _ :: r -> { z with right = r }
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
```h
let delete_left z = match z.left with
| [] -> invalid_arg "delete_left"
| _ :: l -> { z with left = l }
```

We can similarly write an operation `delete_right` that removes the element immediately to the right of the cursor.

## 9.2 Zippers on Trees

The technique of the zipper also applies to other data structures, such as trees. Consider, for example, the case of polymorphic binary trees defined by the following type:

```h
type 'a tree = E | N of 'a tree * 'a * 'a tree
```

What is the analog of the cursor in case of a tree? It is the designation of a particular node together with the possibility of navigating within the tree, that is, moving up, or moving down along the left or right subtrees.

Consider the tree in figure 9.4, in which we place the cursor on the node containing 5. The zipper contains the position of the node in question and the subtree at that position. The former is represented by the path from the root to the node, that is, the sequence of moves indicating at each step whether we descend along the left or the right subtree. For example, the path to reach node 5 is the sequence `[left; right]`.

<figure><img src="image_113.png" /><figcaption>Figure 9.4: A cursor on a tree.</figcaption></figure>

We introduce the following type to represent such paths. The constructor `Top` represents the empty path, that is, the position of the root of the tree. The

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
constructors `Left` and `Right` indicate, respectively, a move along the left or right subtrees.

<figure><img src="image_114.png" /></figure>

For each move, we keep track of the subtree along which we did not descend. This is necessary in order to be able to move back up the tree. The zipper on binary trees is then defined by the following type:

<figure><img src="image_115.png" /></figure>

The field `tree` contains the subtree at the cursor. The field `path` contains the path leading to this subtree. Just as in case of the zipper on lists, where we chose to represent the left list in reverse order, we choose here to store the path of the node *towards* the root, rather than *from* the root. Figure 9.5 gives the zipper for the preceding example.

<figure><img src="image_116.png" /><figcaption>Figure 9.5: A zipper on a tree.</figcaption></figure>

The first constructor of the path is `Right` because the cursor designates a right subtree. The next step of the path is `Left` because node 3 is the root of a left subtree. Finally, the path reaches the root of the tree, 8.

The zipper with the cursor at the root of a tree is constructed using the empty path.

<figure><img src="image_117.png" /></figure>

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Let us now write the navigation functions. We begin with a function `down_left` to move down the left subtree. For this, we examine the form of the tree at the cursor. We fail if it is empty.

```hml
let down_left z = match z.tree with
  | E -> invalid_arg "down_left"
```

Otherwise, we return the zipper where the subtree at the cursor is the left subtree, and the path is extended with the constructor `Left` to signify that we have moved down the left.

```hml
| N (l, x, r) -> { path = Left (z.path, x, r); tree = l }
```

The node `x` and its right subtree `r` are kept in the path, as arguments of the constructor `Left`. We similarly write a function `down_right` that moves down the right subtree.

The function to move up the tree examines the current path. We fail if the path is empty.

```hml
let up z = match z.path with
  | Top -> invalid_arg "up"
```

Otherwise, we must reconstruct the node of the tree situated above the cursor. For this, we use the information contained in the path as well as the subtree contained in the field `tree`. If we move down the left subtree, the path would be of the form `Left (p, x, r)`, and the reconstructed node would have `z.tree` as its left subtree, `x` as its root, and `r` as its right subtree.

```hml
| Left (p, x, r) ->
  { path = p; tree = N (z.tree, x, r) }
```

If we move down the right subtree, we proceed symmetrically.

```hml
| Right (l, x, p) ->
  { path = p; tree = N (l, x, z.tree) }
```

The function that reconstructs the entire tree from a given zipper moves the cursor up the tree until it reaches the top. We stop when the path is `Top` and return the tree contained in the field `tree`.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 74 — A zipper structure on a binary tree.

```h
type 'a tree = E | N of 'a tree * 'a * 'a tree

type 'a path =
  | Top
  | Left of 'a path * 'a * 'a tree
  | Right of 'a tree * 'a * 'a path

type 'a zipper = { path: 'a path; tree: 'a tree }

let of_tree t = { path = Top; tree = t }

let down_left z = match z.tree with
  | E -> invalid_arg "down_left"
  | N (l, x, r) -> { path = Left (z.path, x, r); tree = l }

let down_right z = match z.tree with
  | E -> invalid_arg "down_right"
  | N (l, x, r) -> { path = Right (l, x, z.path); tree = r }

let up z = match z.path with
  | Top -> invalid_arg "up"
  | Left (p, x, r) -> { path = p; tree = N (z.tree, x, r) }
  | Right (l, x, p) -> { path = p; tree = N (l, x, z.tree) }

let rec to_tree z =
  if z.path = Top then z.tree else to_tree (up z)
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
```h
let rec to_tree z =
    if z.path = Top then z.tree else to_tree (up z)
```

The complete code of the zipper on binary trees is given in program 74.

Note: You may wonder about the commonalities between zippers on lists and those on trees. Indeed, the zipper on lists seems to be a completely symmetric structure, its two fields `left` and `right` being of the same type. However, you may also understand it as being asymmetrical, with `left` playing the same role as `path`, and `right`, that of `tree`. Thus, we could have defined the following type to represent a position in a list:

```h
type 'a path = Top | Right of 'a * 'a path
```

Since this type is clearly identical to `'a list`, we did not consider it worthwhile to define a new type. Zippers on lists and trees thus embody the very same idea.

### Application: Comparing Two Binary Search Trees

In chapter 5, we presented binary search trees, which can represent sets of elements equipped with a total order. If we wish to construct a set of sets, it is therefore necessary to equip binary search trees themselves with a total order.

One solution is to construct the list of elements for the two trees in infix order and then compare the two lists, for example, lexicographically. This is nevertheless somewhat naive, since the two trees may contain many elements but may differ rapidly, in which case we would have constructed the two lists in vain.

A more efficient solution is preferable, working directly on the trees. However, this is not straightforward since two binary search trees can contain the same elements without having the same structure, as illustrated in figure 9.6. One solution consists in using the zipper structure to execute a *simultaneous* infix traversal of the two trees.

We begin by writing a function that descends all the way to the bottom left of a tree and returns the zipper corresponding to this position. This amounts to iterating the function `down_left` as long as possible.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
<figure><img src="image_118.png" /><figcaption>Figure 9.6: Two binary search trees containing the same elements.</figcaption></figure>

Concretely, we write a function `leftmost` by passing it a zipper as accumulator. If the tree is empty, we return the zipper passed as argument:

```h
let rec leftmost z = function
  | E -> z
```

Otherwise, we move down the left subtree, accumulating the visited node in the zipper.

```h
| N (l, x, r) -> leftmost (Left (z, x, r)) l
```

If we consider the two trees of figure 9.6, `t1` and `t2`, then the zippers obtained with `leftmost` Top `t1` and `leftmost` Top `t2` have the following form:

```h
leftmost Top t1 = Left (Left (..., 3, ...), 1, ...)
leftmost Top t2 = Left (Left (..., 4, ...), 1, ...)
```

Figure 9.7 depicts the two zippers. We can verify that both start off with the smallest value, namely, 1.

Next, we have to write a function `compare` to perform the actual comparison. This function takes as argument two zippers and a function `cmp` to compare the elements. If the two zippers are `Top`, the comparison is over, and we return 0 to indicate equality.

```h
let rec compare cmp z1 z2 = match z1, z2 with
  | Top, Top ->
    0
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
<figure><img src="image_119.png" /><figcaption>Figure 9.7: The two zippers for the trees of figure 9.6.</figcaption></figure>

If the two zippers are of the form `Left`, we compare the elements that they designate using the function `cmp`. If they differ, the comparison is over, and the result given by `cmp` is returned.

| Left (z1, x1, r1), Left (z2, x2, r2) ->
  let c = cmp x1 x2 in
  if c <> 0 then c

If they are equal, however, the comparison must continue. We look for the following elements, by calling `leftmost` on the two right subtrees, `r1` and `r2`. We then call `compare` recursively.

else compare cmp (leftmost z1 r1) (leftmost z2 r2)

If one of the two zippers equals `Top` and the other `Left`, this means that one of the two traversals has reached the end, but not the other. We therefore return `-1` or `1`, as the case may be.

| Top, Left _ ->
  -1
| Left_, Top ->
  1

Finally, by construction, the case of a zipper of the form `Right` cannot occur. We dispense with it as follows:

| Right_, _ |_, Right _ ->
  assert false

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 75 — Lexicographic comparison of binary trees.

```pascal
let rec leftmost z = function
  | E -> z
  | N (l, x, r) -> leftmost (Left (z, x, r)) l

let rec compare cmp z1 z2 = match z1, z2 with
  | Top, Top ->
    0
  | Left (z1, x1, r1), Left (z2, x2, r2) ->
    let c = cmp x1 x2 in
    if c <> 0 then c
    else compare cmp (leftmost z1 r1) (leftmost z2 r2)
  | Top, Left _ ->
    -1
  | Left_, Top ->
    1
  | Right_, _ |_, Right _ ->
    assert false

let compare_tree cmp t1 t2 =
  compare cmp (leftmost Top t1) (leftmost Top t2)
```

To solve the initial problem of how to compare two trees t1 and t2, it suffices to construct the two zippers with the function `leftmost`, and then compare them using the function `compare`.

```pascal
let compare_tree cmp t1 t2 =
  compare cmp (leftmost Top t1) (leftmost Top t2)
```

The complete code is given in program 75. Note that the two functions `leftmost` and `compare` are tail recursive.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Program 76 — Cursor structure.

```pascal
type 'a enum
val start: 'a t -> 'a enum
val step: 'a enum -> 'a * 'a enum (* raise Exit when done *)
```

## 9.3 Cursors

When discussing the problem of comparing binary trees, we saw the importance of traversing a data structure step by step. This method of traversal is different from that of higher-order iterators like `iter` and `fold`, with which we would not have been able to solve that problem.

In general, for an arbitrary data structure of type `'a t'`, we would like a `cursor` interface in order to traverse its elements one by one. Such an interface is given in program 76.

The interface contains an abstract polymorphic type, `'a enum'`, for the cursor, where the variable `'a` is the type of the elements. The function `start` returns a new cursor. The function `step` returns the current element and the cursor corresponding to the next element, or raises the exception `Exit` if the iteration is complete.

This interface defines persistent cursors, which means that the cursor passed as input to the function `step` is not modified. We may thus reuse the cursor, for example, in an algorithm that uses *backtracking*. Another possibility is that of imperative cursors, where the function `step` returns only the current element and advances the cursor to the following element by a side effect (see exercise 9.12).

We will now implement persistent cursors on lists and trees.

### Cursors on Lists

The zipper on lists is itself a convenient cursor. It is, nevertheless, unnecessarily complex since we do not need to keep track of the elements situated to the left

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
## Program 77 — Cursors on lists.

```pascal
type 'a enum = 'a list

let start l =
    l

let step = function
    [ ] -> raise Exit
    x :: r -> x, r
```

of the cursor. In other words, the cursor for lists is nothing but the right list, that is:

```pascal
type 'a enum = 'a list
```

The function `start` reduces to the identity. The function `step` raises the exception `Exit` if the cursor is empty. Otherwise, it returns the head of the list and the rest of the list as the new cursor.

```pascal
let step = function
    [ ] -> raise Exit
    x :: r -> x, r
```

Note that the persistent nature of lists ensures the persistence of the cursor. The complete code is given in program 77.

## Cursors on Trees

To define a cursor on trees, it is first necessary to choose a traversal order. We arbitrarily choose infix order here. (See exercises 9.8 et 9.9 for other traversal orders.)

Unlike zippers on lists, zippers on trees do not immediately provide us with an operation to produce the next element in the infix traversal. However, we have seen how the zipper permits this operation when dealing with the problem

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by/4.0/) —
Program 78 — Cursors on trees (infix traversal).

```pascal
type 'a enum = Top | Left of 'a * 'a tree * 'a enum

let rec leftmost t e = match t with
  | E -> e
  | N (l, x, r) -> leftmost l (Left (x, r, e))

let start t =
  leftmost t Top

let step = function
  | Top -> raise Exit
  | Left (x, r, e) -> x, leftmost r e
```

of comparing binary trees. Only the two constructors `Top` and `Left` of the zipper are necessary. We therefore define the following type for the cursor:

```pascal
type 'a enum = Top | Left of 'a * 'a tree * 'a enum
```

The function `start` places the cursor on the element situated at the bottom left of the tree. For this, we use the function `leftmost` that was written for program 75.

```pascal
let start t =
  leftmost t Top
```

The function `step` raises the exception `Exit` if the cursor is `Top`. Otherwise, it returns the element `x` contained in the constructor `Left` and constructs a new cursor by calling `leftmost` on the right subtree of `x`.

```pascal
let step = function
  | Top -> raise Exit
  | Left (x, r, e) -> x, leftmost r e
```

This is exactly what we did in program 75 to move to the next element.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
As in case of lists, the immutable nature of the type `enum` guarantees the persistence of the cursor. The complete code is given in program 78 (see previous page).

<figure><img src="image_120.png" /></figure>

**For Further Information**

Zippers were invented by Gérard Huet and presented in the article *The Zipper* [12]. As the author noted, this concept was surely already known to other programmers.

The OCaml standard library provides a module `Seq` that implements delayed lists. This is an alternative to cursors that may be used to traverse data structures step by step.

## 9.4 Exercises

### Zippers on Lists

9.1 Equip the zipper on lists with functions `to_start` and `to_end`, that permit moving to the beginning and end of a list, respectively.

9.2 Rewrite the function `to_list` in the same manner as the function `to_tree`.

9.3 Consider the following problem: Given a non-empty list of integers, determine if there exists a binary tree whose leaves, considered in infix order, are located at the depths given by the list. For example, the list `[1; 3; 3; 2]` corresponds to a binary tree, but not the list `[1; 3; 3]`. We propose the following algorithm to solve this problem:

* If the list contains a single element, the algorithm terminates. The result is positive if and only if the element is 0, which corresponds to the empty tree.

* Otherwise, we search for the first two consecutive elements of the list that are equal. If there are none, we fail. Otherwise, calling `n` the common

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
value of the two integers, we replace the two occurrences of $n$ by the integer $n-1$, which corresponds to an internal node.

* Start again with the resulting list.

Write a function `is_tree: int list -> bool` that implements this algorithm by using the zipper structure so as to avoid starting at the beginning of the list each time. The total complexity must be linear in the length of the list.

9.4 It is possible to program an efficient text editor using the zipper structure. Each line of text is a zipper on a list of characters, as explained earlier. Similarly, a file is viewed as a zipper on the list of lines. We use the following types to represent a line and a text, respectively:

`type line = { left: char list; right: char list }`

`type text = { up: line list; current: line; down: line list }`

The field `current` of the type `text` represents the current line on which the cursor is found. The fields `up` and `down` are the analogs of `left` and `right` for lines. Write the functions `insert_char`, `return`, `backspace`, etc. for such a text editor.

## Zippers on Trees

9.5 Write a function `remove_leaf: 'a zipper -> 'a zipper` that verifies that the zipper passed as argument points to a node that has no children, and removes it. Such a node is of the form `N(E,x,E)`. Conversely, write a function `insert_leaf: 'a -> 'a zipper -> 'a zipper` that verifies that the zipper passed as argument points to a leaf node (that is, `E`), and which inserts in its place a new node containing the value passed as argument.

9.6 Write a function `remove_leftmost: 'a zipper -> 'a zipper` that verifies that the zipper passed as argument points to a node of the tree that has no left child, and removes it. Make sure to retain the right subtree, if it is present.

9.7 Consider $n$-ary trees labeled by strings, defined by the following type:

`type tree = N of string * tree list`

Define a zipper structure on such trees.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
## Cursors

9.8 Write a cursor on binary trees corresponding to a prefix traversal, that is, a traversal in which each node is visited before its two subtrees.

9.9 Write a cursor on binary trees corresponding to a postfix traversal, that is, a traversal in which each node is visited after its two subtrees.

9.10 Rewrite the comparison of binary trees (program 75) using the cursor on trees (program 78).

9.11 Rather than raising an exception `Exit` when the cursor has reached the end of the traversal, the function `step` can return a value of type `option`, that is, `step: 'a enum -> ('a * 'a enum) option`. Rewrite the cursors on lists and trees with this new interface.

9.12 If the persistence of the cursor is not used, we may consider an imperative version in which the function `step` returns the current element only and advances to the next element by a side effect, that is, `step: 'a enum -> 'a`. Rewrite the cursors on lists and trees with this new interface.

9.13 Implement a cursor on arrays.
Part III

# Algorithmic Techniques and Applications

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by/4.0/) —
<figure><img src="image_121.png" /></figure>
## 10.1 Euclid's Algorithm

The most famous of all algorithms is certainly that of Euclid. It calculates the greatest common divisor of two integers, $x$ and $y$, written $x \land y$, and called the "gcd of $x$ and $y$". Given two non-negative integers $x$ and $y$, Euclid's algorithm replaces the pair $(x, y)$ by the pair $(y, x \text{ mod } y)$, continuing until $y$ reaches zero. It then returns the value of $x$, which is the gcd of the initial values of $x$ and $y$. We write this algorithm in the form of a recursive function gcd:

```pascal
let rec gcd x y =
    if y = 0 then x else gcd y (x mod y)
```

The correctness of the algorithm follows from the fact that each iteration preserves the greatest common divisor of $x$ and $y$. When we reach $y = 0$, we return $x$, that is, $x \land 0$, which is then the gcd of the initial values of $x$ and $y$. In other words, we have used the following two properties of the gcd:

$$x \land y = y \land (x \text{ mod } y),$$

$$x \land 0 = x.$$

The termination of the algorithm follows from the fact that $y$ decreases

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
strictly in each iteration, remaining non-negative throughout. Indeed, we have the invariant that $x, y \ge 0$.

The complexity of Euclid's algorithm is given by Lamé's theorem. It states that, if the algorithm performs $s$ iterations with $x > y > 0$, then $x \ge F_{s+1}$ and $y \ge F_s$, where ($F_n$) is the Fibonacci sequence. We deduce from this that Euclid's algorithm is logarithmic, that is, $s = O(\log x)$. In the general case, the complexity is $O(\log(\max(x, y)))$; see exercise 10.1. A detailed analysis may be found in *The Art of Computer Programming* [14, sec. 4.5.3].

Exercise 10.2 generalizes the function `gcd` to integers $x$ and $y$ that are not necessarily non-negative.

## The Extended Euclidean Algorithm

We can easily modify Euclid's algorithm to compute the Bézout coefficients, that is, two integers $u$ and $v$ satisfying the following equation:

$$ux + vy = x \wedge y. \tag{10.1}$$

This is known as the extended Euclidean algorithm. We implement it in the form of a function `extended_gcd` that takes the integers $x$ and $y$ as arguments and returns the triple $(u, v, x \wedge y)$. The code is given in program 79.

It is easy to see that the third component of the returned triple is the `gcd` of $x$ and $y$. Indeed, if we focus on this component, we find the same computations as those performed in the function `gcd` with the variables $x$ and $y$, since $x - \lfloor \frac{x}{y} \rfloor y$ is nothing but $x \bmod y$. To verify the equation 10.1, it suffices to proceed by induction on $y$. The base case $y = 0$ is clear. For the induction step, the induction hypothesis tells us that:

$$uy + v(x - \lfloor \frac{x}{y} \rfloor y) = x \wedge y$$

This is equivalent to:

$$vx + (u - \lfloor \frac{x}{y} \rfloor v)y = x \wedge y,$$

which corresponds to the code of the function `extended_gcd`.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 79 — The extended Euclidean algorithm.

```pascal
let rec extended_gcd x y =
    if y = 0 then
        (1, 0, x)
    else
        let q = x / y in
        let (u, v, g) = extended_gcd y (x - q * y) in
        (v, u - q * v, g)
```

The complexity of `extended_gcd` is the same as that of the function `gcd`. The number of operations performed in each iteration is certainly greater, but remains bounded, and the number of iterations is exactly the same. Therefore, the complexity is still $O(\log x)$.

## 10.2 Exponentiation by Squaring

For $n \in \mathbb{N}$, a naive computation of $x^n$ performs $n-1$ multiplications. The algorithm of exponentiation by squaring consists in computing $x^n$ while performing only $O(\log n)$ multiplications. It makes use of the following identities:

$$x^n = \begin{cases} (x^2)^{n/2} & \text{if } n \text{ is even,} \\ x(x^2)^{(n-1)/2} & \text{if } n \text{ is odd.} \end{cases}$$

Its translation in OCaml, for $x$ of type `int`, is easy. We may write, for example, the code given in program 80. It uses the fact that `n/2` returns $\lfloor n/2 \rfloor$, irrespective of whether $n$ is even or odd.

There are multiple variants. We may, for example, treat $n = 1$ as a special case, but this is not particularly useful. (See also exercise 10.4.) In any case, the central idea remains the following: If we divide the exponent by at least two in each step, the total number of recursive calls is proportionate to $\log(n)$, and so is the number of multiplications.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 80 — Exponentiation by squaring.

```pascal
let rec exp x n =
    if n = 0 then
        1
    else
        let r = exp (x * x) (n / 2) in
        if n mod 2 = 0 then r else x * r
```

More precisely, we can show by induction on $k$ that, if $2^{k-1} \le n < 2^k$, then the function `exp` performs exactly $k$ recursive calls. As each recursive call to `exp` performs one or two multiplications, we deduce that the total number of multiplications is bound by $2 \log n$.

The applications of this algorithm are innumerable, especially since $x$ does not necessarily have to be of type `int`. Once we have a unit and an associative operation, that is, a monoid $M$, we may apply this algorithm to compute $x^n$ for $x \in M$ and $n \in \mathbb{N}$. Exercise 10.6 proposes one such generalization.

## 10.3 Modular Arithmetic

In this section, we assume that we are looking to perform arithmetic modulo $m$ for a fixed value of $m$. This is interesting for several reasons. We may, for instance, wish to compute a very large value but be interested only in its final $k$ digits in base 10. In this case, we will have $m = 10^k$. We could also use the Chinese Remainder Theorem to compute a given value modulo several pairwise relative prime integers and represent in this way integers potentially larger than those provided by the machine. See exercise 10.8.

In the rest of this section, we assume $m$ is a constant of type `int` such that $m > 0$. We will write a certain number of functions that manipulate integers modulo $m$, that is, integers $x$ such that:

$$0 \le x < m. \tag{10.2}$$

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Exercise 10.7 proposes to gather together the various functions in a functor parametrized by $m$. We begin with a function `of_int` that takes an arbitrary integer $x$ and returns its residue modulo $m$. We cannot simply do $x \mod m$ because if $x$ is negative, then $x \mod m$ is also negative. In general, $a \mod b$ has the same sign as $a$. We therefore write:

```h
let of_int x =
    let r = x mod m in if r < 0 then r + m else r
```

Addition modulo $m$ is easily written if we assume that the operands are already integers modulo $m$. To guarantee the property (10.2), it suffices to subtract $m$ if the value of $x+y$ is greater than or equal to $m$.

```h
let add x y =
    let r = x + y in if r >= m then r - m else r
```

The integers $x$ and $y$ being non-negative, the sign is not a problem here. Nevertheless, for the addition to be correct, the operation $x+y$ must not overflow, that is, we must have $x+y \le \max\_int$. An easy way of guaranteeing this is to impose a maximum value on $m$, namely, $m \le \max\_int/2 + 1$. We can do this by placing the following assertion at the beginning of the code:

```h
let () = assert (0 < m && m <= max_int/2 + 1)
```

Subtraction is written just as easily, if we take care to return a non-negative value, as we did above for the function `of_int`.

```h
let sub x y =
    let r = x - y in if r < 0 then r + m else r
```

Arithmetic overflow is not possible here, thanks to the assumption that we made earlier for the addition operation. (It is, nevertheless, worth convincing yourself of this.) The code of these first three operations is given in program 81.

## Multiplication

Multiplication is a subtler affair. For operands that are sufficiently small, that is, less than $\sqrt{m}$, it suffices to perform the multiplication and take its residue:

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 81 — Arithmetic modulo *m* (addition and subtraction).

```pascal
let () = assert (0 < m && m <= max_int/2 + 1)

let of_int x = let r = x mod m in if r < 0 then r + m else r

let add x y = let r = x + y in if r >= m then r - m else r

let sub x y = let r = x - y in if r < 0 then r + m else r
```

```pascal
let mul x y =
  (x * y) mod m
```

However, if we do not wish to restrict the range of the operands, and keep only the restriction already imposed, namely, $m \le \max_{\text{int}}/2 + 1$, we cannot use ordinary multiplication as it could trigger an overflow. We can, however, work around this difficulty by performing the multiplication as one would by hand. Suppose *x* is written in base 2 as:

$$ (x_k x_{k-1} \dots x_1 x_0)_2 $$

where $x_k$ is the most significant digit and $x_0$, the least significant digit. We will then write a loop that calculates $(x_k \dots x_i)_2 \times y$ (mod $m$), for $i$ varying from $k+1$ to 0. We write it as a `for` loop, accumulating the result in a reference `r`. Given the assumption on $m$ and the fact that the OCaml type `int` is represented with `Sys.word_size-1` signed bits (see chapter 3), we can begin with $i = \text{Sys.word_size-4}$.

```pascal
let mul x y =
  let r = ref 0 in
  for i = Sys.word_size - 4 downto 0 do
```

In each iteration of the loop, we begin by multiplying `r` by 2, modulo `m`, which we do with the function `add`.

```pascal
r := add !r !r;
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Program 82 — Arithmetic modulo *m* (multiplication and division).

```pascal
let mul x y =
  let r = ref 0 in
  for i = Sys.word_size - 4 downto 0 do
    r := add !r !r;
    if x land (1 lsl i) <> 0 then r := add !r y
  done;
  !r

let div x y =
  let u,_, g = extended_gcd y m in
  if g <> 1 then invalid_arg "div";
  mul x (of_int u)
```

We then have !r ≡ (x_k...x_i)_2 × 2 × y (mod m). Next, we test the bit x_i and, if it equals 1, we add y to r to re-establish the loop invariant.

```pascal
if x land (1 lsl i) <> 0 then r := add !r y
```

Once out of the loop, we have !r ≡ x_kx_{k-1}...x_1x_0 × y (mod m), that is, !r ≡ x × y (mod m), which is the desired result. The code is given in program 82.

Division

The division of *x* by *y* modulo *m* assumes that *y* and *m* are relatively prime, and returns a result *q* such that *q* × *y* = *x* (mod *m*). We use the extended Euclidean algorithm (program 79). Indeed, if *y* and *m* are relative prime integers, then the extended Euclidean algorithm gives us two integers, *u* and *v*, such that:

$$u \times y + v \times m = 1$$

This follows from equation 10.1. Multiplying this equality by *x*, we get:

$$(x \times u) \times y = x \pmod m$$

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
The result of the division is therefore $x \times u$. The code is given in program 82. It uses `of_int u` because `u` may be negative, but could be made more efficient since Euclid's algorithm guarantees that $|u| < m$, so that the use of `of_int` is unnecessarily complicated.

## 10.4 Matrix Calculus

In this section, we implement some basic matrix operations. For the sake of simplicity, we assume that the matrices have integer coefficients. Exercise 10.9 proposes a generalization to matrices with arbitrary coefficients.

As we explained in section 2.5 *Sieve of Eratosthenes*, a matrix is simply an array of arrays, that is:

```h
type matrix = int array array
```

We recall that we must not write `Array.make 3 (Array.make 4 v)` to create a $3 \times 4$ matrix, but rather `Array.make_matrix 3 4 v`, or equivalently, `Array.init 3 (fun _ -> Array.make 4 v)`.

Other than `Array.make_matrix`, the OCaml standard library does not provide operations on matrices. Let us begin therefore by writing a function `init_matrix`, analog of `Array.init` for matrices. The function takes as arguments the dimensions $n$ and $m$ of the matrix, and a function $f$ to initialize each element.

```h
let init_matrix n m f =
    Array.init n (fun i -> Array.init m (fun j -> f i j))
```

In other words, we return the matrix $M$ of size $n \times m$ such that $M_{i,j} = f_{i,j}$. We derive from this a function `id` to construct the identity matrix of size $n \times n$.

```h
let id n =
    init_matrix n n (fun i j -> if i = j then 1 else 0)
```

Next, let us write a function `size` that returns the dimensions of a matrix, that is, the number of rows and columns, as a pair. Assuming that a matrix always contains at least one row, we may write:

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 83 — Basic matrix operations.

```pascal
type matrix = int array array

let init_matrix n m f =
    Array.init n (fun i -> Array.init m (fun j -> f i j))

let id n =
    init_matrix n n (fun i j -> if i = j then 1 else 0)

let size a =
    (Array.length a, Array.length a.(0))

let add a b =
    let (n, m) as s = size a in
    if size b <> s then invalid_arg "add";
    init_matrix n m (fun i j -> a.(i).(j) + b.(i).(j))

let mul a b =
    let n, p = size a in
    let q, m = size b in
    if q <> p then invalid_arg "mul";
    let product i j =
        let s = ref 0 in
        for k = 0 to p - 1 do s := !s + a.(i).(k) * b.(k).(j) done;
        !s
    in
    init_matrix n m product
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
```pascal
let size a =
  (Array.length a, Array.length a.(0))
```

Let us now consider the addition of two matrices $A$ and $B$. We begin by verifying that they have the same size.

```pascal
let add a b =
  let (n, m) as s = size a in
  if size b <> s then invalid_arg "add";
```

It then suffices to use `init_matrix` to construct the matrix whose general term is $A_{i,j} + B_{i,j}$.

```pascal
init_matrix n m (fun i j -> a.(i).(j) + b.(i).(j))
```

The product of a matrix $A$ of size $n \times p$ and a matrix $B$ of size $p \times m$ is a matrix $C$ of size $n \times m$ whose general term is given by:

$$C_{i,j} = \sum_{k=0}^{k<p} A_{i,k}B_{k,j} \quad \text{for } 0 \le i < n \text{ and } 0 \le j < m. \quad (10.3)$$

We begin by determining the dimensions of the two matrices and by verifying that the number of columns of $A$ is indeed equal to the number of rows of $B$.

```pascal
let mul a b =
  let n, p = size a in
  let q, m = size b in
  if q <> p then invalid_arg "mul";
```

We then write a local function `product` that computes the coefficient $C_{i,j}$ following equation (10.3).

```pascal
let product i j =
  let s = ref 0 in
  for k = 0 to p - 1 do s := !s + a.(i).(k) * b.(k).(j) done;
  !s
```

It only remains to create a matrix of size $n \times m$ with the function `init_matrix` by passing it the function `product` as argument.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
in
init_matrix n m product

The complete code is given above in program 83. An array of arrays is not the only possibility for representing a matrix. There are other representations that are more economical in terms of memory when a large number of elements are equal. We then speak of *sparse matrices*. Exercise 10.13 proposes one such representation.

# 10.5 Exercises

## Euclid's Algorithm

**10.1** The complexity result given in section 10.1 assumes $x > y$. Show that in the general case the complexity is $O(\log(\max(x, y)))$.

**10.2** Euclid's algorithm as presented above assumes $x, y \ge 0$. If $x$ or $y$ is negative, the algorithm could return a negative result. (The operation `mod` returns a value of the same sign as its first argument.) Write a second version of the function `gcd` that always returns a non-negative result irrespective of the sign of its arguments. In which case will the result be zero?

**10.3** Let $x, y, m$ be three positive integers such that $y \land m = 1$. The quotient of $x$ by $y$ modulo $m$ is any integer $w$ such that:

$$0 \le w < m \text{ and } x \equiv yw \pmod{m}.$$

Write a function to compute the quotient of $x$ by $y$ modulo $m$.

## Exponentiation by Squaring

**10.4** Write a variant of the function `exp` based on the following identities:

$$x^n = \begin{cases} (x^{n/2})^2 & \text{if } n \text{ is even,} \\ x((x^{(n-1)/2})^2) & \text{if } n \text{ is odd.} \end{cases}$$

Is there a difference in terms of efficiency?

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
**10.5** Write a tail-recursive version of the function **exp**. This is not very interesting if *n* is a machine integer since log(*n*) is then very small compared to the maximum number of nested calls supported by the stack. However, if *n* is very large, for example *n* = 2<sup>20</sup>, tail recursion is necessary.

**10.6** Write the code of program 80 as a functor computing $x^n$ for $x$ of an arbitrary type $t$, provided that this type $t$ is equipped with a constant **one** and an operation **mul**.

**Modular Arithmetic**

**10.7** Write the code of programs 81 and 82 as a functor parametrized by the value of $m$, that is:

module Modulo(M: sig val m: int end): sig ... end = struct ... end

Why is it useful to write such a functor? Is it useful to have this functor export
an *abstract type* of integers modulo *m*?

**10.8** If $m_1$ and $m_2$ are two relatively prime integers, the Chinese Remainder Theorem allows us to represent every integer $x$ between 0 (inclusive) and $m_1m_2$ (exclusive) by a pair $(x_1, x_2)$, where $x \equiv x_1 \pmod{m_1}$ and $x \equiv x_2 \pmod{m_2}$. This representation is useful when $m_1$ and $m_2$ are representable in the machine but $m_1m_2$ is not (being too large). Write the addition, subtraction, and multiplication operations on this representation of integers. Given the pair $(x_1, x_2)$, how can we recover the integer $x$ that it represents?

More generally, we may consider *k* integers $m_1, m_2, \dots, m_k$ pairwise relatively prime and represent arbitrary integers by *k*-tuples. For more details, consult *The Art of Computer Programming* (vol. 2, sec. 4.3.2).

**Matrix Calculus**

**10.9** Write the code of figure 83 as a functor parametrized by the type of matrix elements, constants **zero** and **one** of this type, and operations **add** and **mul** on this type.

**10.10** Write a function **power**: **t** -> **int** -> **t** that raises a matrix to power *n* using the algorithm of exponentiation by squaring.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
**10.11** The numbers of the Fibonacci sequence ($F_n$) satisfy the following identity:

$$ \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n = \begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix}. $$

In other words, we may compute $F_n$ by raising a $2 \times 2$ matrix to power $n$. Derive a program from the preceding exercise that computes $F_n$ in $O(\log n)$ elementary arithmetic operations.

Note, however, that this does not mean that we are able to compute $F_n$ for large values of $n$. In fact, the terms of the Fibonacci sequence grow exponentially. If we make use of arbitrary-precision integers, the cost of the arithmetic operations themselves must be taken into account, and the complexity will not be $O(\log n)$. Otherwise, the calculation will quickly overflow.

**10.12** Combining the preceding exercise and the modular arithmetic of the preceding section, write a program that computes the final seven digits of $F_{10^6}$. The computation should be instantaneous.

**10.13** When a matrix contains many identical elements (for example, many zero elements), it can be interesting to represent it more compactly than as an array of arrays. One solution consists in representing each row of the matrix by a dictionary that associates some of the column indices with their corresponding elements, the others taking a default value. For example, for matrices with integer coefficients whose zero elements are not represented, we may define the following type:

```python
type matrix = { cols: int; rows: M.t array }
```

where M.t is the type of a dictionary associating integers with integers (for example, obtained using the module `Map`), and where the field `cols` stores the number of columns. Write the functions `id`, `size`, `add`, and `mul` for this type.
— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by/4.0/) —
# Dynamic Programming and Memoization

Dynamic programming and memoization are two closely related techniques based on the idea that it is better to avoid computing the same thing twice.

## 11.1 Basic Idea

We illustrate the two techniques in the simple example of the Fibonacci sequence. Recall that this sequence of integers ($F_n$) is defined by:

$$ egin{cases} F_0 = 0 \\ F_1 = 1 \\ F_n = F_{n-2} + F_{n-1} \text{ if } n \ge 2. \end{cases} $$

It is easy enough to write a recursive function that calculates $F_n$ based on this definition:

```h
let rec fib n =
    if n <= 1 then n else fib (n - 2) + fib (n - 1)
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
But this is also rather simplistic. Indeed, it would take more than a minute to compute $F_{50}$. It is not the recursion that is problematic here, but rather the fact that we compute the same values of $F_n$ several times over. For instance, to compute $F_5$, $F_3$ will be computed twice and $F_2$ three times, as shown in the call graph in figure 11.1.

<figure><img src="image_122.png" /><figcaption>Figure 11.1: Call graph of $F_5$.</figcaption></figure>

It is easy to see that, in general, the computation of $F_n$ requires $F_{n+1} - 1$ additions because the recursion equation is exactly the same as that which defines $F_n$. We can observe this empirically. On a computer circa 2011, it takes 2 seconds to compute $F_{42}$, 3 seconds to compute $F_{43}$, 5 seconds to compute $F_{44}$, etc. We recognize here the terms of the Fibonacci sequence. Extrapolating, we can say that it would take 89 seconds to compute $F_{50}$, and this is correct to within a half-second! As we know that $F_n$ grows exponentially, we cannot hope to go much further in the computation of the sequence using this method.

### 11.1.1 Memoization

The fact that we are computing the same thing several times leads us, quite naturally, to the idea that we could store the intermediate results in a table. Such a table would associate integers *n* with the value $F_n$. Accordingly, we proceed as follows: To compute *fib n*, we begin by checking whether there is an entry in the table for *n*, in which case we return the corresponding value. Otherwise, we calculate *fib n* as *fib (n-2) + fib (n-1)*, that is, recursively. We return the result after having added it to the table. This technique, consisting in using

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by/4.0/) —
a table to store intermediate results, is called *memoization*, a term coined by the researcher Donald Michie in 1968.

Let us implement this idea in a function `fib_memo`. We begin by introducing a hash table to store the intermediate results, choosing an arbitrary initial size for it:

```python
let memo = Hashtbl.create 17
```

The function `fib_memo` begins by checking whether the value of `fib n` is already present in the table `memo`. If it is, it is returned:

```python
let rec fib_memo n =
    try
        Hashtbl.find memo n
```

Otherwise, we compute the result exactly as for the function `fib`, that is, using two recursive calls:

```python
with Not_found ->
    let fn = if n <= 1 then n
              else fib_memo (n-2) + fib_memo (n-1) in
```

Then, we store it in the table `memo` before returning it:

```python
Hashtbl.add memo n fn;
fn
```

This concludes the function `fib_memo`. Its efficiency is far greater than that of `fib`. The computation of $F_{50}$, for example, becomes instantaneous (in contrast to the 89 seconds taken by `fib`). We can show that the complexity of `fib_memo` is linear. Intuitively, we understand that computing $F_n$ now involves computing $F_i$ only once for $i \le n$. The complete code of the function `fib_memo` is given in program 84 below. Note that the table `memo` is defined *outside* the function `fib_memo`, since it must be the same for all recursive calls.

## 11.1.2 Dynamic Programming

Using a hash table to memoize the computation of the Fibonacci sequence may seem unnecessarily costly. Indeed, in the end we will store only the values of

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 84 — Computing $F_n$ via memoization.

```pascal
let memo = Hashtbl.create 17
let rec fib_memo n =
  try
    Hashtbl.find memo n
  with Not_found ->
    let fn =
      if n <= 1 then n else fib_memo (n-2) + fib_memo (n-1) in
    Hashtbl.add memo n fn;
  fn
```

$F_i$ for $i \le n$, for which an array of size $n+1$ would suffice. If we wished to rewrite the function `fib_memo` using this idea, we could, for example, begin by assigning the value $-1$ to each cell of the array, to indicate that the value has not yet been computed. Another possibility would be to fill the array in a specific order. In the present case, we see that it suffices to fill the array in increasing order since the computation of $F_i$ requires the values of $F_{i-2}$ and $F_{i-1}$. This technique, consisting in using a table and progressively filling it with the results of intermediate computations, is called *dynamic programming*, often abbreviated as DP.

Let us implement this idea in a function `fibDP`. We begin by treating the special case $n = 0$. In the general case, when $n > 0$, we allocate an array `f` of size $n+1$, which will contain the values of the $F_i$.

```pascal
let fibDP n =
  if n = 0 then 0 else
  let f = Array.make (n+1) 0 in
```

This array may be allocated inside the function `fibDP` since, unlike `fib_memo`, it will not be recursive. Next, we fill the cells of the array `f` in increasing order, treating `f.(1)` as a special case.

```pascal
f.(1) <- 1;
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
**Program 85** — Computing $F_n$ via dynamic programming.

```pascal
let fibDP n =
    if n = 0 then 0 else
    let f = Array.make (n+1) 0 in
    f.(1) <- 1;
    for i = 2 to n do f.(i) <- f.(i-2) + f.(i-1) done;
    f.(n)
```

```pascal
for i = 2 to n do f.(i) <- f.(i-2) + f.(i-1) done;
```

Once the array is filled, we have only to return the value contained in its last cell.

```pascal
f.(n)
```

This concludes the function `fibDP`. As for `fib_memo`, its complexity is linear. This can be seen easily since the code boils down to a simple `for` loop. The complete code of the function `fib DP` is given in program 85.

## Remark

We have chosen the Fibonacci sequence as an example for pedagogical reasons. Of course, it is very simple to compute the terms of this sequence in linear time without using either memoization or dynamic programming (see exercise 11.1).

## 11.2 Mechanical Memoization

The use of an array in `fib DP`, rather than the hash table in `fibMemo`, may leave us thinking that dynamic programming is simpler to implement than memoization. This is indeed the case in that example.

It is nevertheless important to understand that in order to write `fib DP` we made use of two facts: we had to compute the $F_i$ for all $i \le n$; and we could compute them in increasing order. In general, the inputs of the function that

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
we wish to compute may not necessarily be consecutive integers. They may not be integers at all, or the dependencies between the different values may not be as simple.

Memoization, by contrast, is simpler to implement. It suffices to add a few lines of code to consult and fill the hash table, without modifying the structure of the function. We can do this mechanically.

Better yet, we can implement the memoization procedure *once and for all*, as a higher-order function `memoFun` that takes as argument a function `f` to be memoized. We begin by creating a hash table in which to store the already computed results.

```h
let memoFun f =
    let h = Hashtbl.create 17 in
```

We then return a function that performs the memoization, that is, consults the table and fills it as needed.

```h
fun x ->
    try Hashtbl.find h x
    with Not_found -> let v = f x in Hashtbl.add h x v; v
```

The type of the function `memoFun` is as follows:

```h
memoFun : ('a -> 'b) -> 'a -> 'b
```

As we see, `memoFun` takes as argument a function of an arbitrary type and returns a function of the same type. The polymorphic nature of the function follows from the use of the generic hash tables of OCaml (see exercise 11.2).

In order to memoize the function `fib`, it suffices to write:

```h
let fib_memo2 = memoFun fib
```

If we compute `fib_memo2` 50 twice, the second computation is instantaneous because the result is present in the table. However, the first computation takes just as long as it does without memoization. Indeed, during the first call, we

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
calculate `fib 50` because the value is not present in the table. During this computation, recursive calls are made to the function `fib` and not to its memoized version `fib_memo2`.

We might think that a solution consists in writing the following code:

```pascal
let rec fib_memo3 =
    memoFun (fun n -> if n <= 1 then n
               else fib_memo3 (n-2) + fib_memo3 (n-1))
```

Unfortunately, this still does not work. Each call to `fib_memo3` involves a new call to `memoFun`, which has the effect of allocating a new hash table. Thus the previously calculated results are lost and cannot be found.

To solve this problem, that is, to allocate a single hash table and perform recursive calls to the function being memoized, it is necessary to rewrite the function `memoFun` in a subtler manner. We do this in the form of a function `memo_rec` that may be used as follows:

```pascal
let fib_memo4 =
    memo_rec (fun f n -> if n <= 1 then n else f (n-2) + f (n-1))
```

Here, the single call to `memo_rec` ensures the creation of a single table. The function being memoized is passed as an additional argument `f`.

The function `memo_rec` is written as follows. We begin, as in case of `memoFun`, by creating a table:

```pascal
let memo_rec ff =
    let h = Hashtbl.create 5003 in
```

Next, we define the recursive function `f` that performs the memoization:

```pascal
let rec f x =
    try Hashtbl.find h x
    with Not_found -> let v = ff f x in Hashtbl.add h x v; v
    in
f
```

When the result is not in the table, we compute the result `v` using `ff f x`, that is, by passing the function `f` and the argument `x` to the function `ff` provided by the user. The complete code of `memo_rec` is given in program 86.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 86 — Generic memoization operator.

```pascal
let memo_rec ff =
    let h = Hashtbl.create 5003 in
    let rec f x =
        try Hashtbl.find h x
        with Not_found -> let v = ff f x in Hashtbl.add h x v; v
    in
f
```

As in case of `memoFun`, the type of `memo_rec` is polymorphic, but it is more complex:

`memo_rec: (('a -> 'b) -> 'a -> 'b) -> 'a -> 'b`

The function `ff` has type (`'a -> 'b`), (`'a -> 'b`), which may be read as (`'a -> 'b`), (`'a -> 'b`). That is, it takes as argument a function of type `'a -> 'b` and returns a function of the same type.

## 11.3 Differences between Memoization and Dynamic Programming

In certain situations, dynamic programming may be preferred to memoization. Take the example of the computation of the binomial coefficients $C(n, k)$, the recursive definition of which is as follows:

$$ 
\begin{cases}
C(n, 0) = 1 \\
C(n, n) = 1 \\
C(n, k) = C(n-1, k-1) + C(n-1, k) \text{ if } 0 < k < n.
\end{cases}
$$

Both memoization and dynamic programming may be applied to this definition. In the first case, we will have a hash table indexed by the pair $(n, k)$ and, in the second case, a matrix indexed by $n$ and $k$. However, if we seek to compute

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
$C(n, k)$ for $n = 2 \times 10^5$ and $k = 10^5$, we are probably going to exceed the available memory of the machine (assuming an ordinary desktop computer). In the first case, this would happen when filling the hash table and in the second, when trying to allocate the matrix. The reason is that the complexity in time and in space is $O(nk)$. In the previous example, we would have to store at least 15 billion intermediate results.

Yet, if we look closer, the computation of $C(n, k)$ for a given value of $n$ only requires the values of $C(n-1, k)$. Accordingly, we may carry out the computation for increasing values of $n$, even if not all the values computed are going to be useful in the end. We can visualize the situation by drawing Pascal's triangle as in figure 11.2. We will compute the values line by line, keeping track of only a single line at any given moment.

```
1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 10 10 5 1
...
```

<caption>Figure 11.2: Pascal's triangle.</caption>

Let us implement this idea in a function `comb smart dp`. We begin by allocating an array `row` of size `n+1` that will contain the line of Pascal's triangle being computed:

```pascal
let comb smart dp n k =
    let row = Array.make (n+1) 0 in
```

In fact, an array of size $k+1$ would suffice (see exercise 11.3).

Initially, this array is filled with zeroes. We initialize the first cell with 1.

```pascal
row.(0) <- 1;
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
This value will not change in what follows, because the first column of Pascal's triangle contains only ones. Next, we write a loop to compute the i-th line of Pascal's triangle:

```pascal
for i = 1 to n do
```

This computation will be done in-place within the array `row`, assuming that it contains the values of the `(i - 1)`-th line. To avoid overwriting values that have yet to be used, we will compute the new values from right to left, since these values do not depend on those situated further to the right. We continue with a second loop:

```pascal
for j = i downto 1 do row.(j) <- row.(j) + row.(j-1) done
```

Here, we make use of the fact that `row.(i)` contains 0, which permits us to avoid introducing a special case. We stop when `j = 1` because the value `row.(0)` does not need to be recomputed, as explained above. Once we have exited the double loop, the array `row` contains the n-th line of Pascal's triangle, and we only have to return `row.(k)`:

```pascal
done;
row.(k)
```

The time complexity remains $O(nk)$ but the space complexity is only $O(n)$. In the example given above, with $n = 2 \times 10^5$ and $k = 10^5$, the result is now obtained in a few seconds. Note, however, that it causes an arithmetic overflow (see exercise 11.4).

The conclusion of this little exercise is that dynamic programming, which allows for a finer control of the use of memory, can sometimes be more advantageous than memoization. Conversely, in many situations where this extra measure of control is not necessary, memoization is far simpler to implement.

## 11.4 Hash-consing

Just as we can avoid computing the same thing twice, as we did above, we may similarly avoid constructing the same data several times, to save on memory. We can do this safely as long as we are dealing with values that will not be

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
subsequently modified. If, for example, we have already constructed the list
[1; 2; 3], it is pointless to reconstruct anew. We can instead reuse it.

The technique of reusing values that have already been constructed is known
as *hash-consing*. It is quite simply the application of the technique of memo-
ization to the construction of values. The term hash-consing combines precisely
the idea of hashing used in memoization with that of allocation<fnref n="1" />. Here, we
will demonstrate how to implement this technique, taking as an example binary
trees containing characters, that is:

type tree = E | N of tree * char * tree

The idea is to provide an alternative to the constructor `N`, in the form of
a function `node` that behaves exactly as this constructor, but that returns a
structurally identical value, if one has already been constructed. One possible
solution consists in directly memoizing the function `node`, that is, by storing the
values of type `tree` that have already been constructed in a hash table. However,
proceeding in this manner would not be very efficient. Indeed, both computing
the hash keys and comparing values of the type `tree` require traversals that
may be costly.

A much better solution consists in avoiding these traversals by associating
an integer, to be used in the hash function, with each value of type `tree`. We
begin, therefore, by modifying the type `tree` slightly in order to store an integer
as the first argument of the constructor `N`.

type tree = E | N of int * tree * char * tree

We will see below how this integer is computed. We write a function that
returns this integer when it is given the constructor `N`, and the integer 0 other-
wise.

let unique = function
| E -> 0
| N (u,_,_,_) -> u

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —

<footnote n="1">¹The term *cons* stands for allocation in the language Lisp, where the technique of hash-consing was first introduced.</footnote>
We will justify the name of this function below, by explaining why this
integer is in fact unique. For the present, this is not important since we wish
only to define a hash function. Here is the definition of the hash function:

```h
let hash = function
  | E ->
  | O
  | N (, l, c, r) ->
    (19 * (19 * unique l + Char.code c) + unique r) land max_int
```

This function does not make use of the integer stored in the constructor `N`. Indeed, `hash` will be used precisely to determine if a tree with root `c` and subtrees `l` and `r` has already been constructed. It must therefore only depend on the values `l`, `c`, and `r`. To compute the result in constant time, the function `hash` uses the integers stored in the subtrees `l` and `r`, via the function `unique`. The operation `land max_int` is used here to guarantee a non-negative result.

The second idea consists in comparing two values of type `tree` in constant
time by using the fact that their subtrees have already been shared if they are
equal and are therefore physically equal. For this we write a function `equal`
that uses the physical equality == on the subtrees.

```h
let equal t1 t2 = match t1, t2 with
  | E, E -> true
  | N (, l1, c1, r1), N (, l2, c2, r2) ->
    l1 == l2 && c1 == c2 && r1 == r2
  | _ -> false
```

Equipped with these two functions, we can now construct a hash table to
memoize the construction of trees. We can do this using the functor `Hashtbl.Make`
of the OCaml standard library. However, such a table would keep the con-
structed values forever. What we want, by constrast, is for the GC to be able
to recover the space allocated for values that are no longer used, as it usually
does.

For this, we will use another hash table structure included in the OCaml
standard library, namely, the functor `Weak.Make`. Its signature is identical to
that of the functor `Hashtbl.Make`.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
```python
module X = struct
    type t = tree
    let hash = ...
    let equal = ...
end
module W = Weak.Make(X)
```

The resulting module `W` provides a set structure of values of type `tree`, in which the elements disappear when they are not referenced anywhere except by a pointer from this structure. Such pointers, which are not taken into account by the GC, are called *weak references*. Hence the name of the module: `Weak`.

We may now construct the table that will contain the constructed values of type `tree`. We initialize it with an arbitrary size.

```python
let nodes = W.create 5003
```

It remains only to write the function `node` that will perform the actual memoization. It uses a global counter, initialized with the value 1.

```python
let node =
    let cpt = ref 1 in
    fun l c r ->
```

We begin by constructing the desired tree, using the current value of the counter.

```python
let n0 = N (!cpt, l, c, r) in
```

Next, we consult the hash table to determine if such a value has already been constructed. The module `W` conveniently provides a function `merge` that determines if an equal value is already present in the set and, if this is not the case, adds it. In both cases, the returned value is contained in the table.

```python
let n = W.merge nodes n0 in
```

It only remains to increment the value of the counter if a new node has been added to the table, which is easily determined by comparing `n` and `n0` using physical comparison. In all cases, we return the tree `n`.

```python
if n == n0 then incr cpt;
n
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Program 87 — Hash-consing (code).

```c
    type tree = E | N of int * tree * char * tree
```

```perl
let empty =
    E
```

```perl
let unique = function
    | E -> 0
    | N (u,_,_,_) -> u
```

```perl
module X = struct
    type t = tree
    let hash = function
        | E ->
        |   0
        | N (_, 1, c, r) ->
        |   (19 * (19 * unique 1 + Char.code c) + unique r)
        |   land max_int
```

```perl
    let equal t1 t2 = match t1, t2 with
        | E, E -> true
        | N (_, l1, c1, r1), N (_, l2, c2, r2) ->
        |   l1 == l2 && c1 == c2 && r1 == r2
        | _ -> false
```

```perl
end
```

```perl
module W = Weak.Make(X)
let nodes = W.create 5003
```

```perl
let node =
    let cpt = ref 1 in
    fun l c r ->
        let n0 = N (!cpt, l, c, r) in
        let n = W.merge nodes n0 in
        if n == n0 then incr cpt;
        n
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Program 88 — Hash-consing (interface).

```h
    type tree = private E | N of int * tree * char * tree
    val empty: tree
    val node: tree -> char -> tree -> tree
```

The complete code is given in program 87.

Of course, all of this is correct only if two trees that have the same integer are indeed equal. To guarantee it, we must encapsulate the type `tree` behind an interface that prevents users from arbitrarily setting this integer. One possibility is to make `tree` a private type. Such an interface is given in program 88 (see following page). In this way, we can justify a *posteriori* that the integer stored in each tree is unique in the sense that two trees bearing the same integer are equal. However, two equal values can bear a different integer. Indeed, a value can be constructed, then reclaimed by the GC, then reconstructed identically later, with a different integer.

The advantages of hash-consing are numerous. Firstly, it permits us to economize on memory. Certainly, there is an extra cost linked to memoization, but this cost is of constant time. Indeed, we have been careful to write an equality and a hash function that are of constant cost. The second advantage of hash-consing is that we can now use physical equality `==` instead of structural equality. Indeed, two values of type `tree` are equal if and only if they are physically equal. This is guaranteed by the encapsulation of the type `tree` behind a private type. Thirdly, we have a total order on the type `tree`, of constant time, by comparing the values returned by the function `unique`. In particular, we can now construct data structures, for example by using the functors `Set.Make` or `Map.Make`, that will be as efficient as in case of integers. Finally, we have an excellent hash function on the type `tree`, namely, the function `unique`. Exercise 11.6 proposes using it to memoize a recursive function on trees.

It is important to reiterate here that the hash-consing technique is only applicable to structures that can be shared. The technique is thus applicable to persistent structures, but not to mutable ones.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
## 11.5 Exercises

### Dynamic Programming and Memoization

11.1 Modify the function `fib DP` so that it uses only two integers instead of an array.

11.2 Reimplement the function `memo_rec` as a functor parametrized by the type of the arguments of the memoized function, so as to be able to use hash tables other than the generic tables of OCaml.

11.3 Modify the function `comb smart DP` so that it does not compute the values of Pascal's triangle beyond the k-th column.

11.4 Modify the function `comb smart DP` so that it returns an arbitrary-precision integer (type `Num.num` of the OCaml library). Note: The complexity is no longer $O(nk)$ because additions are no longer atomic operations. Their cost depends on the size of the operands, and these grow rapidly in Pascal's triangle.

11.5 The edit distance between two strings is defined as the smallest number of insertions, deletions, and substitutions of characters required to pass from one to the other. For example, the distance between "duck" and "duke" is two (one deletion, one insertion). Write a recursive function that takes as arguments two strings `a` and `b`, and two indices `i` and `j`, and returns the distance between the strings `a[0..i]` and `b[0..j]`. Explain why we end up making the same calls several times. Improve efficiency by using either memoization or dynamic programming.

### Hash-consing

11.6 Write a memoized function `height` that computes the height of a tree making use of the fact that trees have been hash-consed.

11.7 Write a structure of lists of integers making use of the hash-consing technique.

11.8 Generalize the above exercise to a functor that takes as argument the type of the elements, equipped with an equality and a hash function.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
**11.9** How can the application of the constructor N in the function node in program 87 be avoided when the value has already been constructed?

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by/4.0/) —
— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by/4.0/) —
Sorting Algorithms

This chapter presents several sorting algorithms on lists and arrays. We assume that the elements to be sorted are of type `elt`, equipped with a total order given in the form of a function `compare` of type `elt -> elt -> int`. Sorting functions on lists and arrays will respectively have type `elt list -> elt list` and `elt array -> unit`. We assume that we sort lists and arrays in increasing order. For the sake of convenience, we define two functions `lt` and `le` as follows:

let lt x y = compare x y < 0

let le x y = compare x y <= 0

Throughout this chapter, $N$ will denote the number of elements to be sorted. For each sorting algorithm presented, we indicate its complexity not only in terms of the number of comparisons performed, but also in terms of the number of basic allocations in the case of lists and the number of assignments in the case of arrays. For each one, we give the complexity in the best, worst, and average cases. For the average case, we assume that the elements to be sorted are distinct and that the $N!$ possible permutations of the input are all equally likely. Recall that the optimal complexity of a sorting algorithm only performing comparisons of elements is $O(N \log N)$ (see section 12.5).

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
## 12.1 Insertion Sort

Insertion sort consists in successively inserting each element into the set of elements already sorted. This is what we do, for example, when we sort a pack of cards.

### 12.1.1 Lists

Insertion sort on a list consists in traversing it while successively inserting each element in the part already sorted. We begin therefore by writing a function `insert` that inserts an element `x` in a list that is assumed to be sorted. The code is as follows:

```h
let rec insert x = function
| y :: l when lt y x -> y :: insert x l
| l -> x :: l
```

The first pattern matching case corresponds to a recursive insertion when the first element of the list is smaller than `x`. The second pattern matching case corresponds to that in which `x` is smaller than all the elements of `l`, whether `l` is empty or not.

Sorting then consists in traversing the list using `List.fold_left`, with the portion of the list already sorted as accumulator.

```h
let insertion_sort l =
  List.fold_left (fun acc x -> insert x acc) [] l
```

This code may prove problematic for long lists: Since the function `insert` is not tail recursive, it can overflow the stack. It is, however, easy to modify the function `insert` to make it tail recursive. We can do this by adding an accumulator containing the list, in reverse order, of elements smaller than `x` that have already been considered. If we call this accumulator `acc`, we then have the following situation when we try to insert `x` in `l`:

<figure><img src="image_123.png" /></figure>

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Once the position of x is determined, it suffices to reverse the accumulator acc and concatenate the result with 1. The function List.rev_append does exactly this and is tail recursive. The final code of insertion sort is given in program 89.

**Program 89 — Insertion sort on lists.**

```hkl
let rec insert acc x = function
  | y :: l when lt y x -> insert (y :: acc) x l
  | l -> List.rev_append acc (x :: l)

let insertion_sort l =
  List.fold_left (fun r x -> insert [] x r) [] l
```

## Complexity

When the function insert inserts an element x at position k in l, it performs k comparisons and 2k allocations (k to construct the accumulator and k in List.rev_append). At best, k is always 1 and, at worst, it equals the length of the list, which in the end gives the table in figure 12.1.

<table><caption>Figure 12.1: Complexity of insertion sort on lists.</caption><thead><tr><th></th><th>best case</th><th>average</th><th>worst case</th></tr></thead><tbody><tr><td>comparisons</td><td>N</td><td>N<sup>2</sup>/4</td><td>N<sup>2</sup>/2</td></tr><tr><td>allocations</td><td>N</td><td>N<sup>2</sup>/2</td><td>N<sup>2</sup></td></tr></tbody></table>

### 12.1.2 Arrays

Insertion sort on an array a is performed in-place. It consists in successively inserting each element a.(i) in the portion of the array a[0..i-1] that is already sorted, which corresponds to the following situation:

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
<table><thead><tr><th>0</th><th>i-1</th></tr></thead><tbody><tr><td>...</td><td>already sorted ...</td></tr><tr><td></td><td>a.(i)</td></tr><tr><td></td><td>...</td><td>to be sorted ...</td></tr></tbody></table>

We begin with a `for` loop to traverse the array:

```pascal
let insertion_sort a =
    for i = 1 to Array.length a - 1 do
        let v = a.(i) in
```

To insert the element `a.(i)` in the right place, we use a `while` loop that shifts the elements to the right as long as they are greater than `a.(i)`.

```pascal
let j = ref i in
while 0 < !j && lt v a.(!j - 1) do
    a.(!j) <- a.(!j - 1);
decr j
done;
```

Once out of the loop, all that remains is to put `a.(i)` in its place.

```pascal
a.(!j) <- v
done
```

The complete code of insertion sort on arrays is given in program 90.

## Complexity

Note that the function `insertion_sort` performs as many comparisons as assignments. When the `while` loop inserts the element `a.(i)` at position `i - k`, it performs `k + 1` comparisons. At best, `k` equals 0 and, at worst, it equals `i`, which in the end gives the table in figure 12.2.

## 12.2 Quicksort

Quicksort relies on the method of *divide and conquer*: The elements to be sorted are divided into two subsets, all the elements of the first being smaller than those

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 90 — Insertion sort on arrays.

```pascal
let insertion_sort a =
  for i = 1 to Array.length a - 1 do
    let v = a.(i) in
    let j = ref i in
    while 0 < !j && lt v a.(!j - 1) do
      a.(!j) <- a.(!j - 1);
      decr j
    done;
    a.(!j) <- v
  done
```

<table><caption>Figure 12.2: Complexity of insertion sort on arrays.</caption><thead><tr><th></th><th>best case</th><th>average</th><th>worst case</th></tr></thead><tbody><tr><td>comparisons</td><td>N</td><td>N<sup>2</sup>/4</td><td>N<sup>2</sup>/2</td></tr><tr><td>assignments</td><td>N</td><td>N<sup>2</sup>/4</td><td>N<sup>2</sup>/2</td></tr></tbody></table>

of the second. Then, we sort each subset recursively. In practice, the subdivision is performed with the help of an arbitrary element *p* of the set, called the *pivot*. The two subsets are then respectively the elements less than or equal to *p*, and the elements greater than *p*.

### 12.2.1 Lists

To implement quicksort on lists, we begin by writing a function *partition* that partitions a list at a given pivot *p*. To ensure tail recursion, we use an accumulator consisting of a pair of lists. The function *partition* therefore takes the following form:

```pascal
let rec partition ((left, right) as acc) p = function
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
The lists `left` and `right` represent the elements that are, respectively, less than or equal to `p` and greater than `p`. When we reach the empty list, it suffices to return the accumulator.

| [] -> acc

Otherwise, we consider the first element `x` and compare it with the pivot `p`. If `x` is less than or equal to `p`, we add it to `left` and call `partition` recursively on the rest of the list, `s`.

| x :: s when le x p ->
  partition (x :: left, right) p s

On the other hand, if `x` is greater than `p`, we add it to `right` and similarly call `partition` recursively.

| x :: s ->
  partition (left, x :: right) p s

Quicksort takes the form of a recursive function `quicksort`. The base case is that of the empty list.

```hkl
let rec quicksort = function
| [] ->
  []
```

If the list is not empty, we choose its first element as pivot `p` and partition the remaining elements using the function `partition`.

| p :: s ->
  let (left, right) = partition ([],[]) p s in

Finally, we sort `left` and `right` recursively and concatenate the results, remembering to insert the pivot `p` in the middle.

(quicksort left) @ (p :: quicksort right)

The complete code of quicksort on lists is given in program 91.

It is important to note that neither of the two recursive calls of the function `quicksort` is tail recursive. A stack overflow is therefore possible. Exercise 12.4 discusses a solution to this problem.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 91 — Quicksort on lists.

```pascal
let rec partition ((left, right) as acc) p = function
| [] -> acc
| x :: s when le x p ->
    partition (x :: left, right) p s
| x :: s ->
    partition (left, x :: right) p s

let rec quicksort = function
| [] ->
    []
| p :: s ->
    let (left, right) = partition ([],[]) p s in
    (quicksort left) @ (p :: quicksort right)
```

Complexity

The function `partition` makes as many comparisons and allocations as there are elements in the list passed as argument. Let $C(N)$ be the number of comparisons performed in the function `quicksort` for a list of length $N$.

If `partition` returns a list of length $K$ and a list of length $N - 1 - K$, we deduce:

$$C(N) = N - 1 + C(K) + C(N - 1 - K).$$

The worst case corresponds to $K = 0$, which gives $C(N) = N - 1 + C(N - 1)$, and hence $C(N) \sim \frac{N^2}{2}$. The best case corresponds to a list split into two equal halves, that is, $K = N/2$. From this we easily deduce $C(N) \sim N \log N$. For the number of comparisons on average, we consider that the $N$ possible final places

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
for the pivot are all equally probable, which gives:

$$ 
\begin{aligned}
C(N) &= N - 1 + \frac{1}{N} \sum_{0 \le K \le N-1} C(K) + C(N - 1 - K) \\
&= N - 1 + \frac{2}{N} \sum_{0 \le K \le N-1} C(K).
\end{aligned}
$$

After a bit of algebra, left to the reader, we obtain:

$$ \frac{C(N)}{N+1} = \frac{C(N-1)}{N} + \frac{2}{N+1} - \frac{2}{N(N+1)}. $$

This is a telescoping sum, which leads us to conclude that:

$$ \frac{C(N)}{N+1} \sim 2 \log N $$

We therefore have $C(N) \sim 2N \log N$. For the number of allocations $A(N)$, we proceed in a similar fashion, with a slightly different equation:

$$ A(N) = N + K + A(K) + A(N - 1 - K). $$

In the end, we obtain the results given in the table of figure 12.3.

<table><caption>Figure 12.3: Complexity of quicksort on lists.</caption><thead><tr><th></th><th>best case</th><th>average</th><th>worst case</th></tr></thead><tbody><tr><td>comparisons</td><td>$N \log N$</td><td>$2N \log N$</td><td>$N^2/2$</td></tr><tr><td>allocations</td><td>$\frac{3}{2}N \log N$</td><td>$4N \log N$</td><td>$N^2$</td></tr></tbody></table>

### 12.2.2 Arrays

Quicksort on arrays is performed in-place. The partition and sorting functions take as arguments the array and two indices that delimit the portion to be considered. We retain the idea of a function `partition` that organizes the elements around a pivot, and a sorting function that proceeds recursively.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
The function `partition` takes an array `a` as argument, and two indices `l` and `r`, that delimit a portion of the array. This portion corresponds to the elements between the indices `l` (inclusive) and `r` (exclusive). We begin by choosing `a.(1)` as pivot.

```hklm
let partition a l r =
    let p = a.(1) in
```

This choice is arbitrary and can influence performance. Exercise 12.8 proposes a better method of choosing the pivot.

The idea consists in traversing the array from left to right, between indices `l` (inclusive) and `r` (exclusive), with a `for` loop. In each iteration of the loop, the situation is as follows:

<table><tr><td>l</td><td>m</td><td>i</td><td>r</td></tr><tr><td>p</td><td>&le; p</td><td>&gt; p</td><td>?</td></tr></table>

The index `i` of the loop gives the position of the next element to consider. The index `m` partitions the part that has already been traversed. More precisely, it is the index of the last cell containing a value less than or equal to `p`. We represent it with a reference.

```hklm
let m = ref l in
for i = l + 1 to r - 1 do
```

If `a.(i)` is strictly greater than `p`, we leave it in its place. Otherwise, to maintain the loop invariant, it suffices to increment `m` and swap `a.(i)` and `a.(!m)`.

```pascal
if le a.(i) p then begin incr m; swap a i !m end
done;
```

The code of the function `swap` is given in program 92 below. Once all the elements have been traversed, we perform one more swap to put the pivot in its place, that is, at position `!m`, and we return this index.

```pascal
if l <> !m then swap a l !m;
!m
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
We write the main function of quicksort as a recursive function `quick_rec` that takes the same arguments as the function `partition`. If $1 \ge r-1$, the number of elements to be sorted is either zero or one, and there is nothing to do.

```pascal
let rec quick_rec a l r =
    if l < r - 1 then begin
```

Otherwise, we partition the elements between `l` and `r`.

```pascal
let m = partition a l r in
```

After this call, the pivot `a.(m)` is at its final place. We then perform two recursive calls to sort, respectively, the elements smaller than and larger than the pivot.

```pascal
    quick_rec a l m;
    quick_rec a (m + 1) r
end
```

To sort an array, it suffices to call `quick_rec` on all of its elements.

```pascal
let quicksort a = quick_rec a 0 (Array.length a)
```

The complete code is given in program 92.

As with quicksort on lists, a stack overflow could occur. Indeed, the first call to `quick_rec` is not tail recursive and, in the worst case, the number of these calls may be equal to the length of the array. Exercise 12.6 discusses a solution to this problem.

## Complexity

The function `partition` always performs exactly $(r-1) - 1$ comparisons. The complexity in number of comparisons is therefore the same as in case of quicksort on lists. As far as the number of assignments is concerned, we note that the function `partition` performs as many calls to `swap` as increments of `m`. The best case occurs when, at each partition stage, the pivot remains in the first position. Then no assignment is performed. It is important to note that this case does not correspond to the best complexity in terms of comparisons, which is quadratic in this case. Thus, on average, and assuming that all final positions for the

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Program 92 — Quicksort on arrays.

```pascal
let swap a i j =
    let t = a.(i) in
    a.(i) <- a.(j);
    a.(j) <- t

let partition a l r =
    let p = a.(1) in
    let m = ref l in
    for i = l + 1 to r - 1 do
        if le a.(i) p then begin incr m; swap a i !m end
        done;
        if l <> !m then swap a l !m;
        !m

let rec quick_rec a l r =
    if l < r - 1 then begin
        let m = partition a l r in
        quick_rec a l m;
        quick_rec a (m + 1) r
    end

let quicksort a = quick_rec a 0 (Array.length a)
```
pivot are equally likely, we have fewer than $(r-1)+1$ assignments, each call to `swap` performing two assignments. This leads to a calculation analogous to that of the average number of comparisons. We have already done this computation for lists and obtained $2N \log N$. In the worst case, the pivot is always found at position $r-1$. The function `partition` then performs $2(r-1)$ assignments. There is, thus, a total of $N^2$ assignments. In the end, we obtain the results given in the table of figure 12.4.

<table><caption>Figure 12.4: Complexity of quicksort on arrays.</caption><thead><tr><th></th><th>best case</th><th>average</th><th>worst case</th></tr></thead><tbody><tr><td>comparisons</td><td>$N \log N$</td><td>-</td><td>$2N \log N$</td><td>$N^2/2$</td></tr><tr><td>assignments</td><td>-</td><td>0</td><td>$2N \log N$</td><td>$N^2$</td></tr></tbody></table>

To avoid the worst case, to the extent possible, the two optimizations proposed in exercises 12.8 and 12.9 must be implemented.

## 12.3 Merge Sort

As with quicksort, merge sort also applies the principle of divide and conquer. It partitions the elements to be sorted into two parts of the same size, without trying to compare their elements. It then sorts the two parts recursively and merges them; hence the name *merge sort*. Thus, the worst case of quicksort, in which the two parts have disproportionate sizes, is avoided.

### 12.3.1 Lists

We begin by writing a function `split` that partitions a list into two lists of length differing by at most one. To ensure tail recursion, we pass two lists `11` and `12` as accumulators. When we reach the end of the list, we return the pair `(11, 12)`.

```python
let rec split 11 12 = function
    | [] -> (11, 12)
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Otherwise, we consider the first element `x` and place it alternatively in `11` or `12`. The simplest way to achieve this is to exchange `11` and `12` at each call.

```hml
| x :: 1 -> split (x :: 12) 11 1
```

The second step consists in writing a function `merge` that merges two sorted lists. For this, we traverse the two lists simultaneously.

```hml
let rec merge 11 12 = match 11, 12 with
```

If one of the two is empty, we return the other.

```hml
| [], 1 | 1, [] ->
  1
```

Otherwise, we compare the first elements `x1` and `x2` of each list. If `x1` is smaller, it is placed at the head of the result, and we recursively merge the remainder `s1` of `11` with `12`. We proceed symmetrically if `x2` is smaller.

```hml
| x1 :: s1, x2 :: s2 ->
  if le x1 x2 then x1 :: merge s1 12 else x2 :: merge 11 s2
```

Note: We observe that `merge` is not tail recursive. We can fix this easily by adding an additional argument in which we accumulate the smaller elements.

Finally, we write a function `mergesort` that implements merge sort on a list `1`. If the list contains at most one element, we return it directly.

```hml
let rec mergesort l = match l with
```

```hml
| [] | [_] ->
  1
```

Otherwise, we partition the list `l` into two lists using the function `split`:

```hml
| _ ->
  let l1, l2 = split [] [] l in
```

We then sort the two lists `11` and `12` recursively, and merge the results using `merge`.

```hml
merge (mergesort 11) (mergesort 12)
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 93 — Merge sort on lists.

```ocaml
let rec split l1 l2 = function
| [] -> (l1, l2)
| x :: l -> split (x :: l2) l1 l

let rec merge acc l1 l2 = match l1, l2 with
| [], l | l, [] ->
    List.rev_append acc l
| x1 :: s1, x2 :: s2 ->
    if le x1 x2 then merge (x1 :: acc) s1 l2
    else merge (x2 :: acc) l1 s2

let rec mergesort l = match l with
| [] | [] ->
    l
| _ ->
    let l1, l2 = split [] [] l in
    merge [] (mergesort l1) (mergesort l2)
```

Unlike the function `merge`, the function `mergesort` does not risk triggering a stack overflow: Since the list `l` is partitioned into two equal halves, the size of the stack is logarithmic. Program 93 below contains the complete code, with the function `merge` in its tail-recursive version.

We can avoid allocating the lists `l1` and `l2` by cleverly reusing the list `l` itself. Exercise 12.10 proposes such an optimization. This idea is notably used in the function `List.sort` of the OCaml standard library.

Complexity

The function `split` makes $N$ allocations and no comparisons. The function `merge` makes twice as many allocations as comparisons: In the second branch

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
of the pattern matching, there are as many comparisons as allocations, and in the first branch, the number of allocations performed by `List.rev_append` is equal to the length of `acc`, itself equal to the number of allocations/comparisons performed in the recursive part. Let $C(N)$ (resp. $f(N)$) be the total number of comparisons performed by `mergesort` (resp. `merge`). We then have the following equation:

$$C(N) = 2C(N/2) + f(N)$$

Indeed, the two recursive calls are made on lists of the same length $N/2$. In the best case, the function `merge` examines only the elements of one of the two lists, since they are all smaller than those of the other list. In this case, $f(N) = N/2$ and hence $C(N) \sim \frac{1}{2}N \log N$. In the worst case, `merge` examines every element, so that $f(N) = N - 1$; hence $C(N) \sim N \log N$. The analysis on the average case is subtler (see [15, ex 2 p. 646]) and gives $f(N) = N + O(1)$; hence $C(N) \sim N \log N$ as well. The number of assignments is easily deduced from the number of comparisons. In the end, we obtain the results given in the table of figure 12.5.

<table><caption>Figure 12.5: Complexity of merge sort on lists.</caption><thead><tr><th></th><th>best case</th><th>average</th><th>worst case</th></tr></thead><tbody><tr><td>comparisons</td><td><sup>1</sup>&frasl;<sub>2</sub> N log N</td><td>N log N</td><td>N log N</td></tr><tr><td>allocations</td><td>2N log N</td><td>3N log N</td><td>3N log N</td></tr></tbody></table>

## 12.3.2 Arrays

The idea of merge sort on arrays is the same as for lists: We partition the elements to be sorted into two equal halves. We sort the two halves, and then merge them. We delimit the portion to be sorted by two indices `l` (inclusive) and `r` (exclusive). For the partition, it suffices to compute the middle index $m = \frac{l+r}{2}$. We then recursively sort the two parts delimited by `l` and `m` on the one hand, and `m` and `r` on the other. We have only to merge them. This proves extremely difficult to do in-place. It is simpler to use a second array, allocated once and for all at the start of the sort.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
We begin by writing the function `merge` that implements the merging. It takes as arguments two arrays, `a1` and `a2`, and the three indices `l`, `m`, and `r`. The portions `a1[l..m]` and `a1[m..r]` are sorted. The objective is to merge them in `a2[l..r]`. For this, we traverse the two portions of `a1` (with two references `i` and `j`) and the portion of `a2` that is to be filled (with a `for` loop).

```pascal
let merge a1 a2 l m r =
    let i = ref l in
    let j = ref m in
    for k = l to r - 1 do
```

In each iteration of the loop, the situation is therefore as follows:

<table><thead><tr><th></th><th>l</th><th>m</th><th>r</th></tr></thead><tbody><tr><td>a1</td><td>sorted</td><td>sorted</td><td></td></tr><tr><td></td><td>↑<sub>i</sub></td><td>↑<sub>j</sub></td><td></td></tr><tr><td>a2</td><td>sorted</td><td></td><td></td></tr><tr><td></td><td></td><td>↑<sub>k</sub></td><td></td></tr></tbody></table>

We must determine the next value to put in `a2.(k)`. It is the smallest of the two values `a1.(!i)` and `a1.(!j)`. It is important, however, to also treat the case in which there are no more elements in one of the two halves. We determine whether the element is to be taken from the left half with the following test:

```pascal
if !i < m && (!j = r || le a1.(!i) a1.(!j)) then begin
```

This test checks whether the left half is non-empty and its first element is smaller than the first element of the right half, if such an element exists. In both cases, we copy the element into `a2.(k)` and increment the corresponding index.

```pascal
a2.(k) <- a1.(!i); incr i
end else begin
a2.(k) <- a1.(!j); incr j
end
done
```

Next, we implement merge sort as a function `mergesort`. We begin by allocating a temporary array `tmp` obtained by copying the array to be sorted.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
```h
let mergesort a =
    let tmp = Array.copy a in
```

The recursive part of merge sort is materialized via a recursive function `mergesort_rec` that takes as arguments the indices `l` and `r`, that delimit the portion to be sorted.

```h
let rec mergesort_rec l r =
```

If the segment contains at most one element, that is, if $l \ge r-1$, there is nothing more to be done.

```h
if l < r - 1 then begin
```

Otherwise, we divide the interval into two equal halves by computing the middle position `m`, and then recursively sort `a[l..m]` and `a[m..r]`.

```h
let m = (l + r) / 2 in
mergesort_rec l m;
mergesort_rec m r;
```

Note that the computation of $(l+r)/2$ cannot trigger an integer overflow, since the maximum size of an OCaml array is much smaller than the largest representable integer (see section 3.2 *Runtime Model* for further details).

It remains to perform the merge. For this, we copy the entire portion `a[l..r]` into the array `tmp` and then call the function `merge`.

```h
Array.blit a l tmp l (r - l);
merge tmp a l m r
end
```

Finally, we sort the full array `a` by calling `mergesort` on all its elements.

```h
in
mergesort_rec 0 (Array.length a)
```

Program 94 contains the complete code. Its efficiency can be further improved: As with quicksort, we can switch to insertion sort when the portion to be sorted becomes sufficiently small (see exercise 12.7). Another idea consists in avoiding the copy of `a` to `tmp`, which uses `Array.blit` (see exercise 12.12).

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 94 — Merge sort on arrays.

```pascal
let merge a1 a2 l m r =
  let i = ref l in
  let j = ref m in
  for k = l to r - 1 do
    if !i < m && (!j = r || le a1.(!i) a1.(!j)) then begin
      a2.(k) <- a1.(!i); incr i
    end else begin
      a2.(k) <- a1.(!j); incr j
    end
  end
  done
```

```pascal
let mergesort a =
  let tmp = Array.copy a in
  let rec mergesort_rec l r =
    if l < r - 1 then begin
      let m = (l + r) / 2 in
      mergesort_rec l m;
      mergesort_rec m r;
      Array.blit a l tmp l (r - l);
      merge tmp a l m r
    end
  end
  in
mergesort_rec 0 (Array.length a)
```
Complexity

The number of comparisons is exactly the same as in case of lists. The number of assignments is also the same: *N* assignments in the function merge (each element being copied from a1 to a2) and *N* assignments performed by Array.blit.

Let *A(N)* be the total number of assignments in *mergesort*. We then have:

$$A(N) = 2A(N/2) + 2N,$$

that is, $2N \log N$ assignments in total. In the end, we obtain the results of figure 12.6.

<table><caption>Figure 12.6: Complexity of merge sort on arrays.</caption><thead><tr><th></th><th>best case</th><th>average</th><th>worst case</th></tr></thead><tbody><tr><td>comparisons</td><td>1/2 *N* log *N*</td><td>*N* log *N*</td><td>*N* log *N*</td></tr><tr><td>assignments</td><td>2*N* log *N*</td><td>2*N* log *N*</td><td>2*N* log *N*</td></tr></tbody></table>

12.4 Heapsort

Heapsort is based on the use of a priority queue. Such a structure permits the addition of elements and the extraction of its smallest element. The idea behind heapsort is the following: We construct a priority queue containing all the elements to be sorted and then extract the elements one by one. Since the extraction operation returns the smallest element each time, the elements come out in increasing order.

In principle, any priority queue structure could be used. However, a structure in which the addition and extraction operations have cost $O(\log N)$ gives an optimal sorting algorithm in $O(N \log N)$. Both of the heap structures presented in chapter 6 have such complexity.

12.4.1 Lists

We assume we are given a functor `PriorityQueue` that implements a priority queue structure. In principle, we can choose between imperative and persistent

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
structures. However, there is no advantage to choosing a persistent structure here. We therefore use the imperative heap structure described in section 6.3 *Imperative Priority Queues*. Its interface is given again in program 95.

**Program 95 — Minimal signature for imperative priority queues.**

```pascal
type t

type elt

val create : unit -> t

val is_empty : t -> bool

val add : elt -> t -> unit

val get_min : t -> elt

val remove_min : t -> unit
```

We apply the functor `PriorityQueue` by passing as argument a module that compares elements with `compare`. To avoid having to reverse the list at the end, we order elements in decreasing order, that is, the element of highest priority for the priority queue will be the largest element for the order relation `lt`. We therefore write:

```python
module Heap = PriorityQueue(struct
    type t = elt
    let compare x y = compare y x
end)
```

The function `heapsort` takes a list `l` of elements to be sorted and begins by
creating an empty heap `h`.

```python
let heapsort l =
    let h = Heap.create () in
```

The function then adds all the elements of the list `l` to the heap.

```pascal
List.iter (fun x -> Heap.add x h) l;
```

Next, we create a reference `res` that will receive the sorted result. We extract the elements from the heap with the help of a `while` loop.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Program 96 — Heapsort on lists.

```pascal
let heapsort l =
  let h = Heap.create () in
  List.iter (fun x -> Heap.add x h) l;
  let res = ref [] in
  while not (Heap.is_empty h) do
    let x = Heap.get_min h in
    Heap.remove_min h;
    res := x :: !res
  done;
!res
```

```pascal
let res = ref [] in
while not (Heap.is_empty h) do
```

If the heap is not empty, we extract its smallest element x using `Heap.get_min` and then remove it from the heap using `Heap.remove_min`.

```pascal
let x = Heap.get_min h in
Heap.remove_min h;
```

We add x at the head of the list contained in res and then repeat.

```pascal
res := x :: !res
done;
```

Once the heap becomes empty, it suffices to return the contents of res.

```pascal
!res
```

Note that the fact of having reversed the order relation in the heap structure means that we do not have to reverse the list `!res` at the end of `heapsort`. Program 96 contains the complete code.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
Complexity

The complexity of heapsort depends on the complexity of the heap structure used. As we mentioned in the introduction to this section, if the operations `add` and `remove_min` of the heap structure have at most logarithmic cost, then the sort itself has cost $O(N \log N)$ in the worst case, which is optimal. If we now assume that the heap structure used is that of section 6.3 Imperative Priority Queues, we can give a more precise result. Indeed, we then know that `add` and `remove_min` each make a maximum of $2 \log k$ comparisons to add or remove an element from a heap containing $k$ elements, so that the total number of comparisons in the worst case is at most:

$$4(\log 1 + \log 2 + \dots + \log N) \sim 4N \log N$$

In fact, we will see in the next section that the construction of the heap has only linear cost, so that the total is equivalent to $2N \log N$.

12.4.2 Arrays

To implement heapsort on an array `a`, we will construct the heap structure directly inside this array. The organization of the heap in the array is exactly the same as in section 6.3 Imperative Priority Queues: The left and right children of a node stored at index `i` are respectively stored at indices $2i + 1$ and $2i + 2$. As in the previous section, we construct a heap for the reverse order relation, that is, a heap in which the largest element is to be found at the root.

To construct the heap, we consider the elements of the array from right to left. In each iteration, we have the following situation:

<table><thead><tr><th>0</th><th>k</th><th>k+1</th><th>n</th></tr></thead><tbody><tr><td>?</td><td></td><td>heap being built</td><td></td></tr></tbody></table>

The portion `a[k+1..n]` contains the bottom of the heap being constructed, that is, it is a forest of heaps whose roots are situated at indices `i` such that `k < i < 2k + 3`. We then make the value `a.(k)` descend to its place in the heap

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
rooted at k. Once all the elements have been traversed, we have a single heap rooted at 0.

The second step consists in deconstructing the heap. For this, we exchange its root *r* in *a*.(0) with the element *v* in *a*.(n-1). The value *r* is then in its final place. We then reestablish the heap structure on *a*[0..n-1] by making *v* descend to its place in the heap of size n-1 rooted at 0. We repeat the operation for positions n-1, n-2, etc. In each iteration *k*, we have the following situation:

<table><thead><tr><th>0</th><th>k</th><th>n</th></tr></thead><tbody><tr><td>heap</td><td>≤</td><td>sorted</td></tr></tbody></table>

There is a heap in the portion *a*[0..k], all of whose elements are smaller than those of the portion *a*[k..n] that is already sorted.

The two steps of the above algorithm use the same operation consisting in making a value descend until it reaches its position in a heap. This is implemented using a recursive function *move_down* that takes as arguments an array *a*, an index *k*, a value *v*, and an upper bound on the indices, *n*.

let rec move_down a k v n =

We assume that we already have a heap *h*₁, rooted at 2k+1, so that 2k+1 < *n*, and similarly a heap *h*₂, rooted at 2k+2, so that 2k+2 < *n*. The objective is to construct a heap rooted at *k*, containing *v* and all the elements of *h*₁ and *h*₂. We begin by determining whether the heap rooted at *k* will be reduced to a leaf, that is, whether the heap *h*₁ is empty. If this is so, we assign the value *v* to *a*.(k), and we are done.

let r = 2 * k + 1 in
  if r >= n then
    a.(k) <- v

Otherwise, we determine the index *rmax* of the largest of the two roots of *h*₁ and *h*₂, by carefully handling the case in which *h*₂ is empty.

else
  let rmax =
    if r+1 < n then if lt a.(r) a.(r+1) then r+1 else r
    else

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
else r in

If the value `v` is greater than or equal to `a. (rmax)`, the descent is complete, and it suffices to assign `v` to `a.(k)`.

if le a.(rmax) v then a.(k) <- v

Otherwise, we move `a.(rmax)` up and continue the descent of `v` recursively at the position `rmax`.

else begin a.(k) <- a.(rmax); move_down a rmax v n end

The sorting function takes an array `a` as argument:

let heapsort a =

let n = Array.length a in

We begin by constructing the heap bottom-up by calling `move_down`. We avoid unnecessary calls for heaps reduced to single leaves by beginning the loop at $\lfloor \frac{n}{2} \rfloor - 1$. (Indeed, for every strictly greater index `k`, we have $2k+1 \ge n$.)

for k = n/2 - 1 downto 0 do move_down a k a.(k) n done;

Once the heap is fully constructed, we extract its elements one by one, in decreasing order. As explained above, for each index `k`, we swap `a.(0)` with the value `v` in `a.(k)` and then move `v` down to its place in the heap.

for k = n-1 downto 1 do

let v = a.(k) in a.(k) <- a.(0); move_down a 0 v k
done

Note that the specification of `move_down` allows us to avoid assigning `v` to `a.(0)` before beginning the descent. Program 97 below contains the complete code.

## Complexity

Let us first consider the cost of a call to `move_down` a k v n. The number of recursive calls is bounded by $\log n$ since the value of k is doubled at each call. Furthermore, `move_down` performs at most two comparisons and an assignment in each call. Thus, in the worst case, we have a total of $2 \log n$ comparisons and $\log n$ assignments.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
Program 97 — Heapsort on arrays.

```pascal
let rec move_down a k v n =
  let r = 2 * k + 1 in
  if r >= n then
    a.(k) <- v
  else
    let rmax =
      if r+1 < n then if lt a.(r) a.(r+1) then r+1 else r
      else r in
      if le a.(rmax) v then a.(k) <- v
    else begin a.(k) <- a.(rmax); move_down a rmax v n end

let heapsort a =
  let n = Array.length a in
  for k = n/2 - 1 downto 0 do move_down a k a.(k) n done;
  for k = n-1 downto 1 do
    let v = a.(k) in a.(k) <- a.(0); move_down a 0 v k
  done
```

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
For the sorting itself, we can roughly bound the number $C(N)$ of comparisons in each call to `move_down` by $2 \log N$, that is, at worst a total of $3N \log N$ comparisons ($N \log N$ comparisons for the first step and $2N \log N$ for the second). Similarly, the total number of assignments is, at worst, $\frac{3}{2}N \log N$.

In fact, we can be more precise and show that the first step of the algorithm, namely, the construction of the heap, has linear cost. (See for example [7, Sec. 7.3].) It follows that only the second part of the algorithm contributes to the asymptotic complexity so that $C(N) \sim 2N \log N$. For an analysis on average of heapsort, the reader is referred to *The Art of Computer Programming* [15, vol. 3, p. 152].

Note that heapsort executes in constant memory. Indeed, the function `move_down` performs only tail-recursive calls and, moreover, all computations are performed in-place inside the array `a`.

## 12.5 Optimal Complexity

We present here a quick proof of the fact that the optimal complexity of a sorting algorithm that performs only element comparisons is in $O(N \log N)$.

We can visualize such an algorithm as a binary tree. Each internal node symbolizes a comparison that has been performed, the left (resp. right) subtree representing the rest of the algorithm when the test is positive (resp. negative). Each leaf constitutes a possible result, that is, a permutation performed on the initial sequence. If we assume $N$ distinct elements, there are $N!$ possible permutations and, hence, at least $N!$ leaves in this tree. Its height is consequently at least $\log N!$. Now the longest path from the root to a leaf corresponds to the largest number of comparisons performed by the algorithm on an input. There is, therefore, an input for which the number of comparisons is at least $\log N!$. By Stirling's formula, we know that $\log N! \sim N \log N$. For a more detailed proof, the reader is referred to *The Art of Computer Programming* [15, Sec. 5.3].

## 12.6 Experimental Evaluation

In this section, we compare the different sorting algorithms presented in this chapter empirically.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
## Evaluation Protocol

For a given list or array, we evaluate the performance of several algorithms on the same input as follows. We execute each sorting algorithm five times, measuring the computation time in each case. We exclude the smallest and the largest values, and take the average of the three remaining values.

For lists, we consider randomly constructed as well as already sorted lists, with lengths $2^i \times 1000$ for $i$ ranging from 0 to 10. For arrays, we consider randomly filled as already sorted arrays, with lengths $2^i \times 1000$ for $i$ ranging from 0 to 12. It is important to consider the case of data that is already sorted because that represents a realistic situation for which we wish to verify that our sorting algorithms are efficient.

## Evaluation on Lists

We compare the algorithms listed in figure 12.7. Figures 12.8 and 12.9 present the results.

For random lists (figure 12.8), we observe the following: Insertion sort (insertion) can only be used in case of small lists, given its quadratic complexity. (We do not give the computation time for lists longer than 16,000.) Two of the three quicksort algorithms, (quicksort and quicksort (rand)) triggered a stack overflow on lists of length $2^{10}$. The algorithm heapsort (trees), which uses binary trees, is penalized due to the construction of the trees, even if its complexity is optimal. The other sorting algorithms have comparable efficiency, the fastest being that of the OCaml standard library.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
<table><caption>Figure 12.7: Sorting algorithms on lists.</caption><tr><td>insertion</td><td>insertion sort of program 89</td></tr><tr><td>quicksort</td><td>quicksort of program 91</td></tr><tr><td>quicksort (log stack)</td><td>quicksort with the first call made to the smaller of the two lists after partitioning (see exercise 12.4)</td></tr><tr><td>quicksort (rand)</td><td>quicksort with a randomly chosen pivot (see exercise 12.5)</td></tr><tr><td>mergesort</td><td>merge sort of program 93</td></tr><tr><td>mergesort (reuse)</td><td>merge sort with optimization consisting in not constructing intermediate lists when partitioning (see exercise 12.10)</td></tr><tr><td>List.sort</td><td>the function List.sort of the OCaml standard library (version 3.11.2)</td></tr><tr><td>heapsort (arrays)</td><td>heapsort of program 96, using the structure based on arrays presented in section 6.3</td></tr><tr><td>heapsort (trees)</td><td>heapsort of program 96, using the structure based on binary trees presented in section 6.4</td></tr></table>

In case of already sorted lists (figure 12.9), we observe the following: Insertion sort (insertion) has quadratic complexity. (The favorable case would have been that of a list sorted in decreasing order.) The first two quicksort algorithms, quicksort and quicksort (log stack), have quadratic complexity, which is explained by the arbitrary choice of the first element of the list as pivot. (We do not give the computation time for lists of length greater than 16,000.) The third quicksort algorithm, quicksort (rand), behaves well but triggers a stack overflow on lists of size $2^{10} \times 1000$. The algorithm heapsort (trees) remains penalized due to the construction of the trees, as in case of random lists.

We may draw the following conclusions: Merge sort is the one that best combines simplicity and efficiency. Even its simplest version has an efficiency comparable to the optimized version of the OCaml standard library. Heapsort always performs well, but its implementation is considerably more complex than

— openly licensed via CC BY SA 4.0 —
that of merge sort. Finally, we observe that using quicksort on lists is not straightforward because we must, at a minimum, combine the optimizations of exercises 12.4 and 12.5, to avoid both stack overflows and the unfavorable case of quadratic complexity.

## Evaluation on Arrays

We now compare the different techniques for sorting arrays of integers based on algorithms listed in figure 12.10. Figures 12.11 and 12.12 present the results.

For random arrays (figure 12.11), we observe the following: As in case of lists, insertion sort can only be used for small arrays. The two heapsort algorithms, `heapsort` and `Array.sort`, are slightly less efficient than merge sort and quicksort. The two most efficient algorithms are both merge sorts.

For already sorted arrays (figure 12.12), we observe the following: Here, the first three quicksort algorithms, `quicksort`, `quicksort (insertion)`, and `quicksort (log stack)`, are in their most unfavorable case, which may be explained by the arbitrary choice of the left-most element as pivot. (We do not give the computation time for arrays of length greater than 16,000.) The version of quicksort with a randomly chosen pivot achieves a good efficiency, comparable to that of `mergesort`.

Of all the sorting algorithms in $O(N \log N)$, `heapsort` is among the least efficient while `merge sort` is among the most efficient. Here, `insertion sort` (`insertion`) is the most efficient, since we are in the favorable case of linear complexity. (The unfavorable case would have been that of an array sorted in decreasing order.)

We may draw the following conclusions: If we are allowed to allocate an intermediate array, in other words, to temporarily double memory use, then `merge sort` is the best choice. As in case of lists, it combines simplicity and efficiency. If, however, we wish to sort an array in-place, we must either use `heapsort` or `quicksort`, taking care to combine it with the optimizations of exercises 12.6 and 12.8, to avoid both stack overflows and quadratic behavior due to a bad choice of pivot.

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
12.7 Exercises

12.1 Write a function `two Way Sort: bool array -> unit` that sorts a boolean array in-place, performing only swaps of elements, with the convention that `false < true`. The complexity should be proportionate to the number of elements.

12.2 (Dijkstra's Dutch national flag problem) Write a function that sorts an array in-place. The array contains values representing the three colors of the Dutch flag, namely:

type color = Blue | White | Red

Only swaps of elements in the array may be performed. The complexity should be proportionate to the number of elements.

12.3 More generally, we consider the case of an array containing $k$ distinct values. To simplify matters, we assume that we are dealing with integers $0, \dots, k-1$. Write a function that sorts such an array in $O(\max(k, N))$ time, where $N$ is the size of the array.

Quicksort

12.4 One way to avoid a potential stack overflow in the function `quicksort` on lists is by first implementing a non-tail-recursive call to sort the smaller of the two lists, followed by a tail-recursive call to sort the larger one. Hint: For the tail-recursive call, generalize the function `quicksort` so that it takes as additional arguments the two lists `left` and `right`, so that `quicksort left 1 right` returns the concatenation, in this order, of `left`, the sorted list 1, and `right`. Show that the size of the stack remains logarithmic even in the worst case.

12.5 Modify the function `partition` on lists to use a randomly chosen element of the list as pivot. The aim is to avoid the unfavorable case of a sorted list, where quicksort would then have quadratic complexity.

12.6 One way of avoiding a potential stack overflow in the function `quick_rec` on arrays is to begin by performing the recursive call on the smaller of the two

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0/) —
portions. Show that the size of the stack remains logarithmic even in the worst case.

12.7 A classic idea to speed up sorting algorithms consists in switching to insertion sort when the number of elements to be sorted is small, that is, less than a constant fixed in advance (for example, 5). Modify the quicksort implementation on arrays taking this idea into consideration. Use the function `insertion_sort` on arrays (program 90) and generalize it by passing two indices, `l` and `r`, to delimit the portion of the array to be sorted.

12.8 Modify the quicksort implementation on arrays to use a pivot chosen randomly in the segment `a[l..r]`. As with lists, the idea is to avoid the unfavorable case of an already sorted array, where quicksort would then have quadratic complexity. One simple solution is to shuffle the array *before* sorting it, for example using exercise 2.11.

12.9 When numerous elements of the array are equal, the preceding exercise does not always suffice to guarantee a good choice of pivot. Modify the function `partition` so that it separates the elements that are strictly smaller than the pivot (to the left), those equal to the pivot (in the center), and those strictly greater than the pivot (to the right). Instead of two indices `m` and `i` dividing the segment of the array into three parts, as illustrated in figure 12.2.2, we will use three indices to divide the segment of the array into four parts. The new function `partition` must now return two indices. Modify the function `quick_rec` accordingly. Partitioning into three parts in this way is the subject of exercise 12.2.

## Merge Sort

12.10 One source of inefficiency of merge sort on lists is the function `split`. Given a list `l` containing `n` elements, `split l` constructs two completely new lists, `l1` and `l2`, containing all the elements of `l`. We can avoid this construction by using only the list `l`, without allocating any new lists. Rather than taking every other element of `l` to construct `l1` and `l2`, we choose for `l1` the `[n/2]` first elements of `l` and for `l2` the `[n/2]` last elements. To construct `l2`, write a function `chop: int -> 'a list -> 'a list` such that `chop n1 l` returns

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —
the suffix of 1 obtained by removing its n1 first elements. As for `l1`, it suffices to represent it as a prefix of 1, that is, to add an argument `n1` to the function `mergesort`, indicating that the first `n1` elements of 1 must be sorted. Rewrite the function `mergesort` using this idea. Reuse the function `merge` of program 93. Calculate the complexity in number of allocations of the function `mergesort`.

**12.11** One way of optimizing the function `mergesort` on almost sorted arrays is by avoiding the merging when, after both recursive calls, the elements of the left half are all smaller than the elements of the right half. We can easily test this by comparing the right-most element of the left half with the left-most element of the right half. Modify the function `mergesort` using this idea.

**12.12** To avoid copying `a` into `tmp` with `Array.blit` in the function `mergesort` on arrays, one possibility is to sort the two halves of the array `a` while simultaneously shifting them to the array `tmp`, then merging `tmp` back into `a` as is already done. However, sorting the elements of `a` into `tmp` requires, conversely, that the two halves be sorted in-place, then merged into `tmp`. We therefore need two mutually recursive sorting functions. We can make do with only one by passing an additional parameter indicating whether the sorting must be done in-place or in `tmp`. Modify the functions `mergesort` and `mergesort_rec` following this idea.

**12.13** In merge sort, as in quicksort, we can switch to insertion sort when the number of elements to be sorted is small (see exercise 12.7).
<table><thead><tr><th>N/1000</th><th>insertion</th><th>quicksort</th><th>quicksort (log stack)</th><th>quicksort (rand)</th><th>mergesort</th><th>mergesort (reuse)</th><th>List.sort</th><th>heapsort (arrays)</th><th>heapsort (trees)</th></tr></thead><tbody><tr><td>2<sup>4</sup></td><td>6.03</td><td>0.02</td><td><strong>0.01</strong></td><td>0.02</td><td>0.02</td><td><strong>0.01</strong></td><td><strong>0.01</strong></td><td><strong>0.01</strong></td><td>0.02</td></tr><tr><td>2<sup>5</sup></td><td>-</td><td>0.05</td><td><strong>0.03</strong></td><td>0.05</td><td>0.05</td><td><strong>0.03</strong></td><td><strong>0.03</strong></td><td>0.04</td><td>0.05</td></tr><tr><td>2<sup>6</sup></td><td>-</td><td>0.10</td><td>0.07</td><td>0.14</td><td>0.12</td><td>0.08</td><td><strong>0.05</strong></td><td>0.09</td><td>0.13</td></tr><tr><td>2<sup>7</sup></td><td>-</td><td>0.30</td><td>0.17</td><td>0.27</td><td>0.28</td><td>0.19</td><td><strong>0.12</strong></td><td>0.19</td><td>0.37</td></tr><tr><td>2<sup>8</sup></td><td>-</td><td>0.63</td><td>0.38</td><td>0.67</td><td>0.63</td><td>0.45</td><td><strong>0.27</strong></td><td>0.42</td><td>0.91</td></tr><tr><td>2<sup>9</sup></td><td>-</td><td>1.39</td><td>0.82</td><td>1.38</td><td>1.40</td><td>0.99</td><td><strong>0.65</strong></td><td>0.93</td><td>2.08</td></tr><tr><td>2<sup>10</sup></td><td>-</td><td>-</td><td>1.75</td><td>-</td><td>3.09</td><td>2.20</td><td><strong>1.41</strong></td><td>1.92</td><td>4.84</td></tr></tbody></table>

<figure><img src="image_124.png" /><figcaption>Figure 12.8: Comparison of sorting algorithms on random lists.</figcaption></figure>

— openly licensed via [CC BY SA 4.0](https://creativecommons.org/licenses/by-sa/4.0) —