# Statistics of Extended GCD

Consider a coprime pair of integers $a, b.$ As we all know ("Bezout's theorem") there is a pair of integers $c, d$ such that $ac + bd=1.$ Consider the smallest (in the sense of Euclidean norm) such pair $c_0, d_0$, and consider the ratio $\frac{\|(c_0, d_0)\|}{\|(a, b)\|}.$ The question is: what is the statistics of this ratio as $(a, b)$ ranges over all visible pairs in, for example, the square $1\leq a \leq N, 1 \leq b \leq N?$

Experiment shows the following amazing histogram:

EDIT by popular demand: the histogram is for an experiment for $N=1000.$ The $x$ axis is the ratio, the $y$ axis is the number of points in the bin. The total number of points is $1000000/\zeta(2),$ so there are $100$ bins each with around $6000$ points.

But no immediate proof leaps to mind.

-
Can you label your axes? I take it $x$-axis is the ratio, $y$ axis is number of pairs, $N$ is some moderately large number? –  Gerry Myerson Jan 22 '11 at 22:11
@Gerry: yes, you guessed right. $N$ is 1000 in this histogram. –  Igor Rivin Jan 22 '11 at 22:31
It looks as if you've got a bin-width of about 0.003, so about 166 bins and somewhat over 6000 cases per bin. It would be less distracting from the point of your question if you said all this explicitly. –  Michael Hardy Jan 22 '11 at 23:31
Perhaps rewrite it as $|ac-bd|=1$ with everything positive and examine for each $c,d$ which pairs $a,b$ it "works" for. $d/c$ should be a next to the last convergent to $a/b$ so... –  Aaron Meyerowitz Jan 22 '11 at 23:53
Concerning the importance of labeling axes, see xkcd.com/833 –  Gerry Myerson Jan 23 '11 at 5:44

I did a little experiment. Fix $a=29$, let $b=1,2,\dots,28$. So, you get 28 data points. Well, these points are already extremely regularly distributed. Taking just the first half, $1\le b\le14$, and rearranging the ratios in increasing order, they are (to three decimals) $$.034,.069,.103,.138,.172,.207,.242,.275,.310,.345,.379,.414,.448,.483$$ To three decimals, and modulo round-off errors, these are the numbers $1/29,2/29,\dots,14/29$, which is to say they are about as regularly distributed as possible. The ratios for $15\le b\le28$ are essentially the same numbers - in fact, the ratio for $(a,b)$ seems to be pretty nearly the ratio for $(b-a,b)$.

If what's happening for 29 happens in general, I think it would explain the original histogram.

EDIT: So I think I see what's going on. We're looking at the numbers $$\sqrt{c^2+d^2\over a^2+b^2}$$ But $b$ is very close to $-ac/d$ (since $ac+bd=1$), so these numbers are very close to $$\sqrt{c^2+d^2\over a^2+(ac/d)^2}$$ which simplifies to $|d|/a$. For fixed $a$, as $b$ runs through the units modulo $a$, so does $d$, since $bd\equiv1\pmod a$. So our ratios are as uniformly distributed as the fractions $|d|/a$, which is very.

-
The bijection (reversing continued fractions) in my remark explains this. –  David Feldman Jan 23 '11 at 5:35
"Explains" in what sense? Proves that for fixed $a$ you get pretty nearly exactly the numbers $1/a$, $2/a$, etc.? I don't see how reversing continued fractions gets information this specific. –  Gerry Myerson Jan 23 '11 at 5:51
Because denominators stay invariant on reversing the continued fraction (if you do it right). Look at your own example in that notation. –  David Feldman Jan 23 '11 at 5:57
@Gerry: this explains many things! A question: why are $|d|/a$ "very uniformly distributed"? –  Igor Rivin Jan 23 '11 at 15:20
@Gerry: presumably, one can write down an error term (for deviation from uniformity), which probably has much to do with @Bill's mixing rate for horocycle flow... –  Igor Rivin Jan 23 '11 at 15:22

Here's a more geometric formulation of your question:

On the torus $\mathbb R^2/\mathbb Z^2$, consider a long simple closed geodesic $\overline {(0,0)(a,b)}$. It cuts the torus into a thin cylinder; the cylinder is joined to itself by a twist by some angle to form the torus. What is the distribution of the angle of the twist?

From this perspective, perhaps your intuition tells you that the angles of twist should tend toward the uniform distribution, as the homotopy classes of geodesics are chosen uniformly with longer and longer lengths.

To get a rigorous argument, we can think about the problem from the opposite direction. Begin by starting from a long thin annulus of area 1, and ask what are the ways to glue it together to form a torus? You can glue it by any angle; however, the torus you get is not usually isometric to the square torus; but it is isometric to $\mathbb E^2$ modulo some lattice of area 1.

The space of lattices up to similarity in $\mathbb E^2$ together with a choice of positively oriented generators is the Teichmüller space for the torus, and can be identified with the hyperbolic plane, in the upper halfspace model: make the first vector go from $0$ to $1$ along the $x$-axis, and the second vector will be a point $z$ in upper half space. Change of generators acts by fractional linear transformations, and preserves the hyperbolic metric; the quotient is the space of isometry classes of Euclidean tori of area 1, the moduli space of the torus.

Twisting an annulus is a well-studied operation, the horocycle flow, on the moduli space. The action preserves a probably familiar tiling of the hyperbolic plane by ideal triangles whose vertices are at rational points on the bounding line, completed to make $\mathbb {RP}^1$, and whose edges connect pairs of slopes corresponding to your equation. The points in Teichmüller space that give square toruses are the midpoints of the edges. (Actually, the edges are infinitely long, so they don't have an obvious definition of midpoint, but the triangles have altitudes whose feet we can call the midpoints: these are the square toruses).

In any case: if you look at all lattice vectors with length between say $N$ and $2N$, and ask for the distribution of angles among these, this is equivalent to taking all points representing square lattices in a band in Teichmüller space between horocycles that appear in upper half space as a rectangle bounded by horizontal lines at height $1/N$ and $1/2N$ and vertical lines $x = \pm 1/2$; the question is the distribution of $x$-coordinates of points representing square toruses within that band. That this tends to the uniform distribution follows from ergodicity of the horocycle flow, a well-known fact whose history probably predates the sources I'm familiar with, so I won't try to give the attribution.

-
-
@Bill: This is certainly much of the answer, but a couple of questions: ergodicity implies that any reasonable function will have a limiting distribution, but the one I had picked is uniform. This implies that the ratio I defined is the $x$-coordinate of square toruses. This is not immediately obvious to me (but I have not completely absorbed what you said) A second point is that $1000$ is not a very large number, and yet the distribution is essentially perfectly uniform. So, the mixing must be very rapid, which is not an obvious fact. ($N=100$ gives a very similar histogram, by the way...) –  Igor Rivin Jan 23 '11 at 2:00
@Bill: By the way, a brief googling seems to show that the ergodicity of horocycle flow is a lot more recent than one might have thought ('60s-70s). –  Igor Rivin Jan 23 '11 at 2:01
This is a very nice answer. I like it a lot. –  Eric Naslund Jan 23 '11 at 2:50
I think ergodicity of the Horocycle flow goes back to Hedlund in the 1930s. Unique ergodicity (which might be needed here) was shown by Furstenberg in the 1970s (en.wikipedia.org/wiki/Ergodic_theory#Ergodic_flows_on_manifolds) –  Anthony Quas Jan 23 '11 at 3:11
@Igor: You can consider averages of the delta measure of square toruses smeared out into a small ball. Boundary issues work out since averages converge over horocycles, not just over the area, it works out. The horocycle flow mixes particularly quickly, as indicated by unique ergodicity. That's why the bins are so even. It's closely related to self-encoding sequences. @Quas: Thanks. I was thinking Hedlund was likely involved, but I suspect there are earlier versions not phrased as ergodicity, but perhaps in terms of elliptic functions or Poincare series or modular forms or whatnot. –  Bill Thurston Jan 23 '11 at 3:33

Roughly:

Suppose you have a fraction $a/b$ and you expand it as a simple continued fraction $1/c_1+1/c_2\cdots+1/c_{n-1}+1/c_n$. Now truncate the last convergent and collapse $1/c_1+1/c_2\cdots+1/c_{n-1}$ to get, say, $a'/b'$. Now consider the simple continued fraction expansion of $b'/b$. As I recall, this will (usually? perhaps I need a hypothesis to avoid degenerate cases?) equal the reverse of the continued fraction of $a/b$.

For complicated fractions the beginning and end of a continued fraction should be almost uncorrelated. Also for a random fraction the convergents have a known distribution (Gauss-Kuzmin) and reversing the continued fraction doesn't change the distribution of the convergents, so you get the uniform distribution back.

-
This is certainly very interesting. Any reference for the first paragraph? I don't think I have seen this... –  Igor Rivin Jan 23 '11 at 4:37
Igor, I discovered that as an undergraduate, always assumed that I'd rediscovered the wheel and never bothered to look it up so I have no reference. But I believe it follows now that I'm thinking about it afresh just from the looking at continued fractions from the matrix point of view (as products of matrices $\begin{array}{c} c_i & 1 \\ 1 & 0 \end{array}$) and then just taking the adjoint (which reverses products and switches present numerator and historical denominator). –  David Feldman Jan 23 '11 at 4:48
Just to be clear, there are always two ways to end a simple continued fraction, either ...$c$ or ...$c-1, 1$ with $c>1$. These have distinct reversals, one in the interval $[0,.5]$ and the other in the interval $[.5,1]$. So to get a "bijection" you have to pick an interval on the one hand and a normal form on the other. –  David Feldman Jan 23 '11 at 6:26
I remember now that my original interest in continued fraction reversal had to do with music theory! The ratio of two frequencies comes to one's ears as a real number, but its "meaning" in Western music requires interpreting that number as representing a nearby rational. Small denominators make the interpretation unambiguous, but more complex fractions ("intervals") require context in order to acquire a unique interpretation ("function"). Add a bass ($\alpha$ for frequences $a\alpha$ and $b\alpha$ does it, but $\alpha$ may sound so low as not to work in a practical way. –  David Feldman Jan 23 '11 at 22:53
... But juxtaposing two competing interpretations (a "modulation" via a "pivot chord" one has two $\alpha$'s and even in the wrong register, the ratio between the $\alpha$'s has the reverse continued fraction of the ratio between $a$ and $b$ and thus pins down the least audible information that ratio meant to carry. –  David Feldman Jan 23 '11 at 22:55

Here is an attempt to give a somewhat finer grained view of the distribution. The set of ratios $\sqrt{\frac{s^2+t^2}{a^2+b^2}} \subset(0,\frac{1}{2})$ are essentially the values in the first half of the Farey sequence $\lbrace \frac{p}{q} | \gcd(p,q)=1,\ 2p \le q<N\ \rbrace$. This has already been pointed out but I'll give a simple (if less nuanced) justification. Then I'll mention how that sequence is and is not smoothly distributed.

Instead of looking at all the relatively prime pairs $(a,b)$ with $1 \le a,b\le N$ I'll just consider those with $a<b$ since order is irrelevant for the question asked and $(a,b)=(1,1)$ is an extreme outlier. There are non-negative integers $s,t$ with $|as-bt|=1$ and just one such pair with $\sqrt{\frac{s^2+t^2}{a^2+b^2}}<\frac{1}{2}$. This ratio turns out to be very close to $\frac{t}{s}$. Then $(0,0),(t,s),(a,b)$ and $(t+a,s+b)$ are corners of a long thin parallelogram with area 1 and (thus) no integer points on its boundary or interior. Because the sides are very nearly parallel, the ratio $\sqrt{\frac{s^2+t^2}{a^2+b^2}}$ of their lengths is quite close to $\frac{t}{a}$ and even closer to $\frac{s}{b}$ (in fact they are convergents to the continued fraction for that irrational number). So that set of ratios is quite close to the lower half of the set of fractions

LATER For fixed $b$ the discrepency between $\sqrt{\frac{s^2+t^2}{a^2+b^2}}$ and $\frac{s}{b}$ is almost exactly $\frac{1}{b(2b^2+1)}$ for $\frac{1}{b}$ and increases to $\frac{2}{b(4b-1)}$ for $\frac{b-1}{b}$

Let $\mathcal{H}_N=\lbrace \frac {p}{q} |\frac{p}{q}\le \frac{1}{2} ,\gcd(p,q)=1,q \le N \rbrace$ The letter $\mathcal{H}$ is because this is half a Farey sequence. It is known that $P(N)=|\mathcal{H}_N|=\frac{3N^2}{2\pi^2}+O(N\log N)$. How evenly spaced are these? There are $P \approxeq \frac{0.15}{N^2}$ points in an interval of width $1/2$ so perfectly even spacing would put the kth point at $\frac{k}{2P}\approx\frac{3.3k}{N^2}$. However a fraction $\frac{p}{q}$ with $q$ small will be about $\frac{1}{qN}$ from the next nearest points. Hence the largest point other than $\frac{1}{2}$ is $\frac{1}{2}-\frac{1}{N}$ (replace N by N-1 in the even case) and the smallest point is $\frac{1}{N}$ which seems far from $\frac{3.3}{N^2}$ These empty zones force other points closer together, the first few points are only separated by about $\frac{1}{N^2}$. I can't resist an attempt to put in a picture of Ford Circles. A disk of radius $\frac{1}{q^2}$ is centered at $(\frac{p}{q},\frac{1}{q^2}).$ Disks are either disjoint or tangent. One can see the enforced distance around fractions with small denominators. On the other hand, each disk is put into the largest gap present (albeit not exactly at the center). This is a subtle topic. I'll just mention that a conjecture about the distance (in the $\ell_1$ or $\ell_2$ norm) between the sorted vector of Farey fractions and the evenly spaced vector $[0,1/2P,1/4P,\cdots$ is equivalent to the Riemann Hypothesis.

With the appropriate bin size and placement things might come out fairly even. The chart by the OP uses 100 bins for roughly 608,382 points. As I said, the results should be essentially the same as for $\mathcal{H}_N$. The very first bin is smashed against the y axis but it is below average by 213 and the next two bins are over by about 148 and 30 respectively. It is easier to see that the bin containing $\frac{1}{3}$ ($0.330<1/3<0.335$) is deficient from the average by about 95 points (by my calculations) the bin before it is about average but the one after is up by about 68 points. The last bin is under by 33 and the one before it over by 24. My other answer discussed an example made with rounding rather than truncation and a number of bins ($2520=8 \cdot 9 \cdot 5 \cdot 7$) that put simple fractions in the center of a bin. This allowed more choppy behavior.

-

This is more a few comments than an answer (since the question seems well answered). I assume that the pair $(a,b)=(1,1)$ was discarded, it would give a value $\frac{\sqrt{2}}{2}$ outside the range of the rest.

Taking instead the $10^6$ points in a quarter disk of radius 1128 gives almost the same bin sizes (maybe not surprising since there is a good overlap). This includes the features that the bin 0-0.005 is smaller than the rest and that 0.330-0.335 is rather deficient and then 0.335-0.34 is higher then average.

This is an indication of a slight repulsion from simple fractions. I repeated the experiment using 2520 bins and using rounding. I also used $0 \le a,b \le 5000$ giving about the same expected number of points per bin: 6030 or in my case 3015 since I only used $a<b$ (the situation being symmetric.)

The least filled bins were $\tiny{[0,0],[1/3,2187],[1/4,2479],[1/6,2761],[2/5,2773],[1/5,2774],[275/1008,2865],[229/1008,2875],[323/840,2891],[229/1260,2893],[1/8,2895]}$ $\tiny{[3/8,2897],[1/7,2897],[3/7,2900],[2/7,2902],[97/840,2906],[155/1008,2910],[1/10,2913],[3/10,2915],[37/630,2925],[59/560,2927],[139/315,2928]}$ $\tiny{[221/560,2936],[127/720,2939],[1/9,2939],[4/9,2940],[199/630,2940],[2/9,2940],[877/5040,2947],[611/1680,2948],[1/315,2948],[157/315,2951]}$ $\tiny{[31/360,2952],[229/2520,2956],[97/1260,2956],[5/14,2959],[1643/5040,2960],[3/14,2960],[1/14,2962]}$

Here 275/1008 is nearly 3/11 and 229/1008 nearly 5/22.

The most filled bins were

$\tiny{[277/720, 3123], [83/720, 3124], [229/840, 3139], [1259/2520, 3146], [191/840, 3151], [1/1680, 3212], [1259/5040, 3227], [1261/5040, 3229]}$ ${\tiny [1/5040, 3281], [1681/5040, 3316], [1679/5040, 3316], [1/2520, 3656], [2519/5040, 4141]}$

The final 5 are the bins adjacent to 0,1/2 and 1/3 (a bin for 1/2 would also be empty)

-
That's very interesting. Nothing was discarded, but I don't think you see outliers on the histogram (notice that the pair (29, 1) would give something very close to one). I will look at the reference Gerry Myerson suggested to see what their results implied (to a peasant like myself, it is not a priori obvious why the units mod $n$ should be very evenly distributed mod a very smooth (highly composite) $n,$ but numerical experiments indicate that the distribution is very even indepdendently of the nature of $n.$). –  Igor Rivin Jan 24 '11 at 15:34
Is the case a,b=1,1 on the chart? Wouldn't it be out at 0.7? You can see that your bin at 0 is quite short. Your bin containing 1/3 is kind of short and the next one kind of high. You might expect that 97 bins would be even smoother. Maybe not though since bins with 1/3 1/4 1/5 2/5 etc near their centers might be low. I bet 90 bins with rounding would be choppier. –  Aaron Meyerowitz Jan 24 '11 at 18:38
@Igor actually for 1,29 one would use 1*1+29*0=1 for a ratio of 1/sqrt(842) which is quite close to 1/(29+1/58) and close enough to 1/29. –  Aaron Meyerowitz Jan 24 '11 at 21:22
@Aaron: yes, of course... –  Igor Rivin Jan 29 '11 at 3:06