There's this great little book, A Method of Programming by Edsger W. Dijkstra and Wim Feijen, which is dry and precise as a smack of a ruler against collected fingertips. Really, I cannot recommend it highly enough. It includes a little essay, Coordinate transformation
that —for me— illustrates beautifully the methodical, stepwise refining of a program by meaning-preserving transformations.
They start with a naïve version of code for extracting integer square roots:
let isqrt n =
if n <= 0 then failwith "isqrt" else
let a = ref 0
and b = ref 1 in
while !b * !b <= n do b := !b lsl 1 done;
while !b != !a + 1 do
let c = (!a + !b) lsr 1 in
if !c * !c <= n then a := !c else b := !c
done;
!a
The first loop finds the least k such that b² = 4k > n; this takes 1 + ⌊(lg n)/2⌋ iterations. The second iteration performs a binary search for the least c such that c² ≤ n < (c+1)² in the range 1 to b, for 1 + ⌊lg (b - a)⌋ = 1 + ⌊lg (2k)⌋ = 1 + ⌊lg n⌋ iterations. But they note that "each iteration requires a multiplication, however, and since on the whole the general multiplicative operations requires (an order of magnitude) more time than additive operations, halvings and doublings, the question arises whether these general multiplications can be eliminated." To answering this in the affirmative then they set out.
The basic idea is to replace multiplications (more specifically, squarings) by finite differences. It is well-known that the sum of integers in arithmetic progression gives (more or less) a square. Note that b - a is always a power of 2 (by induction), and so the division by shifting is exact. This enables them to introduce a variable d = b - a at every step:
let isqrt n =
if n <= 0 then failwith "isqrt" else
let a = ref 0
and b = ref 1
and d = ref 1 (* invariant: d = b - a *) in
while !b * !b <= n do b := !b lsl 1; d := !d lsl 1 done;
while !b != !a + 1 do
let c = (!a + !b) lsr 1 in
if !c * !c <= n then begin
a := !c;
d := !d lsr 1
end else begin
b := !c;
d := !d lsr 1
end
done;
!a
which readily transforms to:
let isqrt n =
if n <= 0 then failwith "isqrt" else
let a = ref 0
and b = ref 1
and d = ref 1 (* invariant: d = b - a *) in
while !b * !b <= n do b := !b lsl 1; d := !d lsl 1 done;
while !b != !a + 1 do
let c = (!a + !b) lsr 1 in
d := !d lsr 1;
(* a + d = c ∧ c + d = b *)
if !c * !c <= n then a := !c else b := !c
done;
!a
Moving the update to d upwards out of the conditional instead of downards might seem a rabbit pulled out of the hat, but by doing so all four cs can be replaced by a + d in the last alternative statement, as the invariant in the comment explains. This makes c unused and so it can be eliminated:
let isqrt n =
if n <= 0 then failwith "isqrt" else
let a = ref 0
and b = ref 1
and d = ref 1 (* invariant: d = b - a *) in
while !b * !b <= n do b := !b lsl 1; d := !d lsl 1 done;
while !b != !a + 1 do
d := !d lsr 1;
if (!a + !d) * (!a + !d) <= n then a := !a + !d else b := !a + !d
done;
!a
Now in the first loop, a = 0, and so b = d, hence the loop condition can mention d instead of b. Also, note that because d = b - a is kept invariant in the second repetition, the condition b ≠ a + 1 can be rewritten as d ≠ 1, and the second branch of the conditional is trivial. This makes b unused and can be eliminated:
let isqrt n =
if n <= 0 then failwith "isqrt" else
let a = ref 0
and d = ref 1 in
while !d * !d <= n do d := !d lsl 1 done;
while !d != 1 do
d := !d lsr 1;
if (!a + !d) * (!a + !d) <= n then a := !a + !d
done;
!a
This doesn't really bought us anything except replacing two variables by one. However, the square of the binomial in the condition inside the loop can be rewritten as (a + d)² ≤ n ≡ 2⋅a⋅d + d² ≤ n - a². Now, we could keep three variables, p = a⋅d, q = d² and r = n - a² invariantly (hm, invariant variables…) and profit from the simpler relations holding between them:
let isqrt n =
if n <= 0 then failwith "isqrt" else
let a = ref 0
and d = ref 1
and p = ref 0
and q = ref 1
and r = ref n in
while !d * !d <= n do q := !q lsl 2; d := !d lsl 1 done;
while !d != 1 do
q := !q lsr 2;
p := !p lsr 1;
d := !d lsr 1;
if 2 * !a * !d + !d * !d <= n - !a * !a then begin
r := !r - (2 * !a * !d + !d * !d);
p := !p + !q;
a := !a + !d
end
done;
!a
The introduction was straightforward and mechanical, except for the fact that, since a = 0 in the first loop, doubling p was not necessary. Now observe that in it the condition d² ≤ n is equal to q ≤ r since a = 0, and so we may effect the replacement. On the other hand, in the second loop the condition d ≠ 1 equivales to q ≠ 1, and we may replace it. Finally, in the inner conditional, 2⋅a⋅d + d² = 2⋅p + q and n - a² = r by definition, and we can simplify. These substitutions make both d and a unused and they can be left out, since with d = 1 is p = a:
let isqrt n =
if n <= 0 then failwith "isqrt" else
let p = ref 0
and q = ref 1
and r = ref n in
while !q <= !r do q := !q lsl 2 done;
while !q != 1 do
q := !q lsr 2;
p := !p lsr 1;
if 2 * !p + !q <= !r then begin
r := !r - (2 * !p + !q);
p := !p + !q
end
done;
!p
For convenience we could introduce a new constant h = 2⋅p + q:
let isqrt n =
if n <= 0 then failwith "isqrt" else
let p = ref 0
and q = ref 1
and r = ref n in
while !q <= !r do q := !q lsl 2 done;
while !q != 1 do
q := !q lsr 2;
p := !p lsr 1;
let h = 2 * !p + !q in
if h <= !r then begin
r := !r - h;
p := !p + !q
end
done;
!p
Dijkstra's and Feijen's essay stops here, but not before they remark that:
Here we have shown extensively the transformation process — introduction of new variables and invariants, rewriting guards, and then elimination of original variables. Carried out this way it requires a lot of writing. Someone who is used to this method of program development would make the transition from old to new variables in one step, without writing down the version with both variables. Rewriting, which this method of program development entails, remains a drawback, however.
To which I would add, not so much of a drawback with a good text editor with multiple undos and an interactive top-level interpreter with a read-eval-print loop.
No comments:
Post a Comment