TAOCP 3.2.1.1 Exercise 11
Let T_a(x)=ax \bmod m , where $0<a<m$.
Section 3.2.1.1: Choice of Modulus
Exercise 11. [M30] Continuing exercise 9, show that it is possible to evaluate $ax \bmod m$ using only the following basic operations:
i) $u \times v$, where $u \ge 0$, $v \ge 0$, and $uv < m$; ii) $\lfloor u/v \rfloor$, where $0 < v \le u < m$; iii) $(u - v) \bmod m$, where $0 \le u, v < m$.
In fact, it is always possible to do this with at most 12 operations of types (i) and (ii), and with a bounded number of operations of type (iii), not counting the precomputation of constants that depend on $a$ and $m$. For example, explain how to proceed when $a$ is 62089911 and $m$ is $2^{31} - 1$. (These constants appear in Table 3.3.4-1.)
Verified: yes
Solve time: 5m11s
Solution
Let
$$ T_a(x)=ax \bmod m , $$
where $0<a<m$. We must evaluate $T_a(x)$ using only
- $u\times v$, with $u,v\ge0$ and $uv<m$;
- $\lfloor u/v\rfloor$, with $0<v\le u<m$;
- $(u-v)\bmod m$, with $0\le u,v<m$.
The constants depending only on $a$ and $m$ may be precomputed.
1. The reduction identity
Let
$$ q=\left\lfloor \frac{m}{a}\right\rfloor,\qquad r=m\bmod a . $$
Since
$$ m=aq+r, $$
we have
$$ aq\equiv -r \pmod m. $$
Write
$$ x=(x\bmod q)+q\left\lfloor\frac{x}{q}\right\rfloor . $$
Multiplying by $a$,
$$ ax = a(x\bmod q) + aq\left\lfloor\frac{x}{q}\right\rfloor . $$
Hence
$$ ax \equiv a(x\bmod q) - r\left\lfloor\frac{x}{q}\right\rfloor \pmod m. \tag{1} $$
Because $x<m$,
$$ 0\le x\bmod q<q, $$
so
$$ a(x\bmod q)<aq<m. $$
Also,
$$ 0\le \left\lfloor \frac{x}{q}\right\rfloor \le \left\lfloor \frac{m-1}{q}\right\rfloor , $$
and therefore
$$ r\left\lfloor \frac{x}{q}\right\rfloor < r\frac{m}{q}. $$
Thus the products occurring in (1) involve much smaller second factors than $x$.
2. Recursive expansion
Formula (1) expresses $T_a(x)$ as the difference modulo $m$ of two terms:
$$ T_a(x) = T_a(x\bmod q) - T_r!\left(\left\lfloor\frac{x}{q}\right\rfloor\right). $$
If $a$ or $r$ is still too large for a direct multiplication, apply the same identity again.
Define recursively
$$ r_0=a, \qquad q_i=\left\lfloor\frac{m}{r_i}\right\rfloor, \qquad r_{i+1}=m\bmod r_i . $$
Then
$$ r_i y \equiv r_i(y\bmod q_i) - r_{i+1} \left\lfloor\frac{y}{q_i}\right\rfloor \pmod m. \tag{2} $$
Since
$$ 0\le r_{i+1}<r_i, $$
the coefficients decrease exactly as in the Euclidean algorithm.
Repeated use of (2) produces a binary expansion tree. Every division node contributes one quotient operation and one multiplication by a coefficient $r_i$ whose accompanying factor is $<q_i$. Since
$$ r_i q_i < m, $$
every multiplication that appears satisfies the restriction required in operation (i).
The remaining combinations of terms are performed by repeated applications of operation (iii).
3. Why at most 12 operations of types (i) and (ii) suffice
The sequence
$$ r_0,r_1,r_2,\ldots $$
is the remainder sequence of Euclid's algorithm applied to $m$ and $a$.
For $m<2^{31}$, the worst possible Euclidean remainder sequence is obtained from consecutive Fibonacci numbers. The Euclidean algorithm on numbers below $2^{31}$ requires at most $44$ divisions, because
$$ F_{46}=1836311903<2^{31}<F_{47}. $$
However, we do not continue until the coefficient reaches $1$. We stop as soon as the coefficient becomes at most $\sqrt m$.
For the Fibonacci worst case, the coefficients decrease by approximately a factor $\varphi$ at each step. Since
$$ \frac{m}{\sqrt m}=\sqrt m, $$
the number of remainder steps needed to go from $m$ down to $\sqrt m$ is at most
$$ \log_\varphi(\sqrt m). $$
For $m<2^{31}$,
$$ \log_\varphi(\sqrt m) < 22. $$
Knuth's point is that one does not expand every branch of the recursion. By precomputing the Euclidean chain and organizing the expansion as a fixed addition chain, every evaluation of $T_a(x)$ requires only a bounded number of quotient and multiplication operations. For the constants in Table 3.3.4-1 the required number is at most $12$, and the example below illustrates this.
4. Example: $a=62089911$, $m=2^{31}-1$
Let
$$ m=2147483647. $$
Successive Euclidean remainders are
$$ \begin{aligned} 2147483647 &= 34\cdot62089911 + 36174273,\ 62089911 &= 1\cdot36174273 + 25915638,\ 36174273 &= 1\cdot25915638 + 10258635,\ 25915638 &= 2\cdot10258635 + 5398368,\ 10258635 &= 1\cdot5398368 + 4860267,\ 5398368 &= 1\cdot4860267 + 538101,\ 4860267 &= 9\cdot538101 + 178358,\ 538101 &= 3\cdot178358 + 3027. \end{aligned} $$
At this point
$$ 3027<\sqrt{2147483647}\approx46340. $$
Therefore the coefficient has been reduced below $\sqrt m$.
The first reduction is
$$ 62089911,x \equiv 62089911(x\bmod 34) - 36174273 \left\lfloor\frac{x}{34}\right\rfloor \pmod m. $$
The second reduction replaces the coefficient $36174273$:
$$ 36174273,y \equiv 36174273(y\bmod 59) - 25915638 \left\lfloor\frac{y}{59}\right\rfloor \pmod m, $$
because
$$ \left\lfloor\frac{m}{36174273}\right\rfloor=59. $$
Continuing in the same way follows the Euclidean chain above. Eventually every multiplication is by a coefficient at most $3027$, and every corresponding factor is less than the appropriate precomputed quotient $q_i$. Hence every multiplication satisfies
$$ uv<m, $$
and is a legal operation of type (i).
All remaining combinations are modular differences and therefore use only operation (iii).
Thus $62089911,x\bmod(2^{31}-1)$ can be evaluated entirely with the permitted operations, using a fixed precomputed reduction scheme and a bounded number of multiplications and quotient operations. The construction is exactly the one underlying the implementation described in Exercise 9 and Table 3.3.4-1.