Complete Analysis of the Binary GCD Algorithm

[a text written by Cyril Banderier, based on the talk given by Brigitte Valle 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
Dn ~
12 log 2
p2
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
2
Val2(v-u)
 
,v



gcd(u,v)=gcd(v,u)
where Val2(n) is the greatest integer b such 2b divides n, i.e., the dyadic valuation of n. The corresponding binary GCD algorithm is as follows:
Example 1   If the input is (u,v):=(7,61) then b:=Val2(61-7)=1. Thus v:=54/21=27, and the algorithm continues because u<v. Now b:=Val2(27-7)=2. Thus v:=20/22=5. Now the algorithm restarts with (u,v):=(5,7). It leads to v:=(7-5)/21=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:
7
61
=
1
3+
23
1+
21
1+22
.

In general, for each ``inner while loop'', one has
xi=
1
ai+2
ki
 
xi+1
where xi:=u/v (with (u,v) as in the beginning of the loop), xi+1:=u/v (with (u,v) as after the exchange), where ai:=1+2b1+2b1+b2+...+2b1+...+bl-1 and ki:=b1+...+bl-1+bl (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
u
v
=
1
a1+
2
k1
 
...+
2
kr-1
 
ar+2
kr
 
.
Three interesting parameters are: Their average values on the set Wn={(u,v) odd, 0<u v n} are respectively noted En, Pn and Sn. Note that En is also the average number of exchanges in the algorithm, and that Pn 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:
Vs[f](x):=
 
k 1
 
a odd
1 a 2k
1
(a+2k x)s
f


1
a +2k x



.
The average En is easily expressed in term of Vs, with the help of the following definitions:
F(s):=(Id-Vs)-1[Id](1),    G(s):=(Id-Vs)-2 Vs[Id ](1),     
~
z
 
(s):=
 
k odd
1
ks
=


1-
1
2s



z(s).

Proposition 1   En 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,
Vsl[f](1)=
1
~
z
 
(s)
 
(u,v)W[l]
1
vs
f


u
v



.
Summing over all the possible heights (l 0) yields:
(Id -w Vs)-1[f](1)=
 
l 0
wl Vsl[f](1) =
1
~
z
 
(s)
 
(u,v)W[l]
1
vs
f


u
v



.
Differentiating with respect to w, and then choosing f=1 and w=1 yields
En=
1
| Wn |
 
l 0
l | Wn[l]| =
 
l 0
l
 
k n
vk[l]
 
l 0
 
k n
vk[l]
.
The proof is completed by observing that
F(s)=
1
~
z
 
(s)
 
k 1
1
vs
 
l 0
vk[l],     G(s)=
1
~
z
 
(s)
 
k 1
1
vs
 
l 0
l vk[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
  1. F is analytic on the line (s)=s except at s=s;
  2. 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
ak =
A(s)
s G(g+1)
n
s
 
log
g
 
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-Vs)-1[Id](1)=1+
1
2
~
z
 
(s)
 
v odd
v-1
vs
=
1
2






~
z
 
(s-1)
~
z
 
(s)
+1





.
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 Valle'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 Vs 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 H2(D), where D is an open disk containing ]0,1]. This leads to the fact that for s>3/2, Vs has a unique positive dominant eigenvalue, equal to 1 when s=2. In addition Vs has a spectral radius <1 on (s) 2, s 2. Thus (Id-Vs)-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 Vs=Ps+Ns (Ps is the projection of Vs 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 Wn is
En ~
2
p2 f2(1)
log n,
where f2 is the fixed point of the operator V2 that is normalised by 01 f2(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 Valle [7]. This leads to the following two results.
Theorem 3   The average number of total iterations is
Pn ~ A log n    with   A:=
2
p2 f2(1)
 
a odd
1
2
ka
 
F2


1
a



where f2 is defined as above, F2(x):=0x f2(t) dt, F2(1)=1 (where ka is the integer part of log2 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
Sn ~
2
p2 f2(1)




2
 
a odd
1
2
ka
 
F2


1
a







log n.

4   All Roads Lead to Rome

In Brent's paper [1], one can find a different approach which suggests that
Pn ~
1
M
log n    where   M=log 2-
1
2

1


0
log (1-x) g2(x) dx
and where g2 is the fixed point (and normalised as f2) of
B2[f](x):=
 
b 1



1
1+2b x



2 f


1
1+2b x



+
 
b 1



1
x+2b



2 f


x
x+2b



.
Unfortunately, this approach is based on a heuristic hypothesis (exercise 36, section 4.5.2, rated HM49 by Knuth in [5]). Brigitte Valle explored this approach with a Brent operator Bs, without heuristic arguments but providing a spectral conjecture holds, this leads to the following result:
Pn ~ B log n     where   B:=
4
p2 g2(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 Valle 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). -- Gnralisation du thorme de Ikehara. Annales Scientifiques de l'cole Normale Suprieure, 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]
Leme (Charlie) and Valle (Brigitte). -- Analyse des algorithmes du symbole de Jacobi. GREYC, 1998.

[7]
Valle (Brigitte). -- The complete analysis of the binary Euclidean algorithm. In Proceedings ANTS'98. -- 1998.