# Division algorithm

{{#invoke:Hatnote|hatnote}} A division algorithm is an algorithm which, given two integers N and D, computes their quotient and/or remainder, the result of division. Some are applied by hand, while others are employed by digital circuit designs and software.

Division algorithms fall into two main categories: slow division and fast division. Slow division algorithms produce one digit of the final quotient per iteration. Examples of slow division include restoring, non-performing restoring, non-restoring, and SRT division. Fast division methods start with a close approximation to the final quotient and produce twice as many digits of the final quotient on each iteration. Newton-Raphson and Goldschmidt fall into this category.

• Q = Quotient
• N = Numerator (dividend)
• D = Denominator (divisor).

## Division by repeated subtraction

The simplest division algorithm, historically incorporated into a greatest common divisor algorithm presented in Euclid's Elements, Book VII, Proposition 1, finds the remainder given two positive integers using only subtractions and comparisons:

 

while N ≥ D do N := N - D end return N 

The proof that the quotient and remainder exist and are unique, described at Euclidean division, gives rise to a complete division algorithm using additions, subtractions, and comparisons:

 

function divide(N, D) if D == 0 then throw DivisionByZeroException end if D < 0 then (Q,R) := divide(N, -D); return (-Q, R) end if N < 0 then (Q,R) := divide(-N, D) if R = 0 then return (-Q, 0) else return (-Q - 1, D - R) end end // At this point, N ≥ 0 and D > 0 Q := 0; R := N while R ≥ D do Q := Q + 1 R := R - D end return (Q, R) end 

This procedure always produces R ≥ 0. Although very simple, it takes Ω(Q) steps, and so is exponentially slower than even slow division algorithms like long division. It is useful if Q is known to be small (being an output-sensitive algorithm), and can serve as an executable specification.

## Long division

{{#invoke:main|main}} Long division is the standard algorithm used for pen-and-paper division of multidigit numbers expressed in decimal notation. It shifts gradually from the left to the right end of the dividend, subtracting the largest possible multiple of the divisor at each stage; the multiples become the digits of the quotient, and the final difference is the remainder. When used with a binary radix, it forms the basis for the integer division (unsigned) with remainder algorithm below. Short division is an abbreviated form of long division suitable for one-digit divisors.

## Integer division (unsigned) with remainder

The following algorithm, the binary version of the famous long division, will divide N by D, placing the quotient in Q and the remainder in R. All values are treated as unsigned integers.{{ safesubst:#invoke:Unsubst||date=__DATE__ |$B= {{#invoke:Category handler|main}}{{#invoke:Category handler|main}}[citation needed] }}   if D == 0 then throw DivisionByZeroException end Q := 0 initialize quotient and remainder to zero R := 0 for i = n-1...0 do where n is number of bits in N R := R << 1 left-shift R by 1 bit R(0) := N(i) set the least-significant bit of R equal to bit i of the numerator if R >= D then R = R - D Q(i) := 1 end end  ### Example If we take N=1100 and D=100 Step 1: Set R=0 and Q=0 Step 2: Take i=3 (bits in N, -1) Step 3: R=00 (left shifted by 1) Step 4: R=01 (setting R(0) to N(i)) Step 5: R<D , so skip statement Step 2: Set i=2 Step 3: R=010 Step 4: R=011 Step 5: R<D , statement skipped Step 2: Set i=1 Step 3: R=0110 Step 4: R=0110 Step 5: R>=D , statement entered Step 5b: R=10 (R-D) Step 5c: Q=10 (setting Q(i) to 1) Step 2: Set i=0 Step 3: R=100 Step 4: R=100 Step 5: R>=D , statement entered Step 5b: R=0 (R-D) Step 5c: Q=11 (setting Q(i) to 1) end Q=11 and R=0. ## Slow division methods Slow division methods are all based on a standard recurrence equation: $P_{j+1}=R\times P_{j}-q_{n-(j+1)}\times D\,\!$ where: • Pj = the partial remainder of the division • R = the radix • q n − (j + 1) = the digit of the quotient in position n-(j+1), where the digit positions are numbered from least-significant 0 to most significant n − 1 • n = number of digits in the quotient • D = the denominator. ### Restoring division Restoring division operates on fixed-point fractional numbers and depends on the following assumptions: {{ safesubst:#invoke:Unsubst||date=__DATE__ |$B= {{#invoke:Category handler|main}}{{#invoke:Category handler|main}}[citation needed] }}

• D < N
• 0 < N,D < 1.

The quotient digits q are formed from the digit set {0,1}.

The basic algorithm for binary (radix 2) restoring division is:

 

P := N D := D << n * P and D need twice the word width of N and Q for i = n-1..0 do * for example 31..0 for 32 bits P := 2P - D * trial subtraction from shifted value if P >= 0 then q(i) := 1 * result-bit 1 else q(i) := 0 * result-bit 0 P := P + D * new partial remainder is (restored) shifted value end end where N=Numerator, D=Denominator, n=#bits, P=Partial remainder, q(i)=bit #i of quotient 

The above restoring division algorithm can avoid the restoring step by saving the shifted value 2P before the subtraction in an additional register T (i.e., TP << 1) and copying register T to P when the result of the subtraction 2P − D is negative.

Non-performing restoring division is similar to restoring division except that the value of 2*P[i] is saved, so D does not need to be added back in for the case of TP[i] ≤ 0.

### Non-restoring division

Non-restoring division uses the digit set {−1,1} for the quotient digits instead of {0,1}. The basic algorithm for binary (radix 2) non-restoring division is:

P := N
i := 0
while i < n do
if P[i] >= 0 then
q[n-(i+1)] := 1
P[i+1] := 2*P[i] - D
else
q[n-(i+1)] := -1
P[i+1] := 2*P[i] + D
end if
i := i + 1
end while


Following this algorithm, the quotient is in a non-standard form consisting of digits of −1 and +1. This form needs to be converted to binary to form the final quotient. Example:

### SRT division

Named for its creators (Sweeney, Robertson, and Tocher), SRT division is a popular method for division in many microprocessor implementations. SRT division is similar to non-restoring division, but it uses a lookup table based on the dividend and the divisor to determine each quotient digit. The Intel Pentium processor's infamous floating-point division bug was caused by an incorrectly coded lookup table. Five of the 1066 entries had been mistakenly omitted.

## Fast division methods

### Newton–Raphson division

The steps of Newton–Raphson are:

1. Calculate an estimate for the reciprocal of the divisor ($D$ ): $X_{0}$ .
2. Compute successively more accurate estimates of the reciprocal: $(X_{1},\ldots ,X_{S})$ 3. Compute the quotient by multiplying the dividend by the reciprocal of the divisor: $Q=NX_{S}$ .

In order to apply Newton's method to find the reciprocal of $D$ , it is necessary to find a function $f(X)$ which has a zero at $X=1/D$ . The obvious such function is $f(X)=DX-1$ , but the Newton–Raphson iteration for this is unhelpful since it cannot be computed without already knowing the reciprocal of $D$ . Moreover multiple iterations for refining reciprocal are not possible since higher order derivatives do not exist for $f(X)$ . A function which does work is $f(X)=1/X-D$ , for which the Newton–Raphson iteration gives

$X_{i+1}=X_{i}-{f(X_{i}) \over f'(X_{i})}=X_{i}-{1/X_{i}-D \over -1/X_{i}^{2}}=X_{i}+X_{i}(1-DX_{i})=X_{i}(2-DX_{i})$ which can be calculated from $X_{i}$ using only multiplication and subtraction, or using two fused multiply–adds.

From a computation point of view the expressions $X_{i+1}=X_{i}+X_{i}(1-DX_{i})$ and $X_{i+1}=X_{i}(2-DX_{i})$ are not equivalent. To obtain a result with a precision of n bits while making use of the second expression one must compute the product between $X_{i}$ and $(2-DX_{i})$ with double the required precision (2n bits). In contrast the product between $X_{i}$ and $(1-DX_{i})$ need only be computed with a precision of n bits.

$X_{i}={1 \over D}(1+\epsilon _{i})\,$ $\epsilon _{i+1}=-{\epsilon _{i}}^{2}\,.$ Apply a bit-shift to the divisor D to scale it so that 0.5 ≤ D ≤ 1 . The same bit-shift should be applied to the numerator N so that the quotient does not change. Then one could use a linear approximation in the form

$X_{0}=T_{1}+T_{2}D\approx {\frac {1}{D}}\,$ to initialize Newton–Raphson. To minimize the maximum of the absolute value of the error of this approximation on interval $[0.5,1]$ one should use

$X_{0}={48 \over 17}-{32 \over 17}D\,.$ {{ safesubst:#invoke:Unsubst||date=__DATE__ |\$B=

{{#invoke:Category handler|main}}{{#invoke:Category handler|main}}[citation needed] }}

Using this approximation, the error of the initial value is less than

$\vert \epsilon _{0}\vert \leq {1 \over 17}\approx 0.059\,.$ Since for this method the convergence is exactly quadratic, it follows that

$S=\left\lceil \log _{2}{\frac {P+1}{\log _{2}17}}\right\rceil \,$ steps is enough to calculate the value up to $P\,$ binary places. This evaluates to 3 for IEEE single precision and 4 for both double precision and double extended formats.

#### Pseudocode

The following computes the quotient of N and D with a precision of P binary places:  

Express D as M × 2e where 1 ≤ M < 2 (standard floating point representation) D' := D / 2e+1 // scale between 0.5 and 1, can be performed with bit shift / exponent subtraction N' := N / 2e+1 X := 48/17 - 32/17 × D' // precompute constants with same precision as D repeat $\left\lceil \log _{2}{\frac {P+1}{\log _{2}17}}\right\rceil \,$ times // can be precomputed based on fixed P X := X + X × (1 - D' × X) end return N' × X 

For example, for a double-precision floating-point division, this method uses 10 multiplies, 9 adds, and 2 shifts.

### Goldschmidt division

Goldschmidt (after Robert Elliott Goldschmidt) division uses an iterative process to repeatedly multiply both the dividend and divisor by a common factor Fi to converge the divisor, D, to 1 as the dividend, N, converges to the quotient Q:

$Q={\frac {N}{D}}{\frac {F_{1}}{F_{1}}}{\frac {F_{2}}{F_{2}}}{\frac {F_{\ldots }}{F_{\ldots }}}.$ The steps for Goldschmidt division are:

1. Generate an estimate for the multiplication factor Fi .
2. Multiply the dividend and divisor by Fi .
3. If the divisor is sufficiently close to 1, return the dividend, otherwise, loop to step 1.

Assuming N/D has been scaled so that 0 < D < 1, each Fi is based on D:

$F_{i+1}=2-D_{i}.$ Multiplying the dividend and divisor by the factor yields:

${\frac {N_{i+1}}{D_{i+1}}}={\frac {N_{i}}{D_{i}}}{\frac {F_{i+1}}{F_{i+1}}}.$ The Goldschmidt method is used in AMD Athlon CPUs and later models.

### Binomial theorem

The Goldschmidt method can be used with factors that allow simplifications by the binomial theorem. Assuming N/D has been scaled by a power of two such that $D\in ({\tfrac {1}{2}},1]$ . We choose $D=1-x$ and $F_{i}=1+x^{2^{i}}$ . This yields

${\frac {N}{1-x}}={\frac {N\cdot (1+x)}{1-x^{2}}}={\frac {N\cdot (1+x)\cdot (1+x^{2})}{1-x^{4}}}=...=Q'={\frac {N'=N\cdot (1+x)\cdot (1+x^{2})\cdot \cdot \cdot (1+x^{2^{(n-1)}})}{D'=1-x^{2^{n}}\approx 1}}$ $\epsilon _{n}={\frac {Q'-N'}{Q'}}=x^{2^{n}}$ This algorithm is referred to as the IBM method in.

## Large integer methods

Methods designed for hardware implementation generally do not scale to integers with thousands or millions of decimal digits; these frequently occur, for example, in modular reductions in cryptography. For these large integers, more efficient division algorithms transform the problem to use a small number of multiplications, which can then be done using an asymptotically efficient multiplication algorithm such as the Karatsuba algorithm, Toom–Cook multiplication or the Schönhage–Strassen algorithm. It results that the computational complexity of the division is of the same order (up a multiplicative constant) as that of the multiplication. Examples include reduction to multiplication by Newton's method as described above as well as the slightly faster Barrett reduction algorithm. Newton's method is particularly efficient in scenarios where one must divide by the same divisor many times, since after the initial Newton inversion only one (truncated) multiplication is needed for each division.

## Division by a constant

The division by a constant D is equivalent to the multiplication by its reciprocal. Since the denominator is constant, so is its reciprocal (1/D). Thus it is possible to compute the value of (1/D) once at compile time, and at run time perform the multiplication N·(1/D) rather than the division N/D. In floating point arithmetic the use of (1/D) presents little problem, but in integer arithmetic the reciprocal will always evaluate to zero (assuming |D| > 1).

It is not necessary to use specifically (1/D); any value (X/Y) that reduces to (1/D) may be used. For example, for division by 3, the factors 1/3, 2/6, 3/9, or 194/582 could be used. Consequently, if Y were a power of two so the division step reduces to a fast right bit shift. The effect of calculating N/D as (N·X)/Y replaces a division with a multiply and a shift. Note that the parentheses are important, as N·(X/Y) will evaluate to zero.

However, unless D itself is a power of two, there is no X and Y that satisfies the conditions above. Fortunately, it is not necessary for (X/Y) to be exactly equal to 1/D, but only that it is "close enough" so that the error introduced by the approximation is in the bits that are discarded by the shift operation.

As a concrete example, for 32 bit unsigned integers, division by 3 can be replaced with a multiply by 2863311531 / 233, a multiplication by 2863311531 followed by a 33 right bit shift. This value is equal to 1/2.999999999650754.

In some cases, division by a constant can be accomplished in even less time by converting the "multiply by a constant" into a series of shifts and adds or subtracts. Of particular interest is division by 10, for which the exact quotient is obtained, with remainder if required.

## Rounding error

Template:Expand section Round-off error can be introduced by division operations due to limited precision.