A Short Introduction to the Art of Programming (EWD 316), Chapter 5

5. Programs corresponding to recurrence relations

Theorem 8 mentions successive states connected by a recurrence relation. The meaning of this theorem is twofold: it can be used to prove assertions about a given program, but also —and this, I think, is more important— it suggests to us, when faced with the task of making a program, the use of a while-clause in the case of a problem that in its mathematical formulation presents itself as the evaluation of a recurrence relation. We are going to illustrate this by a number of examples.

Consider the sequence of pairs a_i, c_i given by

for i = 0	a0 = 1 c0 = 1 – b, with 0 < b < 2 (i.e. abs(c0) < 1)	(1)
for i > 0	ai = (1 + ci-1) * ai-1 ci = (ci-1)2     .	(2)
Then	lim ai = 1/b     . i →∞

Exercise. Prove the last formula. (This has nothing to do with programming, it is secondary school algebra. The clue of a proof can be found in the relation

1		1 + ci-1
-------	=	-------	.)
1 – ci-1		1 - ci

It is requested to use this recurrence relation to approximate the value of 1/b; obviously we cannot compute infinitely many elements of the sequence a₀, a₁, a₂, … but we can accept a_k as a sufficiently close (how close?) approximation of 1/b when c_k is less (in absolute value) than a given, small, positive tolerance named "eps". (This example is of historical interest; it has been taken from the subroutine library for EDSAC 1, the world's first stored program controlled automatic computer. The order code of this computer did not comprise a divide instruction and one of the methods used with this computer to compute quotients was based on the above recurrence relation.)

Theorem 8 talks about "a part, s, of the state space" and the loop

while B(s) do s := S(s)

asserts that after the initial state s₀ , the states s_i after the i-th execution of the repeatable statement will satisfy

si = S(si-1)

(3)

Our recurrence relations (2) are exactly of the form (3) if we identify the state s_i with the value pair a_i, c_i. That is, to span the part s of the state space we have to introduce two variables, for the purpose of this discussion called A and C, and we shall denote their values after the i-th execution of the repeatable statement A_i and C_i respectively. We associate the state s_i (as given by the values A_i and C_i ) with the value pair a_i, c_i by the relations

Ai = ai Ci = ci

(4)

(Remember: on the left-hand sides the subscript "i" means "the value of the variable after the i-th execution of the repeatable statement", on the right-hand side the subscript "i" refers to the recurrent sequences as given by (1) and (2). It would have been usual to call the two variables "a" and "c" instead of "A" and "C", i.e. not to distinguish between the quantities defined in the mathematical formulation on the one hand and the associated variables in the program on the other hand. As this association is the very subject of this discussion, it would have been fatal not to distinguish between them.)

Within the scope of a declaration "real A, C" it is now a straightforward task to write the piece of program:

A := 1; C := 1 – b;
while abs(C ) ≥ eps do
begin A := (1 + C ) * A;
         C := C * C
end     .

The first line has to create the proper state s₀ and does so in accordance with (4) and (1), the repeatable statement has the form, symbolically denoted by "s := S(s)" —see the Note below— in accordance with (4) and (2), and the condition guarantees that after termination

(A_k = ) A = a_k

will hold with the proper value of k.

Exercise. Prove that the loop terminates.

Note. The symbolic assignment "s := S(s)" has the form of two assignments

A := (1 + C ) * A;
C := C * C     .

With the initial condition A = a_i-1 , C = c_i-1 the first assignment is equivalent to

A := (1 + c_i-1) * a_i-1

and after the first, but before the second assignment we have —on account of (2)—

A = a_i, C = c_i-1    .

We have the complete pair A = a_i, C = c_i only after the second assignment. Thanks to the explicit occurrence of the subscripts, the order of the two relations comprising (2) is immaterial, this in contrast to the two assignment statements composing the repeatable statement, whose order is vital.

Exercise. In the same EDSAC 1 subroutine library the next scheme is used. Consider the sequence of pairs a_i, c_i, given by

for i = 0	a0 = b c0 = 1 – b, with 0 < b < 2 (i.e. abs(c0) < 1)
for i > 0	ai = (1 + .5 * ci-1) * ai-1 ci = (ci-1)2 * (.75 + .25 * ci-1)     .
Then	lim ai = b.5     . i→∞

Prove the last formula and make a program using it for the approximation of the square root. The clue of a proof can be found in the relation

(1 – c_i-1)^-.5 = (1 + .5 * c_i-1) * (1 - c_i)^-.5     .

Prove also the termination of the repetition in the program made.

Exercise. In the same EDSAC 1 subroutine library the next scheme is used. Consider the sequence of triples inc_i, s_i, x_i, given by

for i = 0	inc0 = log 2 s0 = 0 x0 = arg               (with 1 ≤ arg < 2)
for i > 0
for (xi-1)2 < 2	inci = .5 * inci-1 si = si-1 xi = (xi-1)2     .
for (xi-1)2 ≥ 2	inci = .5 * inci-1 si = si-1 + .5 * inci-1 xi = .5 * (xi-1)2
Then	lim si = log(arg)     . i→∞

Prove this relation and make a program using it to approximate the logarithm of a value arg in the interval stated. (In this program "log 2" may be regarded as a known constant which, in fact, determines the base of the logarithm.) The clue of the proof can be found in the relation

log(arg) = s_i + inc_i * log(x_i ) / log 2     .

Our next example is very simple; it is so traditional that we could call it standard. (No self-respecting programming course omits it, it is often the very first example of a loop; Peter Naur uses it in his article "Proof of algorithms by general snapshots", 1966, BIT, 6, pp 310-316.)

Given a sequence of values

a[1], a[2], a[3], … , a[N]                  (with N ≥ 1)

and a variable called "max." Make a piece of program assigning to the variable named "max" the maximum value occurring in the sequence. (As N ≥ 1, the sequence is not empty and therefore the task makes sense; it is not required that any two values in the sequence differ from each other, the maximum value sought may occur more than once in the sequence.) If he welcomes the experience the reader is invited to try to make this piece of program himself before reading on.

How do we define the maximum value occurring in a sequence of length N for general N ≥ 1? If we call "maximum_k" the maximum value occurring among the first k elements a[1], … , a[k], then

1) the answer sought is maximum_N

2) the values maximum_k are given

for k = 1	by the base: maximum1 = a[1] appealing to the knowledge that the maximum element in a sequence of length 1 must be the only element in the sequence	(5)
for k > 1	by the recurrence relation: maximumk = MAX(maximumk-1, a[k]) assuming the knowledge of the function MAX of two arguments.	(6)

The recurrence relation (6) presents us with an additional difficulty because it is not of the form

s_i = S(s_i-1)

because —via "a[k]"— the value k occurs on the right-hand side not exclusively in the subscript "k-1". To overcome this we use a trick that might be called a method. If we call n_k the k-th natural number, then n_k = k; the numbers n_k satisfy the obvious recurrence relation n_k = 1 + n_k−1. We can now rewrite the definition for the sequence of values maximum_k in the form of a definition for the pairs n_k, maximum_k:

for k = 1	n1 = 1 maximum1 = a[1]	(7)
for k > 1	nk = 1 + nk-1 maximumk = MAX(maximumk-1, a[1 + nk -1])	(8)

and now the recurrence relations are of the form s_i = S(s_i-1), the only —trivial— difference being that in Theorem 8 we started with i = 0 and here with k = 1. The trick we called a method shows that we need a second (integer) variable; call it "m". Our state s_i will associate (with k = i + 1)

max_i = maximum_k
m_i = n_k     .

The piece of program now becomes:

max := a[1]; m := 1;
while m < N do begin m := m + 1;
                                   max := MAX(max, a[m])
                          end     .

Again, the order of the two assignment statements is essential.

We have given the above piece of reasoning and the explicit references to the recurrence relation of Theorem 8 because it shows a mechanism leading to the conclusion that the part of the state space on which the repetition operates needs to comprise an additional variable. Even a moderately trained programmer draws this conclusion "intuitively" and from now onwards I shall assume my reader equipped with that intuition. Then —and only then— there is a much shorter path of reasoning that leads to the program we found. It does not consider —"statically" so to speak— the sequence of values s₀, s₁, … in the sense that it bothers about the values of the subscript i in s_i. It appeals directly to Theorems 5 and 6 and works in terms of assertions valid (before and after) any execution of the repeatable statement. The price paid for this is the duty to prove termination separately.

Given the base

k = 1

maximum1 = a[1]

and the step

1 < k ≤ N

maximumk = MAX(maximumk -1, a[k])

the programmer "intuitively" introduces two variables which he calls "maximum" and "k" for short and the relation to be kept invariant is

P :	1 ≤ k ≤ N   and   maximum = maximumk     .

(Here the use of "maximum" and "k" stands for the current value of the variable thus named, while "maximum_k" stands for the value as given by the recurrence relation. This double use of the same names is tricky but programmers do it. I too.)

The program then consists of two parts: establishing the relation P in accordance with the base and repeatedly increasing k under invariance of relation P, i.e. in accordance with the step.

The initialization

"maximum := a[1]; k := 1"

establishes P (with k = 1), the repetition

while k < N do
begin k := k + 1;
          maximum := MAX(maximum, a[k])
end

causes the repeatable statement to be executed under the combined relation "B and P ", i.e.

k < N   and   1 ≤ k ≤ N   and   maximum = maximum_k

which reduces to

1 ≤ k < N   and   maximum = maximumk     .

(9)

In order to show that the execution of the repeatable statement under the initial condition (9) leaves relation P valid, it is desirable to distinguish between the values before and after its execution; now it would be confusing to do so with subscripts (why?), therefore we distinguish the values after execution by primes.

Initially we have relation (9); after the assignment k := k + 1 we have the relation k' = k + 1 and from the first part of (9), i.e. 1 ≤ k ≤ N, follows 2 ≤ k' ≤ N, which implies

1 ≤ k' ≤ N     .

(10)

The second assignment now becomes effectively maximum := MAX(maximum_k, a[k']), resulting in the relation

maximum' = maximumk'     .

(11)

Relations (10) and (11) combine to a replica of P, but now for the primed quantities.

Termination follows from the fact that each execution of the repeatable statement involves an effective increase of the integer valued variable k. After termination we have, according to Theorem 5, "P and non B ", i.e.

1 ≤ k ≤ N   and   maximum = maximum_k    and   non k < N ;

from the first and the last term we conclude k = N and then from the middle part

maximum = maximum_N

which concludes the proof.

Exercise. Make a program effectively assigning "prod := X * Y" with integer X and Y, satisfying X ≥ 0, Y ≥ 0
a) using only addition and subtraction
b) using in addition the boolean function "odd(x)", doubling and halving of a number. (The so-called Egyptian multiplication.)

Exercise. Make a program effectively assigning "rem :=REM(X, Y)" with integer X and Y, X ≥ 0, Y > 0, where the function REM(X, Y) is the remainder after the division of X by Y
a) using only addition and subtraction
b) using in addition doubling and halving of a number. Modify both programs in such a way that in addition "quot :=QUOT(X, Y)" will take place. (The so-called Chinese division.)

We conclude this section by an example of the (standard) circumstance in which a recurrence relation should not be translated blindly into a loop. Given two sequences of values

x[1], x[2], … , x[N]        and
y[1], y[2], … , y[N]                    with N ≥ 0     ;

make a program assigning to the boolean variable named "eq" the value true if x[i] = y[i] for all i satisfying 1 ≤ i ≤ N and the value false if in that range a value for i exists such that x[i] ≠ y[i]. (The sequence may be empty, in that case "eq" should get the value true.)

How do we define equality for sequences of length N for general N? Again by means of a recurrence relation. Let eq_i mean "no difference occurs among the first i pairs"; the sequence of values eq_i is given by

for i = 0	eq0 = true
for i > 0	eqi = eqi-1 and x[i] = y[i]     .

The net effect of the program to be made should be eq := eq_N.

A blind translation into initialization followed by repetition would lead to

eq := true; i := 0;
while i < N do begin i := i + 1; eq := (eq and x[i] = y[i]) end

Although the above program is correct, the following program, besides being equally correct, is on the average more efficient:

eq := true; i := 0;
while i < N and eq do begin i := i + 1; eq := (x[i] = y[i]) end     .

because it terminates the repetition as soon as a difference has been found.

Exercise. Prove the correctness of the second program.

Back to top

Next chapter: 6. A first example of step-wise program composition

---