Rounding to the nearest remainder

The oldest algorithm, drawn as a tiling showed Euclid’s algorithm cutting squares from a rectangle: from an $a \times b$ rectangle remove as many $b \times b$ squares as fit, then repeat on what is left, until nothing is left. Each round is a division with remainder, $a = qb + r$ with $0 \le r < b$, and the number of rounds — the number of division steps — is the algorithm’s cost. That essay measured it: on pairs of numbers up to $N$, about $(12 \ln 2/\pi^2)\ln N$ steps on average, a result of Hans Heilbronn and John Dixon, and at worst about $\log_\varphi N$, reached by consecutive Fibonacci numbers.

There is a choice hidden in the algorithm that the tiling makes invisible. The remainder $r$ is what is left after rounding the quotient $a/b$ down. Round it up instead, take one square too many, and the remainder is $b - r$ with its sign reversed; for finding a greatest common divisor the sign does not matter. So at every step there are two candidates, $r$ and $b - r$, and the usual algorithm always takes the first. This essay takes the smaller of the two, which is the remainder of rounding to the nearest, and measures what that buys.

The least-remainder algorithm #

Rounding to the nearest means taking whichever of $r$ and $b - r$ is smaller, so the new remainder is at most $b/2$. The least-remainder algorithm is the usual one with that single change: replace $(a, b)$ by $(b, \min(r, b - r))$ and repeat until the remainder is nought. It computes the same greatest common divisor, because $b - r$ and $r$ differ from multiples of $b$ by the same amount and so share every common divisor of $a$ and $b$.

The figure measures it on sixty thousand random pairs at each size from a hundred to a hundred million. Both averages grow in a straight line against the number of digits. The usual algorithm adds 1.94 steps for every factor of ten, matching $12 \ln 2 \cdot \ln 10 / \pi^2 = 1.941$. The least-remainder algorithm adds 1.35, matching $12 \ln\varphi \cdot \ln 10/\pi^2 = 1.347$, where $\varphi$ is the golden ratio. The two formulas differ only in the logarithm on top, and the ratio of the rates is

$$\frac{\ln \varphi}{\ln 2} \approx 0.694.$$

Rounding to the nearest saves, asymptotically, about thirty per cent of the steps.

A Fibonacci pair, two ways #

The saving is largest exactly where the usual algorithm is slowest.

For consecutive Fibonacci numbers every quotient is one: $1597 = 1 \times 987 + 610$, $987 = 1 \times 610 + 377$, and so on down the Fibonacci sequence, one term per step, fifteen steps in all. Rounding to the nearest sees that $1597/987 \approx 1.62$ is closer to 2 than to 1, takes two squares where only one fits, and is left with $1597 - 2 \times 987 = -377$: a remainder of 377, two Fibonacci numbers further down. Every subsequent step does the same, and the algorithm finishes in eight steps.

That is the mechanism of the whole saving. A quotient of one is the usual algorithm’s least efficient step, removing a single square and leaving a remainder more than half the divisor. Rounding to the nearest never leaves a remainder above half, so it never makes a step that small, and it removes the quotients of one by merging each with its neighbour. A fraction that never closes described the sequence of quotients as the continued fraction of $a/b$; the least-remainder algorithm computes a different expansion, the nearest-integer continued fraction, whose digits may be negative and are never one in absolute value after the first.

Never slower, and never beaten #

Taking the smaller remainder is a greedy choice: it makes the current number as small as possible without looking ahead. Greedy choices are often wrong, and a smaller remainder now could, in principle, lead to a longer run of steps later. It does not.

The figure checks two things for every pair with the larger number up to 400. First, the least-remainder algorithm never takes more steps than the usual one: the savings run from nought to five, and no pair has a negative saving. Second, and more strongly, no mixture of rounding up and rounding down does better. The fewest steps over all such mixtures is computed by a table over every pair — the cost of a pair is one plus the cheaper of the two pairs it can lead to — and it equals the least-remainder count for all 80,200 pairs.

Leopold Kronecker proved this in his lectures on number theory, published in 1901: of all the Euclidean algorithms that choose, at each step, between the two remainders, the one choosing the smaller is the shortest, for every input. Earlier, Karl Vahlen had studied the expansions that the free choices produce. A small case shows what is being compared. For the pair 13 and 8 the usual algorithm passes through 5, 3, 2, 1 and 0, five steps. Rounding to the nearest at the first step gives $13 = 2 \times 8 - 3$, remainder 3, then $8 = 3 \times 3 - 1$, remainder 1, then nought: three steps. Every other mixture of choices — round down first and to the nearest after, or the reverse — gives three, four or five steps, never two, and the table confirms the same for each of the eighty thousand pairs. The greedy rule wins because a smaller remainder never makes the rest of the computation longer: the number of steps needed from a pair is, roughly, a measure of how slowly the ratio of the pair can be approximated by fractions, and making the second number smaller can only make that ratio easier.

The slowest inputs: Fibonacci, and Pell #

For the usual algorithm, the slowest pair up to any size is a pair of consecutive Fibonacci numbers, which is Gabriel Lamé’s theorem of 1844. For the least-remainder algorithm the slowest pairs are different.

The record pairs for the least-remainder algorithm, found by checking every pair up to ten thousand, are $3/2$, $7/5$, $17/12$, $41/29$, $99/70$, $239/169$, $577/408$, $1393/985$, $3363/2378$ and $8119/5741$ — the convergents of $\sqrt2$, built from the Pell numbers $1, 2, 5, 12, 29, 70, 169, \ldots$, each twice the last plus the one before. Its worst case grows like $\log_{1 + \sqrt2} N$ rather than $\log_\varphi N$. The result goes back to Dupré in 1846, two years after Lamé.

The pattern has a clean explanation. Each algorithm is a way of expanding a ratio as a continued fraction, and its worst case is the ratio whose expansion is as slow as possible: for the usual algorithm, all quotients one, which is the golden ratio, $1 + 1/(1 + 1/(1 + \cdots))$; for the least-remainder algorithm, whose quotients cannot be one, all quotients two, which is $1 + \sqrt2 = 2 + 1/(2 + 1/(2 + \cdots))$. The Pell numbers are to $\sqrt2$ what the Fibonacci numbers are to the golden ratio, and the essays on one solution that makes all the others met them as the solutions of Pell’s equation $x^2 - 2y^2 = \pm 1$.

The tiling, with squares taken away #

The rectangle picture of the first essay can be redrawn for the new algorithm, and the redrawing shows what rounding up means geometrically. The usual algorithm cuts from an $a \times b$ rectangle the $q$ squares of side $b$ that fit, leaving a strip of width $r$. Rounding up cuts $q + 1$ squares, which do not fit: the last one overhangs the rectangle by $b - r$. The overhang is a strip of width $b - r$ lying outside, and the algorithm continues on the rectangle $b \times (b - r)$ as if the overhang were a hole to be filled.

So the least-remainder algorithm tiles a rectangle with squares that are sometimes added and sometimes subtracted — a signed tiling — and its rule is to overshoot whenever the overshoot is smaller than the shortfall. The total of the signed squares’ areas is still $ab$, and the sequence of signed quotients is the nearest-integer continued fraction. The player who meets the first long run found in Euclid’s game a reason to care which remainder is left at every step; there the player chooses how many squares to remove and the choice decides the winner, and here a fixed rule chooses and the choice decides the cost.

The signed picture also explains why the worst case changes. A Fibonacci rectangle is the usual algorithm’s slowest because every step removes exactly one square and leaves a rectangle of the same shape, scaled down by the golden ratio. In the signed tiling, one step on a Fibonacci rectangle takes two squares, overshooting, and lands two Fibonacci numbers down, so the shape repeats after a shrinkage of $\varphi^2$ per step rather than $\varphi$. The shape that repeats after the least shrinkage in the signed tiling is the one with quotient two at every step, the rectangle of ratio $1 + \sqrt2$, which loses two squares and leaves a rectangle of the same shape, scaled down by $1 + \sqrt 2$.

The same choice in polynomials #

The choice between rounding down and to the nearest has an analogue for polynomials, where it disappears. The remainders that count the roots ran Euclid’s algorithm on polynomials, dividing one by another and keeping the remainder, and used the sequence of remainders — with their signs reversed, in Sturm’s version — to count real roots. For polynomials there is no rounding: the quotient is determined exactly, a polynomial of the right degree, and the remainder has strictly smaller degree. The algorithm’s length is fixed by the degrees, and no choice could shorten it.

What survives the change of setting is the sign. Sturm’s sequence negates each remainder, which is exactly the rounding-up choice made at every step, and for polynomials it changes nothing about the length and everything about what the sequence can be used for: the sign changes along the sequence count roots. For whole numbers the negation of a remainder is a choice about cost, and for polynomials it is a choice about meaning, and what the signs allow followed that second use as far as Descartes’s rule.

Why the averages have those constants #

The average-case constants come from the same picture, applied to typical rather than extreme ratios. A random ratio’s continued fraction digits follow the Gauss–Kuzmin distribution, under which the average of the logarithm of the denominators’ growth per step is $\pi^2/(12 \ln 2)$ — Lévy’s constant. Dividing $\ln N$ by that growth per step gives the average number of steps, $(12 \ln 2/\pi^2)\ln N$, which is where the usual algorithm’s constant comes from.

The least-remainder algorithm’s denominators grow faster per step, because it never takes the small steps, and the corresponding constant for the nearest-integer continued fraction, which Adolf Hurwitz studied in 1889, is $\pi^2/(12 \ln \varphi)$. The golden ratio enters through the map that generates the expansion — take the fractional part to the nearest, invert, repeat — whose long-run distribution of values has a density built from $\varphi$, just as the usual expansion’s Gauss–Kuzmin density is built from $\ln 2$. The nearest-integer expansion of a number can also be obtained from its ordinary expansion by merging every quotient of one into its neighbour, which is how the two densities are related. So the golden ratio is the worst case for one algorithm and the source of the average for the other, and $\sqrt2$, the other algorithm’s worst case, does not appear in either average.

How the counts spread #

The averages describe the centre of a distribution that is itself worth seeing.

For pairs up to a million the usual algorithm takes $11.71$ steps on average with a spread of $2.70$; the least-remainder algorithm takes $8.52$ with a spread of $1.70$. Both distributions are bell-shaped. Doug Hensley proved in 1994 that the usual algorithm’s step count is asymptotically normal, with a variance growing in proportion to $\ln N$, and the same holds for the least-remainder version with different constants. A step count is a sum of many nearly independent contributions, one per quotient, and a bell curve assembled out of coin flips is the general reason sums of that kind spread into bells. The nearest-rounding bell is narrower because it is a sum of fewer terms: fewer steps, less room for fluctuation.

The limit arrives slowly #

The ratio of the two averages tends to $\ln\varphi/\ln 2 = 0.694$, but not quickly.

At a hundred the least-remainder algorithm takes 79 per cent of the usual algorithm’s steps; at a hundred million, 72 per cent; extended to $10^{16}$ by the two straight-line fits, 71. Each average is a straight line in $\ln N$ with its own intercept, and the intercepts differ — the least-remainder algorithm’s is larger, which comes from how the first few steps of the two expansions differ — so the ratio approaches its limit only as $\ln N$ grows large enough to swamp them. For every size at which anybody computes, the saving is a little under the asymptotic thirty-one per cent.

Why computers usually do not bother #

The least-remainder algorithm is never slower in steps, and yet the usual algorithm is the one in every textbook and most libraries. The reason is the cost of a step. Computing $\min(r, b - r)$ needs a comparison and possibly a subtraction on top of the division, and on most hardware division dominates the cost of a step anyway; saving thirty per cent of the divisions saves roughly thirty per cent of the time, which is worth having but rarely decisive. Binary versions of the algorithm, which replace division by shifts and subtractions, change the accounting again: there the cost of a step is tiny and the number of steps is governed by the binary length of the numbers rather than by their continued fractions, and the question of which remainder to keep becomes a question about which of two odd numbers has more factors of two in its difference. For very large numbers, the fastest methods use neither: they work on the leading digits only, deciding several quotients at once, and the choice between rounding down and to the nearest dissolves into the half-gcd methods of Lehmer, Knuth and Schönhage.

The idea mattered most where steps were expensive and done by hand. A method that is allowed to miss found it inside the chakravala, the Indian method for Pell’s equation, whose choice of the multiplier nearest to $\sqrt D$ is a nearest-integer continued fraction centuries before the name; and the same rounding appears in reducing quadratic forms and in finding short lattice vectors, where Gauss’s reduction is Euclid’s algorithm on vectors with the nearest multiple subtracted.

What the counts do and do not establish #

The average curves are sampled — sixty thousand seeded pairs per size — so their constants carry sampling error of a hundredth of a step, small against the slopes measured. The never-slower and never-beaten checks are exhaustive, but only up to 400; Kronecker’s theorem is what extends them to every pair. The worst-case records are exhaustive up to ten thousand, and Dupré’s theorem is what says the Pell pattern continues.

The figures also leave out a variant that complicates the clean picture: algorithms that may round either way by any rule, not just to the nearest — for instance, choosing the remainder that makes the next quotient large. Kronecker’s theorem covers all of these, since it compares against every possible choice, which is what the table of best choices checks.

Still open: the exact constants of the variance #

The average number of steps is known exactly for both algorithms, including the constant term in the usual algorithm, Porter’s constant. The variance is known to grow in proportion to $\ln N$, and its leading constant for the usual algorithm is an explicit but complicated number, Hensley’s constant, computed only numerically and through the spectral theory of a transfer operator associated with the Gauss map; for the least-remainder algorithm and the other variants, the corresponding constants are less well determined. Whether any of them has a closed form in terms of familiar constants — as the averages do, through $\pi$, $\ln 2$ and $\ln\varphi$ — is not known, and the operator whose eigenvalue they come from has resisted every attempt to diagonalise it in closed form.