TAOCP 4.5.2 Exercise 38
Algorithm B is the binary gcd algorithm.
Section 4.5.2: The Greatest Common Divisor
Exercise 38. ▶ [**] [M32] [R. W. Gosper.] Demonstrate how to modify Algorithm B for large numbers, using ideas analogous to those in Algorithm L.
Verified: yes
Solve time: 5m02s
Exercise 4.5.2.38 [M32]
Algorithm B is the binary gcd algorithm. Its advantage is that it replaces division by repeated halving and subtraction. For large multiprecision integers, however, the subtraction step can still be expensive when many iterations are required. The exercise asks for a modification analogous to Algorithm L, which accelerates Euclid's algorithm by using information from the leading digits to perform several quotient steps at once.
The corresponding modification of Algorithm B is to replace long sequences of subtraction steps by a single multiprecision subtraction using a quotient estimated from the leading words of the operands.
Suppose that after the removal of common factors of 2, we have odd positive integers $u$ and $v$. Assume $u>v$.
Write
$$ u=(u_{n-1}\cdots u_0)\beta,\qquad v=(v{m-1}\cdots v_0)_\beta, $$
where $\beta$ is the machine radix and $u_{n-1},v_{m-1}\neq0$.
In Algorithm B, the basic reduction is
$$ u\leftarrow u-v, $$
followed by the removal of powers of $2$. If $u\gg v$, many such reductions may occur before the leading part of $u$ changes appreciably. The analogue of Algorithm L is therefore:
B6'.
Using only the leading one or two words of $u$ and $v$, compute an approximate quotient
$$ q\approx \left\lfloor \frac{u}{v}\right\rfloor . $$
More precisely, let
$$ \hat u=u_{n-1}\beta+u_{n-2},\qquad \hat v=v_{m-1}\beta+v_{m-2}, $$
with missing digits interpreted as $0$. Form
$$ q=\left\lfloor \frac{\hat u}{\hat v}\right\rfloor , $$
and replace $q$ by $1$ if this estimate is $0$.
Then perform
$$ u\leftarrow u-qv . $$
If the result is negative, decrease $q$ until the result is nonnegative, exactly as in the quotient-adjustment stage of Algorithm L.
Afterwards remove all factors of $2$ from the new value of $u$:
$$ u\leftarrow u/2^{\nu_2(u)}, $$
where $\nu_2(u)$ denotes the exponent of the highest power of $2$ dividing $u$.
Finally interchange $u$ and $v$ if necessary so that $u\ge v$, and repeat.
Correctness
The fundamental invariant is
$$ \gcd(u,v)=\gcd(u-qv,v) $$
for every integer $q$. Hence the replacement
$$ u\leftarrow u-qv $$
preserves the gcd.
The subsequent division by powers of $2$ is also valid because $u$ and $v$ are odd at the beginning of each main cycle. Therefore
$$ \gcd(u,v) =\gcd!\left(\frac{u}{2^{\nu_2(u)}},v\right). $$
Thus every iteration preserves the same gcd as Algorithm B.
Termination follows because $0\le u-qv<u$. After the removal of factors of $2$, the new value of $u$ is still strictly smaller than the previous one. Therefore the sequence of positive operands decreases and must eventually reach $0$. At that point the remaining nonzero operand is the gcd.
Relation with Algorithm L
Algorithm L accelerates Euclid's algorithm by using the leading digits of the operands to predict several quotient steps before carrying out a full multiprecision operation. The modification above does exactly the same thing for Algorithm B.
Ordinary Algorithm B repeatedly performs
$$ u\leftarrow u-v,\quad u\leftarrow \frac{u}{2^{\nu_2(u)}}, $$
often many times in succession. The modified version replaces a whole block of such subtraction steps by a single transformation
$$ u\leftarrow u-qv, $$
where $q$ is estimated from the leading words, just as Algorithm L estimates Euclidean quotients from leading digits.
Hence the large-number version of Algorithm B is obtained by grouping many binary-subtraction iterations into one multiprecision step, using leading-word quotient estimates exactly analogous to those used in Algorithm L. This preserves the gcd invariant while substantially reducing the number of full-length multiprecision operations required. ∎