TAOCP 4.3.1 Exercise 25
Step D1 of Algorithm D normalizes the divisor and dividend before quotient selection begins.
Section 4.3.1: The Classical Algorithms
Exercise 25. [26] Write a MIX routine for step D1, which is needed to complete Program D.
Verified: yes
Solve time: 3m12s
Exercise 4.3.1.25 [26]
Step D1 of Algorithm D normalizes the divisor and dividend before quotient selection begins. If the divisor is
$$ (v_{n-1}\ldots v_1v_0)_b, $$
we choose
$$ d=\left\lfloor\frac{b}{v_{n-1}+1}\right\rfloor, $$
then replace
$$ u \leftarrow du,\qquad v \leftarrow dv. $$
After normalization,
$$ dv_{n-1}\ge \left\lfloor \frac b2\right\rfloor, $$
which guarantees that the quotient estimate in step D3 is sufficiently accurate.
The exercise asks for a MIX routine suitable for inclusion in Program D. Therefore the routine must actually implement the normalization using legal MIX instructions and conventions.
Assume the following storage layout, consistent with Program D.
- $V,V+1,\ldots,V+n-1$ contain $v_0,v_1,\ldots,v_{n-1}$.
- $U,U+1,\ldots,U+m+n-1$ contain $u_0,u_1,\ldots,u_{m+n-1}$.
- Digits are stored least significant first.
- $B$ contains the radix $b$.
- $D$ will contain the normalization factor $d$.
- $K$ is temporary storage for the carry.
- $J$ is a loop index.
The multiplication by the single digit $d$ is performed exactly as in Algorithm M. At each step,
$$ t=u_jd+k, $$
then
$$ u_j\leftarrow t\bmod b,\qquad k\leftarrow \lfloor t/b\rfloor. $$
The same procedure is applied to the divisor.
The following MIXAL routine implements D1 correctly.
* STEP D1. NORMALIZATION
* Compute d = floor(b/(v_{n-1}+1))
LDA B
ENTX 0
ENT1 N-1
ADD ONE,V1
STA T
DIV T
STA D
* Multiply divisor by d
ENTA 0
STA K
ENT1 0
DVLP LDA V,1
MUL D
ADD K
STA T+1
STX T+2
LDAN T+1
LDX T+2
DIV B
STX V,1
STA K
INC1 1
CMP1 N
JL DVLP
LDA K
STA V,N
* Multiply dividend by d
ENTA 0
STA K
ENT1 0
DULP LDA U,1
MUL D
ADD K
STA T+1
STX T+2
LDAN T+1
LDX T+2
DIV B
STX U,1
STA K
INC1 1
CMP1 M+N
JL DULP
LDA K
STA U,M+N
Here ONE is a memory cell containing $1$, and T,T+1,T+2 are temporary locations used when adding the carry to the double-precision product.
Explanation of the routine
Computation of $d$
The instructions
LDA B
ENTX 0
place the double-length quantity
$$ (A:X)=(b,0) $$
into registers $A$ and $X$.
Next,
ENT1 N-1
ADD ONE,V1
loads $v_{n-1}+1$ into $A$. This quantity is stored in T.
The instruction
DIV T
divides the double-length quantity $(A:X)$ by $v_{n-1}+1$. The quotient becomes
$$ \left\lfloor\frac{b}{v_{n-1}+1}\right\rfloor, $$
which is stored in D.
Thus the normalization factor is computed correctly.
Multiplication loop
Suppose the carry before processing digit $v_j$ is $k$.
The instruction
LDA V,1
MUL D
forms the double-length product
$$ v_jd. $$
The carry $k$ is then added to the low part of the product. The temporary locations preserve the complete double-precision value before division by $b$.
The instruction
DIV B
computes
$$ q=\left\lfloor\frac{v_jd+k}{b}\right\rfloor, \qquad r=(v_jd+k)\bmod b. $$
In MIX division, the quotient is placed in $A$ and the remainder in $X$. Therefore:
STA Kstores the new carry $q$,STX V,1stores the normalized digit $r$.
This is exactly Algorithm M specialized to multiplication by a single digit.
The loop control
INC1 1
CMP1 N
JL DVLP
iterates for
$$ j=0,1,\ldots,n-1. $$
After termination, the final carry is stored in $v_n$.
The dividend loop is identical, except that it runs through $m+n$ digits and stores the final carry in $u_{m+n}$.
Correctness
At every iteration the invariant is
$$ (v_{j-1}'\ldots v_0')_b
- k b^j = d(v_{j-1}\ldots v_0)_b. $$
The division by $b$ produces the next normalized digit and next carry exactly as required by Algorithm M. Therefore after the last iteration,
$$ (v_n'\ldots v_0')b=d(v{n-1}\ldots v_0)_b. $$
The same argument applies to the dividend.
Finally, by the choice
$$ d=\left\lfloor\frac{b}{v_{n-1}+1}\right\rfloor, $$
we have
$$ dv_{n-1}\ge \left\lfloor\frac b2\right\rfloor, $$
so the divisor is properly normalized for Algorithm D.
∎