\documentstyle{article}
\begin{document}
\title{Integer Division}
\author{George E. Collins}
\maketitle

Analagous to the program DPR that is needed for multiplication,
SACLIB has an assembly language program DQR (digit quotient
and remainder) that takes as input a two-$\beta$-digit
dividend and a $\beta$-digit divisor, and produces as output
two $\beta$-digits, a quotient and a remainder.  DQR
computes the quotient one bit at a time by subtraction
and shifting.  Its computing time on ren is approximately
0.94 microseconds.

The SACLIB program IQR divides an arbitrary integer $A$ by an
arbitrary nonzero integer $B$. obtaining a quotient $Q$ and
remainder $R$ such that $Q$ is the floor function of $A / B$ 
and $R = A - B Q$.  If $B$ is a $\beta$-digit, IQR uses
the special subprogram IDQR, which we will not discuss.
Otherwise, $|A|$ and $|B|$ are copied into arrays using
IFLTA and the method described in Section~4.3.1, Algorithm~D,
of \cite{Knu97} is used.  Note that if $Q'$ is the integral
part of $|A| / |B|$ and $R'$ = $|A| - |B| Q'$ then $Q = Q'$
if sign$(A) = $sign$(B)$ and $Q = - Q'$ otherwise.  Also, 
sign$(R) = $sign$(A)$.

Step 4 of IQR detects cases where the quotient is obviously 0.
Theorem A, page 271, of \cite{Knu97} shows that if the
leading $\beta$-digit of the divisor is at least $\beta / 2$,
then dividing the two high-order digits of the divident by
the leading digit of the divisor will produce a number that
is greater than or equal to the leading digit of the quotient.
Theorem B, immediately following, shows that if the leading
$\beta$-digit is at least $\beta / 2$, then the quotient digit
obtained will be at most 2 larger than the correct quotient
digit.  Therefore Step 5 of IQR "normalizes" the divisor,
multiplying it by the appropriate power of 2 to achieve this
condition.  Notice that here we have an application of PDLOG2,
that the multiplication is done by shifting, and that the
dividend must also be multiplied by the same power of 2.
The quotient digit approximation determined by this two-digit
by one-digit division is computed in the first 10 lines of
Step 10 of IQR.  The division is performed by PNDQR, not DQR. 
On a prior architecture PNDQR was faster than DQR, but on ren
the times of the two programs are nearly the same. 
 Next the method of Step~D3 of Algorithm~D of
\cite{Knu97} is employed.  This method adjusts this first
aproximation to become the integer part of the three highest
order $\beta$-digits of the dividend divided by the two
highest $\beta$-digits of the divisor.  The adjustment is
made without division since Theorem~B guarantees that the
adjusted quotient is larger by at most two.  Now we know that
our approximate quotient digit is either correct or too large
by only one.  Exercise~21 at the end of Section~4.3.1 shows
that it is too large with a probability of at most $2 / \beta$.

Step 11 multiplies the divisor by the quotient digit approximation
$q$ and subtracts the product from the dividend, carrying out the
multiplication and subtraction as a single process.  If the
result is negative, Step 12 adds the divisor back and decreases
$q$ by 1.  After all quotient digits have been computed,
Step~14 "denormalizes", dividing the remainder by the power
of 2 by which we multiplied at the beginning, and imposing the
correct signs on quotient and remainder.  Also, the quotient
and remainder are converted to lists.

Table~\ref{table-iqr} shows times T in microseconds for using IQR
to divide a $2 n$ $\beta$-digit integer by an $n$ $\beta$-digit
divisor.

For nonzero integers $A$ and $B$ of lengths $m$ and $n$ respectively,
with $m \geq n$, the computing time of IQR is codominant with 
$n(m - n + 1)$.  If $m < n$, its computing time is codominant
with $n$.

The program IEQ (integer exact quotient) can be used in lieu of
IQR if it is known in advance, as it often is, that $A$ is exactly
divisible by $B$.  The method used by IEQ is due to T. Jebelean,
\cite{Jeb93}.  In this method the low-order digits of the quotient 
are computed first.  Here we will describe the method for the case
that the number system radix is a power of 2, as in SACLIB, but
it can be used more generally.  First we arrange for $B$ to be
odd, as follows.  If $B$ is divisible by $2^k$, with k maximal,
then $A$ is also divisible by $2^k$ since $A | B$.  So we divide
both $A$ and $B$ by $2^k$ (by shifting).  This does not change
the quotient $Q$, and now $B$ is odd.  Let $b_0$ be the low order
$\beta$-digit of $B$.  $b_0$ is odd so it is relatively to $\beta$,
and has an inverse, $b_0'$ modulo $\beta$.  Let $a_0$ be the low
order digit of $A$ and let $c_0 = a_0 b_0'$ modulo $\beta$.  $c_0$ 
must be the low order digit of $Q$, and if we subtract $c_0 B$ from
$A$, the result has 0 as low order digit, and is divisible by
$\beta$.  If we call this result $A'$, we then have $A = A' \beta
+ c_0 B$.  The next digit of $Q$ can be obtained in the same way
with $A'$ in place of $A$.  The method as described thus far has
a computing time that is the same, to within codominance, as IQR.
But it can be made faster if we observe the following.  When we
multiply $B$ by $c_0$ and subtract from $A$, we need not, in general,
carry out the process to completion as we proceed from low order
digits to higher order digits.  We only need to proceed through
$a_{m-n}$ because the quotient will have at most $m - n + 1$ digits.
Any higher order digits of $A$ would never be used in the algorithm
so they need not be computed.  However, if $m - n + 1 \geq n$ then
no saving is achieved in computing $c_0$.  Jebelean shows that the
number of digit multiplications used by this method is

\begin{tabular}{ll}
$(m - n + 1)(m - n + 2)/2$ & if $m + 1 \leq 2 n$, \\
$m n - 3 n(n - 1)/2$ & otherwise.
\end{tabular}

Table~\ref{table-ieq} shows the times required by IEQ to divide
a $2$ $n \beta$-digit number by an $n$ $\beta$-digit number.  As
expected, these times are approximately one half of the
corresponding times for IQR.

\begin{table}
\caption{IQR Times}
\label{table-iqr}
\begin{center}
\begin{tabular}
{| r      |   r    |} \hline
\multicolumn{1}{|c|}{n} &
\multicolumn{1}{|c|}{T} \\ \hline
   1    &  2.59  \\ \hline
   2    &  6.32  \\ \hline
   3    & 10.3   \\ \hline
   4    & 14.5   \\ \hline
   5    & 17.2   \\ \hline
  10    & 47.2   \\ \hline
  20    &141     \\ \hline
  40    &478     \\ \hline
\end{tabular}
\end{center}
\end{table}

\begin{table}
\caption{IEQ Times}
\label{table-ieq}
\begin{center}
\begin{tabular}
{| r      |   r    |} \hline
\multicolumn{1}{|c|}{n} &
\multicolumn{1}{|c|}{T} \\ \hline
   1    &  2.25  \\ \hline
   2    &  3.97  \\ \hline
   3    &  5.41  \\ \hline
   4    &  7.12  \\ \hline
   5    &  9.67  \\ \hline
  10    & 23.5   \\ \hline
  20    & 66.9   \\ \hline
  40    &228.    \\ \hline
\end{tabular}
\end{center}
\end{table}

\begin{thebibliography}{99}

\bibitem{Jeb93}
Tudor Jebelean.
An Algorithm for Exact Division.
J. Symbolic Computation,
Vol. 15 (1993), pp.169-180.

\bibitem{Knu97}
Donald E. Knuth.
The Art of Computer Programming, Vol. 2, 3rd Ed.
Addison-Wesley, 1997

\end{thebibliography}

\end{document}
