Complete Analysis of the Binary GCD Algorithm

Complete Analysis of the Binary GCD Algorithm

[a text written by Cyril Banderier, based on the talk given by Brigitte Vallée at the Algoriths Seminar, April 27, 1998]

A properly typeset version of this document is available in postscript and in pdf.

1 Introduction

The analysis of the classical Euclidean algorithm has been performed by Heilbronn [4] and Dixon [3], using different approaches. For a random pair of rational numbers, the average number of divisions is

D_n ~

12 log 2

p²

log n.

Here, we will analyse the binary Euclidean algorithm, which uses only subtractions and right binary shifts. This ``binary GCD algorithm'' takes as input a pair of odd integers (u,v) from the set W={(u,v) odd, 0<u £ v}. Then the GCD is recursively defined by

ì
ï
í
ï
î

gcd(u,v)=gcd

æ
ç
ç
ç
è

v-u

Val₂(v-u)

ö
÷
÷
÷
ø

gcd(u,v)=gcd(v,u)

where Val₂(n) is the greatest integer b such 2^b divides n, i.e., the dyadic valuation of n. The corresponding binary GCD algorithm is as follows:

while u¹ v do
- while u<v do
  - b:=Val₂(v-u);
  - v:=(v-u)/2^b;
- end;
- exchange u and v;
end;
return u.

Example 1 If the input is (u,v):=(7,61) then b:=Val₂(61-7)=1. Thus v:=54/2¹=27, and the algorithm continues because u<v. Now b:=Val₂(27-7)=2. Thus v:=20/2²=5. Now the algorithm restarts with (u,v):=(5,7). It leads to v:=(7-5)/2¹=1 and therefore one restarts with (u,v):=(1,5) which leads to v=1=u so the algorithm stops and returns u, namely 1 (as expected since 7 and 61 are coprime). One can write:

2³

2¹

1+2²

In general, for each ``inner while loop'', one has

x_i=

a_i+2

k_i

x_i+1

where x_i:=u/v (with (u,v) as in the beginning of the loop), x_i+1:=u/v (with (u,v) as after the exchange), where a_i:=1+2^b_¹+2^b_¹^+b_²+...+2^b_¹^+...+b_^l-1 and k_i:=b₁+...+b_l-1+b_l (the sum of all the b's obtained during the i-th inner while loop). The algorithm thus produces the following binary continued fraction expansion

a₁+

k₁

...+

k_r-1

a_r+2

k_r

Three interesting parameters are:

r, the depth of the continued fraction or equivalently the number of outer loops performed;
å_i=1^r n(a_i), the number of subtractions (where n(w) is the number of 1's in the binary expansion of the integer w);
å_i=1^r k_i, number of rights shifts performed or equivalently inner loop executions.

Their average values on the set W_n={(u,v) odd, 0<u £ v £ n} are respectively noted E_n, P_n and S_n. Note that E_n is also the average number of exchanges in the algorithm, and that P_n is the average number of operations that are necessary to obtain the expansion.

2 A Ruelle Operator for a Tauberian Theorem

In order to establish that these three parameters have averages that are asymptotic to log n, we introduce the following Ruelle operator:

V_s[f](x):=

k ³ 1

a odd

1£ a £ 2^k

(a+2^k x)^s

æ
ç
ç
è

a +2^k x

ö
÷
÷
ø

The average E_n is easily expressed in term of V_s, with the help of the following definitions:

F(s):=(Id-V_s)^-1[Id](1), G(s):=(Id-V_s)^-2 ° V_s[Id

](1),

(s):=

k odd

k^s

æ
ç
ç
è

2^s

ö
÷
÷
ø

z(s).

Proposition 1 E_n is a ratio of partial sums of the two Dirichlet series z^~(s) F(s) and z^~(s) G(s).

Proof. Let W^[l] be the subset of W for which the algorithm performs exactly l exchanges. Then,

V_s^l[f](1)=

(s)

(u,v)ÎW^[l]

v^s

æ
ç
ç
è

ö
÷
÷
ø

Summing over all the possible heights (l³ 0) yields:

(Id

-w V_s)^-1[f](1)=

l ³ 0

w^l V_s^l[f](1) =

(s)

(u,v)ÎW^[l]

v^s

æ
ç
ç
è

ö
÷
÷
ø

Differentiating with respect to w, and then choosing f=1 and w=1 yields

E_n=

| W_n |

l³ 0

l | W_n^[l]| =

l³ 0

k£ n

v_k^[l]

l³ 0

k£ n

v_k^[l]

The proof is completed by observing that

F(s)=

(s)

k ³ 1

v^s

l³ 0

v_k^[l], G(s)=

(s)

k ³ 1

v^s

l³ 0

l v_k^[l].

The key is now to prove that the following theorem may be used:

Theorem 1 [Tauberian theorem] If F(s) is a Dirichlet series with non-negative coefficients that is convergent for Â(s)>s>0 and if

F is analytic on the line Â(s)=s except at s=s;
F(s)=A(s)/(s-s)^g+1+C(s) where A,C are analytic at s (with A(s) ¹ 0);

then one has, as n ® ¥,

k £ n

a_k =

A(s)

s G(g+1)

log

n (1+e(n)),

where e(n) ® 0.

Proof. See Delange [2].

Lemma 1 The Tauberian theorem applies to F with s=2 and g=0.

Proof. Indeed

F(s):=(Id-V_s)^-1[Id](1)=1+

(s)

v odd

v-1

v^s

æ
ç
ç
ç
ç
ç
è

(s-1)

(s)

ö
÷
÷
÷
÷
÷
ø

The last member of the equality clearly satisfies the conditions of the Tauberian theorem, and the same holds for z^~F with s=2 and g=0.

Lemma 2 The Tauberian theorem applies to G with s=2 and g=1.

Proof. Here lies the complex part of Brigitte Vallée's proof. It is impossible to conclude as quickly as in lemma 1, indeed, this time we need to find an appropriate functional space on which V_s is a compact operator. A mixture of various functional analysis theorems (Fejer-Riesz' inequality, Gabriel's inequality, Krasnoselsky's theorem and other works by Shapiro and Grothendieck) show that it is the case on the Hardy space H²(D), where D is an open disk containing ]0,1]. This leads to the fact that for s>3/2, V_s has a unique positive dominant eigenvalue, equal to 1 when s=2. In addition V_s has a spectral radius <1 on Â(s)³ 2, s¹ 2. Thus (Id-V_s)^-1 is regular on the domain D and condition 1 of the Tauberian theorem is fulfilled. Condition 2 is proved by means of perturbation theory applied to V_s=P_s+N_s (P_s is the projection of V_s on the dominant eigensubspace), in a neighbourhood of s=2. See [7] for a detailed proof.

This implies the following fundamental result:

Theorem 2 The average number of exchanges of the binary Euclidean algorithm on W_n is

E_n ~

p² f₂(1)

log n,

where f₂ is the fixed point of the operator V₂ that is normalised by ò₀¹ f₂(t) dt =1.

3 The Other Two Parameters

In order to study the other two parameters (total number of subtractions, total number of shifts) one still uses the Tauberian theorem but with a more intricate Ruelle operator, see Vallée [7]. This leads to the following two results.

Theorem 3 The average number of total iterations is

P_n ~ A log n with A:=

p² f₂(1)

a odd

k_a

F₂

æ
ç
ç
è

ö
÷
÷
ø

where f₂ is defined as above, F₂(x):=ò₀^x f₂(t) dt, F₂(1)=1 (where k_a is the integer part of log₂ a).

Theorem 4 The average number of the sum of exponents of 2 used in the numerators of the binary continued fraction expansions, i.e., average total number of right shifts is

S_n ~

p² f₂(1)

æ
ç
ç
ç
è

a odd

k_a

F₂

æ
ç
ç
è

ö
÷
÷
ø

ö
÷
÷
÷
ø

log n.

4 All Roads Lead to Rome

In Brent's paper [1], one can find a different approach which suggests that

P_n ~

log n where M=log 2-

ó
õ

log (1-x) g₂(x) dx

and where g₂ is the fixed point (and normalised as f₂) of

B₂[f](x):=

b ³ 1

æ
ç
ç
è

1+2^b x

ö
÷
÷
ø

² f

æ
ç
ç
è

1+2^b x

ö
÷
÷
ø

b ³ 1

æ
ç
ç
è

x+2^b

ö
÷
÷
ø

² f

æ
ç
ç
è

x+2^b

ö
÷
÷
ø

Unfortunately, this approach is based on a heuristic hypothesis (exercise 36, section 4.5.2, rated HM49 by Knuth in [5]). Brigitte Vallée explored this approach with a Brent operator B_s, without heuristic arguments but providing a spectral conjecture holds, this leads to the following result:

P_n ~ B log n where B:=

p² g₂(1)

The miracle holds and, after numerical experiments, A=1/M=B=1.0185.... But nobody has proved these equalities. We can also note that a similar method was used by Brigitte Vallée and one of her students to analyse the Jacobi symbol algorithm [6]. Finally, the binary Euclidian algorithm is only a slight variation on one of the oldest known algorithms but there is still some unknown territories in its ``complete'' analysis!

References

[1]: Brent (Richard P.). -- Analysis of the binary Euclidean algorithm. In Algorithms and complexity, pp. 321--355. -- Academic Press, New York, 1976. Proceedings of a Symposium held at Carnegie-Mellon University, 1976.
[2]: Delange (Hubert). -- Généralisation du théorème de Ikehara. Annales Scientifiques de l'École Normale Supérieure, vol. 71, n°3, 1954, pp. 213--242.
[3]: Dixon (John D.). -- The number of steps in the Euclidean algorithm. Journal of Number Theory, vol. 2, 1970.
[4]: Heilbronn (H.). -- On the average length of a class of finite continued fractions. In Number Theory and Analysis (Papers in Honor of Edmund Landau), pp. 87--96. -- Plenum, New York, 1969.
[5]: Knuth (Donald E.). -- The Art of Computer Programming. -- Addison-Wesley, 1997, third edition, vol. 2.
[6]: Lemée (Charlie) and Vallée (Brigitte). -- Analyse des algorithmes du symbole de Jacobi. GREYC, 1998.
[7]: Vallée (Brigitte). -- The complete analysis of the binary Euclidean algorithm. In Proceedings ANTS'98. -- 1998.