6  Two-dimensional analysis

As already mentioned, averaging the number of divisions in the Euclidean algorithm over both inputs allows us to prove much more than we can in the one-dimensional analysis. There are various equivalent ways to formulate results: let

\[ \begin{aligned} U & = \{1,\ldots,N\} \times \{1,\ldots,N\} \\ U_{1} & = \{(a,b) \in U : \gcd(a,b) = 1\} \\ \Omega & = \{(a,b) \in U : a > b\} \\ \Omega_{1} & = \{(a,b) \in U_{1}: a > b\} \end{aligned} \]

then define random variables

\[ Y : U \to \mathbb{N}, \quad Y_1 : U_1 \to \mathbb{N}, \quad X : \Omega \to \mathbb{N}, \quad X_1 : \Omega_1 \to \mathbb{N} \]

by letting each be equal to \(T\) restricted to its sample space (\(Y(a,b) = T(a,b)\) for each \((a,b) \in U\), and so on), with \((a,b)\) being chosen uniformly at random from the sample space.

6.1 Average-case analysis

Thus, for instance,

\[ \mathbb{E}[Y] = \frac{1}{N^2} \sum_{1 \le a,b \le N} T(a,b) \tag{6.1}\]

is the average number of divisions performed by Euclid’s algorithm, over all pairs \((a,b)\) satisfying \(1 \le a,b \le N\). (We’ll discuss the restriction \(a > b\) shortly.) Knuth (Knuth 1997) gave a heuristic suggesting \(\mathbb{E}[Y] = \lambda\log N + O(1)\) (\(\lambda\) as in (4.4)), and numerical evidence suggesting the \(O(1)\) is actually around \(0.06\) plus a small error. This was first proved by Norton (Norton 1990) (whose proof hinged on Theorem 5.1 of Porter (Porter 1975)).

Theorem 6.1 (Norton, 1990) For all \(N \ge 2\) and any \(\epsilon > 0\),

\[ \mathbb{E}[Y] = \lambda \log N + \nu + O_{\epsilon}\left(N^{\epsilon - 1/6}\right), \tag{6.2}\]

with \(\lambda\) as in (4.4) and \(\nu\) as in (4.6).

As one might expect, if Conjecture 5.1 is true, the \(O\)-term in Norton’s theorem can be replaced by \(O_{\epsilon}\left(N^{\epsilon - 1/2}\right)\), and as we will see, the numerical evidence seems to support this.

By Norton’s result and the next proposition, the following estimates all hold and are in fact equivalent: for \(N \ge 2\) and \(\epsilon > 0\),

\[ \begin{aligned} \mathbb{E}[Y] & = \lambda \log N + \nu + O_{\epsilon}\left(N^{\epsilon - 1/6}\right) \nonumber \\ \mathbb{E}[Y_1] & = \lambda \log N + \nu_1 + O_{\epsilon}\left(N^{\epsilon - 1/6}\right) \\ \mathbb{E}[X] & = \lambda \log N + \nu - {\textstyle \frac{1}{2}} + O_{\epsilon}\left(N^{\epsilon - 1/6}\right) \\ \mathbb{E}[X_1] & = \lambda \log N + \nu_1 - {\textstyle \frac{1}{2}} + O_{\epsilon}\left(N^{\epsilon - 1/6}\right) \end{aligned} \tag{6.3}\]

with \(\lambda\), \(\nu\), and \(\nu_1\) as in (4.4) and (4.6), respectively. Again, we believe the \(O\)-term in all four estimates can be replaced by \(O_{\epsilon}(N^{\epsilon - 1/2})\) (and if such an \(O\)-term holds for one estimate, it holds in all four).

Proposition 6.1 (a) For \(N \ge 2\) we have

\[ \begin{aligned} \mathbb{E}[Y] = \mathbb{E}[X] + {\textstyle \frac{1}{2}} + O\left(N^{-1}\log N\right) \quad \text{and} \quad \mathbb{E}[Y_1] = \mathbb{E}[X_1] + {\textstyle \frac{1}{2}} + O\left(N^{-1}\log N\right). \end{aligned} \tag{6.4}\]

(b) The following statements are equivalent. (i) There exist constants \(\lambda\) and \(\nu_1\) such that, for \(N \ge 2\) and \(\epsilon > 0\),

\[ \mathbb{E}[Y_1] = \lambda \log N + \nu_1 + O_{\epsilon}\left(N^{\epsilon - 1/6}\right). \]

(ii) There exist constants \(\lambda\) and \(\nu_1\) such that, for \(N \ge 2\) and \(\epsilon > 0\),

\[ \mathbb{E}[Y] = \lambda \log N + \nu_1 + \lambda\frac{\zeta'(2)}{\zeta(2)}+ O_{\epsilon}\left(N^{\epsilon - 1/6}\right). \]

Proof. (a) For the second equation in (6.4), note that

\[ \begin{aligned} \#U_1 \mathbb{E}[Y_1] & = \sum_{(a,b) \in U_1} T(a,b) \\ & = \sum_{\substack{(a,b) \in U_1 \\ a > b}} T(a,b) + \sum_{\substack{(a,b) \in U_1 \\ b > a}} T(a,b) + \sum_{\substack{(a,b) \in U_1 \\ a = b}} T(a,b) \\ & = \sum_{\substack{(a,b) \in U_1 \\ a > b}} T(a,b) + \sum_{\substack{(a,b) \in U_1 \\ a > b}} T(b,a) + T(1,1) \\ & = \sum_{\substack{(a,b) \in U_1 \\ a > b}} T(a,b) + \sum_{\substack{(a,b) \in U_1 \\ a > b}} \left[T(a,b) + 1\right] + 1 \\ & = 2\sum_{\substack{(a,b) \in U_1 \\ a > b}} T(a,b) + \sum_{\substack{(a,b) \in U_1 \\ a > b}} 1 + 1 \\ & = 2\#\Omega_1 \mathbb{E}[X_1] + \#\Omega_1 + 1. \end{aligned} \]

Since \(\#U_1 = 2\#\Omega_1 + N\) (if this isn’t clear, replace \(T(a,b)\) by \(1\) in the above equations), this gives

\[ \#U_1\mathbb{E}[Y_1] = \#U_1 \mathbb{E}[X_1] - N\mathbb{E}[X_1] + {\textstyle \frac{1}{2}} \left(\#U_1 - N\right) + 1, \]

which upon dividing by \(\#U_1\) becomes

\[ \mathbb{E}[Y_1] = \mathbb{E}[X_1] + {\textstyle \frac{1}{2}} - \frac{N}{\# U_1}\left(\mathbb{E}[X_1] + {\textstyle \frac{1}{2} + \frac{1}{N}}\right). \]

All of this is valid for \(N \ge 1\). Proposition 3.3(c) implies that \(\mathbb{E}[X_1] \ll \log N\) for \(N \ge 2\). Also, a classical result of Mertens is that for \(N \ge 2\),

\[ \sum_{\substack{1 \le a,b \le N \\ \gcd(a,b) = 1}} 1 = \frac{N^2}{\zeta(2)} + O\left(N \log N\right). \tag{6.5}\]

We’ll use this in the proof of part (b), but at this point we only need that \(\#U_1 \gg N^2\). Combining gives the second equation in (6.4). The verification of the first equation in (6.4) is almost identical.

(b) We show that (i) implies (ii). Since \(\gcd(c,d) = 1\) if and only if \(\gcd(gc,gd) = g\), and since \(T(gc,gd) = T(c,d)\) (Proposition 3.1(c)),

\[ \begin{aligned} N^2 \mathbb{E}[Y] & = \sum_{1 \le a,b \le N} T(a,b) \\ & = \sum_{g = 1}^{N} \sum_{\substack{1 \le a,b \le N \\ \gcd(a,b) = g}} T(a,b) \\ & = \sum_{g = 1}^{N} \sum_{\substack{1 \le c,d \le N/g \\ \gcd(c,d) = 1}} T(gc,gd) \\ & = \sum_{g = 1}^{N} \sum_{\substack{1 \le c,d \le N/g \\ \gcd(c,d) = 1}} T(c,d) \\ & = \sum_{1 \le g \le N/2} \sum_{\substack{1 \le c,d \le N/g \\ \gcd(c,d) = 1}} T(c,d) + \sum_{N/2 < g \le N} T(1,1) \\ & = \sum_{1 \le g \le N/2} \sum_{\substack{1 \le c,d \le N/g \\ \gcd(c,d) = 1}} T(c,d) + O(N). \end{aligned} \tag{6.6}\]

We now assume (i), which is tantamount to assuming that, for \(M \ge 2\) and \(\epsilon > 0\),

\[ \sum_{\substack{1 \le c,d \le M \\ \gcd(c,d) = 1}} T(c,d) = \sum_{\substack{1 \le c,d \le M \\ \gcd(c,d) = 1}} \left(\lambda \log M + \nu_1 + O_{\epsilon}\left(M^{\epsilon - 1/6}\right)\right), \]

which, in view of (6.5), is equivalent to

\[ \sum_{\substack{1 \le c,d \le M \\ \gcd(c,d) = 1}} T(c,d) = \frac{M^2}{\zeta(2)}\left(\lambda \log M + \nu_1\right) + O_{\epsilon}\left(M^{2 + \epsilon - 1/6}\right). \]

Putting this into (6.6), we obtain

\[ \begin{aligned} \mathbb{E}[Y] & = \sum_{1 \le g \le N/2} \left[\left(\frac{(1/g)^2}{\zeta(2)}\right)\left(\lambda \log(N/g) + \nu_1\right) + O_{\epsilon}\left(N^{\epsilon - 1/6}/g^{2 + \epsilon - 1/6}\right)\right] + O\left(N^{-1}\right)\\ & = \lambda \log N \left(\frac{1}{\zeta(2)} \sum_{1 \le g \le N/2} \frac{1}{g^2}\right) + \lambda\left(\frac{1}{\zeta(2)} \sum_{1 \le g \le N/2} \frac{-\log g}{g^2}\right) + \nu_1 \left(\frac{1}{\zeta(2)} \sum_{1 \le g \le N/2} \frac{1}{g^2}\right) + O_{\epsilon}\left(N^{\epsilon - 1/6} \right). \end{aligned} \tag{6.7}\]

where we have used the fact that

\[\sum_{g = 1}^{\infty} \frac{1}{g^{2 + \epsilon - 1/6}} \ll 1\]

in the \(O\)-term. Finally,

\[ \sum_{1 \le g \le N/2} \frac{1}{g^2} = \sum_{g = 1}^{\infty} \frac{1}{g^2} + O\left(N^{-1}\right) = \zeta(2) + O\left(N^{-1}\right) \] and \[ \sum_{1 \le g \le N/2} \frac{-\log g}{g^2} = \sum_{g = 1}^{\infty} \frac{\log g}{g^2} + O\left(N^{-1}\log N\right) = \zeta'(2) + O\left(N^{-1}\log N\right), \]

and making these substitutions in (6.7) yields

\[ \mathbb{E}[Y] = \lambda\log N + \lambda \frac{\zeta'(2)}{\zeta(2)} + \nu_1 + O_{\epsilon}\left(N^{\epsilon - 1/6} \right). \]

Hence (i) implies (ii). The converse can be obtained similarly: we begin by noting that (\(\mu\) denoting the Möbius function)

\[ \sum_{\substack{1 \le a,b \le N \\ \gcd(a,b) = 1}} T(a,b) = \sum_{1 \le a,b \le N} T(a,b) \sum_{\substack{\delta \mid a \\ \delta \mid b}} \mu(\delta) = \sum_{1 \le \delta \le N} \mu(\delta) \sum_{1 \le c,d \le N/\delta} T(\delta c,\delta d) = \sum_{1 \le \delta \le N} \mu(\delta) \sum_{1 \le c,d \le N/\delta} T(c,d). \]

We leave the remainder of the proof as an exercise for the reader.

Naturally, the next questions are related to variance, and, ultimately, distribution.

6.2 Variance and distribution

Dixon (Dixon 1970) showed that \(|T(a,b) - \lambda\log a| < (\log a)^{1/2 + \epsilon}\) (with \(\lambda\) as in (4.4)) for all but at most \(N^2\exp\left(-c(\log N)^{\epsilon/2}\right)\) (\(c\) a constant) of the pairs \((a,b)\) in the range \(1 \le b < a \le N\). This may hint at a variance, but it was Hensley (Hensley 1994) who, in a tour de force of functional analysis, first attained a sharp estimate for the variance of \(X_1\) (and \(X\)), and also established that the distribution of \(X\) (and of \(X_1\)) is asymptotically normal.

Theorem 6.2 (Hensley, 1994) There exist constants \(\lambda\) and \(\eta\) (as in (4.4) and (4.7)), such that the following holds for all sufficiently large \(N\). Let \(\mu = \lambda \log N\), and \(\sigma = \sqrt{\eta \log N}\). Uniformly for \(|s - \mu| \le \frac{1}{2}(\log N)^{1/2}(\log\log N)^2\),

\[ \mathbb{P}(X_1 = s) = \left(1 + O\left((\log N)^{-1/24}\right)\right) \frac{\exp\left(-\frac{1}{2}\left(\frac{s - \mu}{\sigma}\right)^2\right)}{\sigma \sqrt{2\pi}}. \tag{6.8}\]

The same is true with \(X\) in place of \(X_1\).

A unifying framework for analysing a general class of Euclidean (and other) algorithms has been developed, principally, by Vallée. The methods rely on properties of transfer operators suitably adapted from dynamical systems theory. The following is a special case of a very deep result of Baladi and Vallée (Baladi and Vallée 2005).

Theorem 6.3 (Baladi and Vallée, 2005) There exist constants \(\lambda\), \(\nu_1\), \(\eta\) (as in (4.4), (4.6), and (4.7)), and \(\kappa_1\) such that, for all \(N \ge 2\) and any \(\epsilon > 0\),

\[ \mathbb{E}[X_1] = \lambda\log N + \nu_1 + O_{\epsilon}\left(N^{\epsilon - 1/6}\right) \tag{6.9}\]

and

\[ \mathrm{Var}(X_1) = \eta\log N + \kappa_1 + O\left(N^{-c}\right) \tag{6.10}\]

(for some positive constant \(c\)). Moreover, \(X_1\) is asymptotically normal: for all \(N \ge 2\) and all \(z\),

\[ \mathbb{P}\left(\frac{X_1 - \lambda\log N}{\sqrt{\eta\log N}} \le z\right) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^z e^{-t^2/2} dt + O\left(\left(\log N\right)^{-1/2}\right). \tag{6.11}\]

The same is true of \(X\), with suitable constants \(\nu\) and \(\kappa\) in place of \(\nu_1\) and \(\kappa_1\) (the constants \(\lambda\) and \(\eta\) remain the same).

As we mentioned following Theorem 6.1 of Norton, if Conjecture 5.1 is true, the \(O\)-term in (6.9) can be replaced by \(O_{\epsilon}\left(N^{\epsilon - 1/2}\right)\). We also believe that the \(O\)-term in (6.10) can be replaced by \(O_{\epsilon}\left(N^{\epsilon - 1/2}\right)\).

In the variance (6.10), Hensley’s constant \(\eta\) (see (4.7)) does not admit a simple closed form (it can be described in terms of the spectrum of a certain transfer operator), but was shown by Lhote (Lhote 2005) to be polynomial-time computable, and the first seven digits (\(\eta = 0.5160524\ldots\)) were given by Lhote op. cit. The constant \(\kappa_1\) is “even more obscure” (in the words of Baladi and Vallée (Baladi and Vallée 2005)), and to our knowledge has not been investigated numerically. What we can say is that, assuming our back-of-the-envelope calculation is correct (we invite the reader to verify this along the lines of Proposition 6.1(b)), \(\mathrm{Var}[X] - \mathrm{Var}[X_1]\) is (up to an \(O\)-term)

\[ \kappa - \kappa_1 = \eta\frac{\zeta'(2)}{\zeta(2)} + \lambda^2\left(\frac{\zeta''(2)}{\zeta(2)} - \left(\frac{\zeta'(2)}{\zeta(2)}\right)^2\right) = 0.3340\ldots \]

(see (4.8)). Our numerical data leads us to “guesstimate” the numerical value of \(\kappa\) to be somewhere around \(-0.1\).

Naturally, \(Y\) and \(Y_1\) are asymptotically normal as well, and analogously to (6.4) (another exercise for the reader), we have

\[ \begin{aligned} \mathrm{Var}[Y] = \mathrm{Var}[X] + {\textstyle \frac{1}{4}} + O\left(N^{-1}(\log N)^2\right) \quad \text{and} \quad \mathrm{Var}[Y_1] = \mathrm{Var}[X_1] + {\textstyle \frac{1}{4}} + O\left(N^{-1}(\log N)^2\right). \end{aligned} \tag{6.12}\]

Baladi, Viviane, and Brigitte Vallée. 2005. “Euclidean Algorithms Are Gaussian.” Journal of Number Theory 110 (2): 331–86. https://doi.org/10.1016/j.jnt.2004.08.008.

Dixon, John D. 1970. “The Number of Steps in the Euclidean Algorithm.” Journal of Number Theory 2 (4): 414–22. https://doi.org/10.1016/0022-314X(70)90044-2.

Hensley, Doug. 1994. “The Number of Steps in the Euclidean Algorithm.” Journal of Number Theory 49 (2): 142–82. https://doi.org/10.1006/jnth.1994.1088.

Knuth, Donald E. 1997. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. 3rd ed. Addison-Wesley.

Lhote, Loïck. 2005. “Efficient Computation of a Class of Continued Fraction Constants.” In Algorithms Seminar 2002–2004, edited by Frédéric Chyzak.

Norton, G. H. 1990. “On the Asymptotic Analysis of the Euclidean Algorithm.” Journal of Symbolic Computation 10 (1): 53–58. https://doi.org/10.1016/S0747-7171(08)80036-3.

Porter, J. W. 1975. “On a Theorem of Heilbronn.” Mathematika 22 (1): 20–28. https://doi.org/10.1112/S0025579300004459.