Twenty years’ analysis of the Binary Euclidean Algorithm^{∗}

Dedicated to Tony Hoare on the occasion of his departure from Oxford

Richard P. Brent

Computing Laboratory Oxford University

15 September 1999

∗Presented at a Symposium in Celebration of the Work of Tony Hoare, Oxford, 13–15 September 1999.

Copyright \begin{matrix}c \\ ◯ \end{matrix}1999, R. P. Brent.  rpb183t

Abstract

The binary Euclidean algorithm is a variant of the classical Euclidean algorithm. It avoids divisions and multiplications, except by powers of two, so is potentially faster than the classical algorithm on a binary machine.

I will describe the binary algorithm and consider its average case behaviour. In particular, I will describe the recent discovery of an error which I made in 1976, discuss some recent results of Brigitte Vallée, and describe a numerical computation which verifies a conjecture of Vallée to 40 decimal places.

Outline

• Introduction

• The binary Euclidean algorithm – First formulation (Algorithm B)

– History

• Continuous model – Conjectured/empirical results

• Another formulation (Algorithm V)

• Useful operators – Relation between the operators – Fixed points

– Vallée’s conjectures and results • Numerical results

• Tracking down an error

• Moral

• Conclusion and open problems

Introduction – Quicksort etc

Why a paper on average-case analysis of algorithms ?

Tony Hoare and I both have an interest in Euclidean algorithms, and, since the publication of his 1962 paper on Quicksort, Tony can hardly deny an interest in the average-case analysis of algorithms. Although Tony proved the correctness of Quicksort in his usual thorough and convincing manner, the real significance of Quicksort is its excellent average case behaviour. Tony was well aware of this: a large part of his 1962 paper is devoted to the average case analysis of quicksort.

Personally, I prefer Heapsort, since (unlike Quicksort) it has guaranteed worst case time bound O(n log n) to sort n items. After the Proceedings paper was written, it emerged in a tea-room conversation that Tony was closely involved in the discovery of Heapsort. He and Williams were in the same room; Williams found heapsort then Tony found a quicker sort!

The Topic ⋯

In any case, given Tony’s interest in average case analysis, my conclusion was that I should consider the average case behaviour of Euclidean algorithms. Fortunately, thanks to Don Knuth and Brigitte Vallée, I had something (relatively) new to say about the binary (not the classical) Euclidean algorithm.

Notation

lg(x) denotes log_{2}(x).

N, n, u, v are positive integers.

Val_{2}(u) denotes the dyadic valuation of the positive integer u, i.e. the greatest integer j such that 2^{j} | u.

The Binary Euclidean Algorithm

The idea of the binary Euclidean algorithm is to avoid the “division” operation r ← m mod n of the classical algorithm, but retain O(log N) worst (and average) case.

We assume that the algorithm is implemented on a binary computer so division by a power of two is easy. In particular, we assume that the “shift right until odd” operation

u ← u/2^{Val_{2}(u)}

or equivalently

while even(u) do u ← u/2

can be performed in constant time, although time O(Val_{2}(u)) would be sufficient.

Definition

There are several almost equivalent ways to define the algorithm. It is easy to take account of the largest power of two dividing the inputs, so for simplicity we assume that u and v are odd positive integers.

Following is a simplified version of the algorithm given in Knuth, §4.5.2.

Algorithm B

B1. t ← |u − v|; if t = 0 terminate with result u

B2. t ← t/2^{Val_{2}(t)}

B3. if u ≥ v then u ← t else v ← t; go to B1.

History

The binary Euclidean algorithm is attributed to Silver and Terzian (unpublished, 1962) and Stein (1967). However, it seems to go back almost as far as the classical Euclidean algorithm. Knuth (§4.5.2) quotes a translation of a first-century AD Chinese text Chiu Chang Suan Shu on how to reduce a fraction to lowest terms:

If halving is possible, take half.

Otherwise write down the

denominator and the numerator, and subtract the smaller from the greater.

Repeat until both numbers are equal.

Simplify with this common value.

This looks very much like Algorithm B !

A Heuristic Continuous Model

To analyse the expected behaviour of Algorithm B, we can follow what Gauss did for the classical algorithm. This was first attempted in my 1976 paper^{1} and there is a summary in Knuth (Vol. 2, third edition, §4.5.2).

Assume that the initial inputs u_{0}, v_{0} to Algorithm B are uniformly and independently distributed in (0, N), apart from the restriction that they are odd. Let (u_{n}, v_{n}) be the value of (u, v) after n iterations of step B3.

Let

x_{n} = \frac{min(u_{n}, v_{n})}{max(u_{n}, v_{n})}

and let F_{n}(x) be the probability distribution function of x_{n} (in the limit as N → ∞). Thus F_{0}(x) = x for x ∈ [0, 1].

1R. P. Brent, Analysis of the Binary Euclidean Algorithm, New Directions and Recent Results in Algorithms and Complexity, (J. F. Traub, editor), Academic Press, New York, 1976, 321–355.

Plausible Assumption

We make the plausible (but not proved) assumption^{2} that Val_{2}(t) takes the value k with probability 2^{−k} at step B2.

It is plausible because Val_{2}(t) at step B2 depends on the least significant bits of u and v, whereas the comparison at step B3 depends on the most significant bits, so one would expect the steps to be (almost) independent.

2Vallée does not need to make this assumption. Her results will be mentioned later.

The Recurrence for F_{n}

Consider the effect of steps B2 and B3. We can assume that u > v so t = u − v. If Val_{2}(t) = k then X = v/u is transformed to

X′ = min (\frac{u − v}{2^{k}v} , \frac{2^{k}v}{u − v} )

= min (\frac{1 − X}{2^{k}X} , \frac{2^{k}X}{1 − X} ) .

It follows that X′ < x iff

X < \frac{1}{1 + 2^{k}/x} or X > \frac{1}{1 + 2^{k}x} .

Thus, the recurrence for F̃_{n}(x) = 1 − F_{n}(x) is

F̃_{n+1}(x) =

∑_{k≥1} 2^{−k} (F̃_{n} (\frac{1}{1 + 2^{k}/x}) − F̃_{n} (\frac{1}{1 + 2^{k}x}))

and F̃_{0}(x) = 1 − x for x ∈ [0, 1].

The Recurrence for f_{n}

Differentiating the recurrence for F̃_{n} we obtain (formally) a recurrence for the probability density f_{n}(x) = F′_{n}(x) = −F̃′_{n}(x):

f_{n+1}(x) = ∑_{k≥1} (\frac{1}{x + 2^{k}} )^{2} f_{n} (\frac{x}{x + 2^{k}} )

+ ∑_{k≥1} (\frac{1}{1 + 2^{k}x})^{2} f_{n} (\frac{1}{1 + 2^{k}x}) .

Operator Notation

The recurrence for f_{n} may be written as

f_{n+1} = B_{2}f_{n},

where the operator B_{2} is the case s = 2 of a more general operator B_{s} which will be defined later.

Conjectured and Empirical Results

In my 1976 paper I gave numerical and analytic evidence that F_{n}(x) converges to a limiting distribution F (x) as n → ∞, and that f_{n}(x) converges to the corresponding probability density f(x) = F′(x) (note that f = B_{2}f so f is a “fixed point” of the operator B_{2}).

Assuming the existence of F, it is shown in my 1976 paper that the expected number of iterations of Algorithm B is ∼ K lg N as N → ∞, where K = 0.705 … is a constant defined by

K = ln 2/E_{∞} ,

and

E_{∞} = ln 2 +

∫_{0}^{1} (∑_{k=2}^{∞} (\frac{1 − 2^{−k}}{1 + (2^{k} − 1)x}) − \frac{1}{2(1 + x)}) F (x) dx .

A Simplification

We can simplify the expression for K to obtain

K = 2/b ,

where

b = 2 − ∫_{0}^{1} lg(1 − x)f(x) dx .

Using integration by parts we obtain an equivalent expression

b = 2 + \frac{1}{ln 2} ∫_{0}^{1} \frac{1 − F (x)}{1 − x} dx .  (1)

For my direct proof of (1), see Knuth, third edition, §4.5.2.

Another Formulation – Algorithm V

It will be useful to rewrite Algorithm B in the following equivalent form (using pseudo-Pascal):

Algorithm V { Assume u ≤ v }

while u ≠ v do begin while u < v do

begin

j ← Val_{2}(v − u); v ← (v − u)/2^{j}; end;

u ↔ v;

end; return u.

Continued Fractions

Vallée (Algorithmica, 1998) shows a connection between Algorithm V and continued fractions of a certain form:

\frac{u}{v} = 1/a_{1} + 2^{k_{1}}/a_{2} + 2^{k_{2}}/ … /a_{r} + 2^{k_{r}} ,

where a_{j} is odd, k_{j} > 0, and 0 < a_{j} < 2^{k_{j}}.

Some Useful Operators

Operators B_{s}, U_{s}, Ũ_{s}, V_{s}, useful in the analysis of the binary Euclidean algorithm, are defined on suitable function spaces by

U_{s}[f](x) = ∑_{k≥1} (\frac{1}{1 + 2^{k}x})^{s} f (\frac{1}{1 + 2^{k}x}) , (2)

Ũ_{s}[f](x) = (\frac{1}{x})^{s} U_{s}[f] (\frac{1}{x}) ,  (3)
B_{s} = U_{s} + Ũ_{s},

V_{s}[f](x) = ∑_{k≥1} ∑_{{a odd, \atop 0<a<2^{k}}} (\frac{1}{a + 2^{k}x})^{s} f (\frac{1}{a + 2^{k}x}) .

(4)

In these definitions s is a complex variable, and the operators are called Ruelle operators. They are linear operators acting on certain function spaces.

The case s = 2 is of particular interest. B_{2} encodes the effect of one iteration of the inner “while” loop of Algorithm V, and V_{2} encodes the effect of one iteration of the outer “while” loop.

History and Notation

B_{2} (denoted T) was introduced in my 1976 paper and was generalised to B_{s} by Vallée. V_{s} was introduced by Vallée. We shall call

• B_{2} the binary Euclidean operator and

• V_{s} Vallée’s operator.

Relation Between the Operators

The operators are closely related, as the following results show.

Lemma 1

V_{s} = V_{s}Ũ_{s} + U_{s}.

An algebraic proof of Lemma 1 is given in the Proceedings.

Algorithmic interpretation

Algorithm V gives an interpretation of Lemma 1 in the case s = 2. If the input density of x = u/v is f(x) then execution of the inner “while” loop followed by the exchange of u and v transforms this density to V_{2}[f](x). However, by considering the first iteration of this loop (followed by the exchange if the loop terminates) we see that the transformed density is given by

V_{2}Ũ_{2}[f](x) + U_{2}[f](x),

where the first term arises if u < v without an exchange, and the second arises if an exchange occurs.

Consequence of Lemma 1

The following Theorem gives a simple relationship between B_{s}, V_{s} and U_{s}. The proof is immediate from Lemma 1 and the definitions of the operators.

Theorem 1

(V_{s} − I)U_{s} = V_{s}(B_{s} − I) .

Fixed Points

It follows immediately from Theorem 1 that, if

g = U_{2}f,

then

(V_{2} − I)g = V_{2}(B_{2} − I)f.

Thus, if f is a fixed point of the operator B_{2}, then g is a fixed point of the operator V_{2}. From a recent result of Vallée we know that V_{2}, acting on a certain Hardy space H^{2}(D), has a unique positive dominant simple eigenvalue 1, so g must be (a constant multiple of) the corresponding eigenfunction (provided g ∈ H^{2}(D)). Also, from the definitions of B_{2} and U_{2}, we have

λ = f(1) = 2g(1)

which is useful for proving the consistency of two of the expressions for K given below.

Some Recent Results of Vallée

Using her operator V_{s}, Vallée recently proved that

K = \frac{2 ln 2}{π^{2}g(1)} ∑_{{a odd, \atop a>0}} 2^{−⌊lg a⌋}G (\frac{1}{a})

where g is a nonzero fixed point of V_{2} (i.e. g = V_{2}g ≠ 0) and G(x) = ∫_{0}^{x} g(t) dt . This is the only expression for K which has been rigorously proved.

Because V_{s} can be proved to have nice spectral properties, the existence and uniqueness (up to scaling) of g can be proved rigorously.

Warning: Knuth uses G(x) for our F̃(x) = 1 − F (x) ! Unfortunately Knuth and Vallée use incompatible notation. We have followed Vallée, with the exception that we do not assume the normalisation G(1) = 1.

A Conjecture of Vallée

Let λ = f(1), where f is the limiting probability density (conjectured to exist) as above. Vallée (see Knuth, third edition, §4.5.2(61)) conjectured that

\frac{λ}{b} = \frac{2 ln 2}{π^{2}} ,

or equivalently that

K = \frac{4 ln 2}{π^{2}λ} .  (5)

Vallée proved the conjecture under the assumption that the operator B_{s} satisfies a certain spectral condition.

Numerical Results

Using an improvement of the “discretization method” of my 1976 paper, and the MP package with the equivalent of more than 50 decimal places (50D) working precision, we computed the limiting probability density f, then K, λ = f(1), and Kλ. The results were

K = 0.7059712461 0191639152 9314135852 8817666677
λ = 0.3979226811 8831664407 6707161142 6549823098
Kλ = 0.2809219710 9073150563 5754397987 9880385315

These are believed to be correctly rounded values.

Vallée’s conjecture (5) is that

Kλ = 4 ln 2/π^{2} .

The computed value of Kλ agrees with 4 ln 2/π^{2} to 40 decimals!

Correcting an Error

In my 1976 paper I claimed that, for all n ≥ 0 and x ∈ (0, 1],

F_{n}(x) = α_{n}(x) lg(x) + β_{n}(x) ,  (6)

where α_{n}(x) and β_{n}(x) are analytic and regular in the disk |x| < 1. However, this is incorrect, even in the case n = 1.

The error appeared to go unnoticed until 1997, when Don Knuth was revising Volume 2 in preparation for publication of the third edition. Knuth computed the constant K using recurrences for the analytic functions α_{n}(x) and β_{n}(x), and I computed K directly using the defining integral and recurrences for F_{n}(x). Our computations disagreed in the 14th decimal place ! Knuth found

K = 0.70597 12461 019\underline{45 99986} ⋯

but I found

K = 0.70597 12461 019\underline{16 39152} ⋯

Some Detective Work

After a flurry of emails we tracked down the error. It was found independently, and at the same time (within 24 hours), by Flajolet and Vallée.

The source of the error is illustrated by Lemma 3.1 of my 1976 paper, which is wrong (and corrected in the solution to ex. 4.5.2.29 of Knuth, third edition).

The Mellin transform of a function g(x) is defined by

g^{∗}(s) = ∫_{0}^{∞} g(x)x^{s−1}dx .

It is easy to see that, if

f(x) = ∑_{k≥1} 2^{−k}g(2^{k}x) ,

then the Mellin transform of f is

f^{∗}(s) = ∑_{k≥1} 2^{−k(s+1)}g^{∗}(s) = \frac{g^{∗}(s)}{2^{s+1} − 1} .

Mellin Inversion

Under suitable conditions we can apply the Mellin inversion formula to obtain

f(x) = \frac{1}{2πi} ∫_{c−i∞}^{c+i∞} f^{∗}(s)x^{−s}ds .

Applying these results to g(x) = 1/(1 + x), whose Mellin transform is g^{∗}(s) = π/ sin πs when 0 < Rs < 1, we find

f(x) = ∑_{k≥1} \frac{2^{−k}}{1 + 2^{k}x}

as a sum of residues of

(\frac{π}{sin πs}) \frac{x^{−s}}{2^{s+1} − 1}

for Rs ≤ 0. This gives

f(x) = 1+x lg x+ \frac{x}{2} +xP (lg x)− \frac{2}{1}x^{2}+ \frac{4}{3}x^{3}−⋯,

where

P (t) = \frac{2π}{ln 2} ∑_{n=1}^{∞} \frac{sin 2nπt}{sinh(2nπ^{2}/ ln 2)} .

The “Wobbles” Caused by P (t)

P (t) is a very small periodic function:

|P (t)| < 7.8 × 10^{−12}

for real t. In Lemma 3.1 of my 1976 paper, the term xP (lg x) was omitted.

Essentially, the poles off the real axis at

s = −1 ± \frac{2πin}{ln 2} , n = 1, 2, …

were ignored^{3}.

Because the residues at the non-real poles are tiny (thanks to the sinh term in the denominator) numerical computations performed using single-precision floating-point arithmetic did not reveal the error.

3In fact, the incorrect result was obtained without using Mellin transforms. If I had used them I probably would have obtained the correct result!

The Moral of this Tale

Computing results to (unreasonably) high precision by more than one method can be useful.

Small discrepancies should be explained, not ignored !

Conclusion and Open Problems

Since Vallée’s recent work, analysis of the average behaviour of the binary Euclidean algorithm has a rigorous foundation. However, some interesting open questions remain. For example, does the binary Euclidean operator B_{2} have a unique positive dominant simple eigenvalue 1? Vallée has proved the corresponding result for her operator V_{2}. In order to estimate the speed of convergence of f_{n} to f (assuming f exists), we need more information on the spectrum of B_{2}. What can be proved ? Preliminary numerical results indicate that the sub-dominant eigenvalue(s) are a complex conjugate pair:

λ_{2} = \overline{λ}_{3} = 0.1735 ± 0.0884i ,

with |λ_{2}| = |λ_{3}| = 0.1948 to 4D. Vallée has proved related results for some other algorithms (variants of the Euclidean algorithm, algorithms for computing the Jacobi symbol), but many analogous questions remain open.

Acknowledgements

Thanks to

• Don Knuth for encouraging me to correct and extend my 1976 results for the third edition of Seminumerical Algorithms,

• Brigitte Vallée for discussions and correspondence regarding her unpublished results and conjectures,

• Philippe Flajolet for his excellent notes on Mellin transforms.

Further details of the new results presented in this lecture will appear in a Technical Report.

References

[1] Richard P. Brent, Analysis of the Binary Euclidean Algorithm, New Directions and Recent Results in Algorithms and Complexity, (J. F. Traub, editor), Academic Press, New York, 1976, 321–355.

[2] C. A. R. Hoare, Quicksort, Comp. J. 5, 1 (1962), 10–15.

[3] Donald E. Knuth, The Art of Computer Programming, Volume 2: Seminumerical Algorithms (third edition). Addison-Wesley, Menlo Park, 1997.

[4] Brigitte Vallée, The complete analysis of the Binary Euclidean Algorithm, Proc. ANTS’98, Lecture Notes in Computer Science 1423, Springer-Verlag, 1998, 77–94.

[5] Brigitte Vallée, Dynamics of the binary Euclidean algorithm: functional analysis and operators, Algorithmica 22 (1998), 660–685.