In words, given an array of integers having at least two elements, compute the sum of squares of the difference between all pairs of elements. (Following the convention of the guarded command language, function application is written f.x, and an array is seen as a function from indices to values.)

It is not hard to quickly write up a O(N²) program using nested loops, which, I have to confess, is what I would do before reading Kaldewaij’s book and realised that it is possible to do the task in linear time using one loop. Unfortunately, not many students managed to come up with this solution, therefore I think it is worth some discussion.

Quantifiers

Before we solve the problem, let us review the “Dutch style” quantifier syntax and rules. Given a commutative, associative binary operator ⊕ with unit element e, if we informally denote the (integral) values in the interval [A .. B) by i₀, i₁, i₂ ... in, the quantified expression:

denotes F.i₀ ⊕ F.i₁ ⊕ F.i₂ ⊕ ... ⊕ F.in. We omit the i in R.i and F.i when there can be no confusion.

A more formal characterisation of the quantified expression is given by the following rules:

(⊕ i : false : F.i) = e

(⊕ i : i = x : F.i) = F.x

(⊕ i : R : F) ⊕ (⊕ i : S : F) = (⊕ i : R ∨ S : F) ⊕ (⊕ i : R ∧ S : F)

(⊕ i : R : F) ⊕ (⊕ i : R : G) = (⊕ : R : F ⊕ G)

(⊕ i : R.i : (⊕ j : S.j : F.i.j)) = (⊕ j : S.j : (⊕ i : R.i : F.i.j))

Rules 1 and 3 give rise to a useful rule "split off n": consider i such that 0 ≤ i < n + 1. If n > 0, the set of possible values of i can be split into two subsets: 0 ≤ i < n and i = n. By rule 3 (reversed) and 1 we get:

(⊕ i : 0 ≤ i < n + 1 : F.i) = (⊕ i : 0 ≤ i < n : F.i) ⊕ F.n

Expressions quantifying more than one variables can be expressed in terms of quantifiers over single variables:

(⊕ i,j : R.i ∧ S.i,j : F.i.j) = (⊕ i : R.i : (⊕ j : S.i.j : F.i.j))

If ⊗ distributes into ⊕, we have an additional property:

x ⊗ (⊕ i : R : F) = (⊗ i : R : x ⊗ F)

As a convention, (+ i : R : F) is often written (Σ i : R : F).

Computing the Sum of Squares of Differences

The first step is to turn the constant N to a variable n. The main worker of the program is going to be a loop, in whose invariant we try to maintain:

P ≣ r = (Σ i,j : 0 ≤ i < j < n : (a.i - a.j)²)

In the end of the loop we increment n, and the loop terminates when n coincides with N:

{ Inv: P ∧ 2 ≤ n ≤ N , Bound: N - n}
do n ≠ N → ... ; n := n + 1 od

We shall then find out how to update r before n := n + 1 in a way that preserves P.

Assume that P and 2 ≤ n ≤ N holds. To find out how to update s, we substitute n for n + 1 in the desired value of r:

This is where most people stop the calculation and start constructing a loop computing (Σ i : 0 ≤ i < n : (a.i - a.n)²). One might later realise, however, that most computations are repeated. Indeed, the expression above can be expanded further:

Another “One Loop” Solution

Among those students who did come up with a program, most of them resorted to a typical two-loop, O(N²) solution. Given that this 9-hour course is, for almost all of them, their first exposure to program derivation, I shall perhaps be happy enough that around 3 to 4 out of 38 students came up with something like the program above.

The program uses only one loop, but is still O(N²) — on a closer inspection one realises that it is actually simulating the inner loop manually. Still, I’d be happy if the student could show me a correctness proof, with a correct loop invariant and a bound, since both of them are more complex than what I expected them to learn. Unfortunately, in the answer handed in, the program, the invariant, and the bound all contain some bugs. Anyone wants to give it a try?

If you think you know everything you need to know about binary search, but have not read Netty van Gasteren and Wim Feijen’s note The Binary Search Revisited, you should.

In the FLOLAC summer school this year we plan to teach the students some basics of Hoare logic and Dijkstra’s weakest precondition calculus — a topic surprisingly little known among Taiwanese students. Binary search seems to be a topic I should cover. A precise specification will be given later, but basically the requirement is like this: given a sorted array of N numbers and a value to search for, either locate the position where the value resides in the array, or report that the value does not present in the array, in O(log N) time.

Given that everybody should have learned about binary search in their algorithm class, you would not expect it to be a hard programming task. Yet in his popular book Programming Pearls, Jon Bentley noted that surprisingly few professions programmers managed to implement the algorithm without bugs at their first attempt. I tried it myself and, being very careful (and having learned something about program construction already), I produced a program which I believe to be correct, but rather inelegant. I tried it on my students and they did produce code with some typical bugs. So it seems to be a good example for a class about correct program construction.

The Van Gasteren-Feijen Approach

The surprising fact that van Gasteren and Feijen pointed out was that binary search does not apply only to sorted lists! In fact, the usual practice comparing binary search to searching for a word in a dictionary is, according to van Gasteren and Feijen, a major educational blunder.

Van Gasteren and Feijen considered solving a more general problem: let M and N be integral numbers with M < N, and let Φ be a relation such that M Φ N (where Φ is written infix), with some additional constraints to be given later. The task is to find l such that

Notice first that the loop guard l+1 ≠ r, if satisfied, guarantees that l and r are not adjacent numbers, therefore assigning m := (l + r) / 2 establishes l < m < r, and thus the bound r - l is guaranteed to decrease. The if statement clearly maintains the invariant, if at least one of the guards are always satisfied:

l Φ r ∧ l < m < r ⇒ l Φ m ∨ m Φ r(*)

For Φ satisfying the condition above, at the end of the loop we will find some l such that l Φ (l+1).

What relations satisfy (*)? Examples given by van Gasteren and Feijen include:

i Φ j = a[i] ≠ a[j] for some array a. Van Gasteren and Feijen suggested using this as the example when introducing binary search.

i Φ j = a[i] < a[j],

i Φ j = a[i] × a[j] ≤ 0,

i Φ j = a[i] ∨ a[j], etc.

Searching for a Key in a Sorted Array

To search for a key K in an ascending-sorted array a, it seems that we could just pick this Φ:

i Φ j = a[i] ≤ K < a[j]

and check whether a[i] = K after the loop. There is only one problem, however -- we are not sure we can establish the precondition a[l] ≤ K < a[r]!

Van Gasteren and Feijen's solution is to add to imaginary elements to the two ends of the array. That is, for a possibly empty array a[0..N), we imagine two elements a[-1] such that a[-1] ≤ x and a[N] such that x < a[N] for all x. I believe this is equivalent to using this Φ:

Do not worry about the idea "adding" elements to a. The invariant implies that -1 < m < N, thus a[-1] and a[N] are never accessed, and the array a needs not be actually altered. They are just there to justify the correctness of the program. It also enables us to handle possibly empty arrays, while the loop body seems to be designed for the case when the range [l..r) is non-empty.

I would like to be able to derive this program in class, since this appears to be the more popular version. Apart from the presence of break, which I do not yet know of a easy variation of Hoare logic that helps to derive it, to relate the test a[m] < K to l := m + 1 I will have to bring in the fact that a is sorted in an earlier stage of the development. Thus it is harder to put it in a more general picture.

For several reasons I used to believe that Bentley's program could be preferred, for example, it seems to shrink the range more effectively, assigning l and r to m + 1 and m - 1, rather than m. On a second thought I realised that it might not be true. Variable l can be assigned m + 1 because the possibility of a[m] = K is covered in another case with an early exit, and r is assigned m - 1 because this algorithm represents an array segment with an inclusive right bound, as opposed to the previous algorithm.

The two algorithms do not solve exactly the same problem. With multiple occurrences of K in the array, Bentley's algorithm is non-deterministic about which index it returns, while the van Gasteren-Feijen algorithm, enforced by the specification, always returns the largest index. When K does not appear in the array, van Gasteren and Feijen's program could be more efficient because it needs only one comparison in the loop, rather than two as in Bentley's case (I am assuming that the last comparison is a catch-all case and need not be implemented). What if K does present in the array? An analysis by Timothy J. Rolfe concluded that a single-comparison approach is still preferable in average -- benefit of the early exit does not outweigh the cost of the extra comparison in the loop.

On Computing the Middle Index

There are some other interesting stories regarding the assignment m := (l + r) / 2. Joshua Bloch from Google noted that for large arrays, adding l and r may cause an overflow, and Bloch was not picking on Bentley -- the bug was reported by Sun. Bloch suggests one of the following:

Since the publication of the blog post, there have been numerous discussions on whether it should be considered a bug in the binary search algorithm or the integer datatype, and some more machine dependent issues like whether one may have an array so large that cannot be indexed by an int, etc.

Exercise?

Among the exercises suggested by van Gasteren and Feijen, this one caught my eye: let array a[0,N), with 0 < N, be the concatenation of a strictly increasing and a strictly decreasing array. Use binary search to find the maximum element. (For this problem I think it is reasonable to assume that the two sub-arrays could be empty, while a is non-empty.) This happens to a sub-routine I needed for an O(N log N) algorithm for the maximum segment density problem (there are linear-time algorithms for this problem, though), and I do remember I started off treating it as an unimportant sub-routine but had a hard time getting it right. I am glad that now I know more about it.