Sum of Squares of Differences

In the final exam of the Program Construction course in FLOLAC ’10, I gave the students this problem (from Kaldewaij’s book):

|[ con N {N ≥ 2}; a : array [0..N) of int;
   var r : int;
   S
   { r = (Σ i,j : 0 ≤ i < j < N : (a.i - a.j)²) }
]|

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:

   (⊕ i : A ≤ i < B : F.i)

informally denotes F.i₀ ⊕ F.i₁ ⊕ F.i₂ ⊕ ... ⊕ F.in. More generally, if all values satisfying predicate R can be enlisted i₀, i₁, i₂ ... in, the expression

   (⊕ i : R.i : F.i)

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:

1. (⊕ i : false : F.i) = e
2. (⊕ i : i = x : F.i) = F.x
3. (⊕ i : R : F) ⊕ (⊕ i : S : F) = (⊕ i : R ∨ S : F) ⊕ (⊕ i : R ∧ S : F)
4. (⊕ i : R : F) ⊕ (⊕ i : R : G) = (⊕ : R : F ⊕ G)
5. (⊕ 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:

   (Σ i,j : 0 ≤ i < j < n : (a.i - a.j)²)[n+1 / n]
 = (Σ i,j : 0 ≤ i < j < n + 1 : (a.i - a.j)²)
 =   { split off j = n }
   (Σ i,j : 0 ≤ i < j < n : (a.i - a.j)²) + 
   (Σ i : 0 ≤ i < n : (a.i - a.n)²)
 =   { P }
   r + (Σ i : 0 ≤ i < n : (a.i - a.n)²)

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:

   r + (Σ i : 0 ≤ i < n : (a.i - a.n)²)
 =   { (x - y)² = x² - 2xy + y² }
   r + (Σ i : 0 ≤ i < n : a.i² - 2 × a.i × a.n + a.n²)
 =   { Rule 4 }
   r + (Σ i : 0 ≤ i < n : a.i²) 
     - (Σ i : 0 ≤ i < n : 2 × a.i × a.n) 
     + (Σ i : 0 ≤ i < n : a.n²)
 =   { a.n is a constant, multiplication distributes into addition }
   r + (Σ i : 0 ≤ i < n : a.i²) 
     - 2 × (Σ i : 0 ≤ i < n : a.i) × a.n 
     + (Σ i : 0 ≤ i < n : a.n²)
 =   { simplifying the last term }
   r + (Σ i : 0 ≤ i < n : a.i²) 
     - 2 × (Σ i : 0 ≤ i < n : a.i) × a.n + n × a.n²

which hints at that we can store the values of (Σ i : 0 ≤ i < n : a.i²) and (Σ i : 0 ≤ i < n : a.i) in two additional variables:

  Q₀  ≣  s = (Σ i : 0 ≤ i < n : a.i²)
  Q₁  ≣  t = (Σ i : 0 ≤ i < n : a.i)

It merely takes some routine calculation to find out how to update s and t. The resulting code is:

|[ con N {N ≥ 2}; a : array [0..N) of int;
   var r, s, t, n : int;

   r, s, t, n := (a.0 - a.1)², a.0² + a.1², a.0 + a.1, 2
   { Inv: P ∧ Q₀ ∧ Q₁ ∧ 2 ≤ n ≤ N , Bound: N - n }
 ; do n ≠ N → 
      r := r + s - 2 × t × a.n + n × a.n²;
      s := s + a.n²;
      t := t + a.n;
      n := n + 1 
   od
   { r = (Σ i,j : 0 ≤ i < j < N : (a.i - a.j)²) }
]|

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.

One student, however, delivered a program I did not expect to see:

|[ con N {N ≥ 2}; a : array [0..N) of int;
   var r, i, j : int;

   r, i, j := 0, 0, 0
   { Inv: ... ∧ 0 ≤ i ≤ j ∧ 0 ≤ j ≤ N, Bound: ? }
 ; do j ≠ N →
       if i < j → r := r + (a.i - a.j)²;  i := i + 1
        | i = j → i, j := 0, j + 1
       fi
   od
]|

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?