Division algorithm
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, nonperforming restoring, nonrestoring, 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.
Discussion will refer to the form , where
 N = Numerator (dividend)
 D = Denominator (divisor)
is the input, and
 Q = Quotient
 R = Remainder
is the output.
Contents
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 error(DivisionByZero) 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
return divide_unsigned(N, D)
end
function divide_unsigned(N, D)
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 outputsensitive algorithm), and can serve as an executable specification.
Long division
Long division is the standard algorithm used for penandpaper 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 onedigit divisors. Chunking (also known as the partial quotients method or the hangman method) is a lessefficient form of long division which may be easier to understand.
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. In the following code, all values are treated as unsigned integers.
if D = 0 then error(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  Leftshift R by 1 bit
R(0) := N(i)  Set the leastsignificant 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_{2} (12_{10}) and D=100_{2} (4_{10})
Step 1: Set R=0 and Q=0
Step 2: Take i=3 (one less than the number of bits in N)
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_{2} (3_{10}) and R=0.
Slow division methods
Slow division methods are all based on a standard recurrence equation
 ,
where:
 R_{j} is the jth partial remainder of the division
 B is the radix (base, usually 2 internally in computers and calculators)
 q _{n − (j + 1)} is the digit of the quotient in position n−(j+1), where the digit positions are numbered from leastsignificant 0 to most significant n−1
 n is number of digits in the quotient
 D is the divisor
Restoring division
Restoring division operates on fixedpoint fractional numbers and depends on the following assumptions:^{[citation needed]}
 D < N
 0 < N, 0 < D
The quotient digits q are formed from the digit set {0,1}.
The basic algorithm for binary (radix 2) restoring division is:
R := N
D := D << n  R and D need twice the word width of N and Q
for i := n − 1 .. 0 do  For example 31..0 for 32 bits
R := 2 * R − D  Trial subtraction from shifted value (multiplication by 2 is a shift in binary representation)
if R ≥ 0 then
q(i) := 1  Resultbit 1
else
q(i) := 0  Resultbit 0
R := R + D  New partial remainder is (restored) shifted value
end
end
 Where: N = Numerator, D = Denominator, n = #bits, R = Partial remainder, q(i) = bit #i of quotient
The above restoring division algorithm can avoid the restoring step by saving the shifted value 2R before the subtraction in an additional register T (i.e., T = R << 1) and copying register T to R when the result of the subtraction 2R − D is negative.
Nonperforming restoring division is similar to restoring division except that the value of 2R is saved, so D does not need to be added back in for the case of R < 0.
Nonrestoring division
Nonrestoring division uses the digit set {−1, 1} for the quotient digits instead of {0, 1}. The algorithm is more complex, but has the advantage when implemented in hardware that there is only one decision and addition/subtraction per quotient bit; there is no restoring step after the subtraction. This lets it be executed faster. The basic algorithm for binary (radix 2) nonrestoring division of nonnegative numbers is:
R := N
D := D << n  R and D need twice the word width of N and Q
for i = n − 1 .. 0 do  for example 31..0 for 32 bits
if R >= 0 then
q[i] := +1
R := 2 * R − D
else
q[i] := −1
R := 2 * R + D
end if
end
 Note: N=Numerator, D=Denominator, n=#bits, R=Partial remainder, q(i)=bit #i of quotient.
Following this algorithm, the quotient is in a nonstandard form consisting of digits of −1 and +1. This form needs to be converted to binary to form the final quotient. Example:
Convert the following quotient to the digit set {0,1}:  
Start:  
1. Form the positive term:  
2. Mask the negative term*:  
3. Subtract: :  
*.( Signed binary notation with One's complement without Two's Complement ) 
If the −1 digits of are stored as zeros (0) as is common, then is and computing is trivial: perform a one's complement (bit by bit complement) on the original .
Q := Q − bit.bnot(Q) * Appropriate if −1 Digits in Q are Represented as zeros as is common.
Finally, quotients computed by this algorithm are always odd, and the remainder in R is in the range −D ≤ R < D. For example, 5 / 2 = 3 R −1. To convert to a positive remainder, do a single restoring step after Q is converted from nonstandard form to standard form:
if R < 0 then
Q := Q − 1
R := R + D  Needed only if the Remainder is of interest.
end if
The actual remainder is R >> n. (As with restoring division, the loworder bits of R are used up at the same rate as bits of the quotient Q are produced, and it is common to use a single shift register for both.)
SRT division
Named for its creators (Sweeney, Robertson, and Tocher), SRT division is a popular method for division in many microprocessor implementations.^{[1]}^{[2]} SRT division is similar to nonrestoring 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 floatingpoint division bug was caused by an incorrectly coded lookup table. Five of the 1066 entries had been mistakenly omitted.^{[3]}^{[4]}
Fast division methods
Newton–Raphson division
Newton–Raphson uses Newton's method to find the reciprocal of and multiply that reciprocal by to find the final quotient .
The steps of Newton–Raphson division are:
 Calculate an estimate for the reciprocal of the divisor .
 Compute successively more accurate estimates of the reciprocal. This is where one employs the Newton–Raphson method as such.
 Compute the quotient by multiplying the dividend by the reciprocal of the divisor: .
In order to apply Newton's method to find the reciprocal of , it is necessary to find a function that has a zero at . The obvious such function is , but the Newton–Raphson iteration for this is unhelpful, since it cannot be computed without already knowing the reciprocal of (moreover it attempts to compute the exact reciprocal in one step, rather than allow for iterative improvements). A function that does work is , for which the Newton–Raphson iteration gives
which can be calculated from using only multiplication and subtraction, or using two fused multiply–adds.
From a computation point of view, the expressions and are not equivalent. To obtain a result with a precision of 2n bits while making use of the second expression, one must compute the product between and with double the given precision of (n bits).^{[citation needed]} In contrast, the product between and need only be computed with a precision of n bits, because the leading n bits (after the binary point) of are zeros.
If the error is defined as , then:
This squaring of the error at each iteration step — the socalled quadratic convergence of Newton–Raphson's method — has the effect that the number of correct digits in the result roughly doubles for every iteration, a property that becomes extremely valuable when the numbers involved have many digits (e.g. in the large integer domain). But it also means that the initial convergence of the method can be comparatively slow, especially if the initial estimate is poorly chosen.
For the subproblem of choosing an initial estimate , it is convenient to apply a bitshift to the divisor D to scale it so that 0.5 ≤ D ≤ 1; by applying the same bitshift to the numerator N, one ensures the quotient does not change. Then one could use a linear approximation in the form
to initialize Newton–Raphson. To minimize the maximum of the absolute value of the error of this approximation on interval , one should use
The coefficients of the linear approximation are determined as follows. The absolute value of the error is . The minimum of the maximum absolute value of the error is determined by the Chebyshev equioscillation theorem applied to . The local minimum of occurs when , which has solution . The function at that minimum must be of opposite sign as the function at the endpoints, namely, . The two equations in the two unknowns have a unique solution and , and the maximum error is . Using this approximation, the absolute value of the error of the initial value is less than
It is possible to generate a polynomial fit of degree larger than 1, computing the coefficients using the Remez algorithm. The tradeoff is that the initial guess requires more computational cycles but hopefully in exchange for fewer iterations of Newton–Raphson.
Since for this method the convergence is exactly quadratic, it follows that
steps are enough to calculate the value up to 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 × 2^{e} where 1 ≤ M < 2 (standard floating point representation)
D' := D / 2^{e+1} // scale between 0.5 and 1, can be performed with bit shift / exponent subtraction
N' := N / 2^{e+1}
X := 48/17 − 32/17 × D' // precompute constants with same precision as D
repeat times // can be precomputed based on fixed P
X := X + X × (1  D' × X)
end
return N' × X
For example, for a doubleprecision floatingpoint division, this method uses 10 multiplies, 9 adds, and 2 shifts.
Variant NewtonRaphson division
The NewtonRaphson division method can be modified to be slightly faster as follows. After shifting N and D so that D is in [0.5, 1.0], initialize with
This is the best quadratic fit to 1/D and gives an absolute value of the error less than or equal to 1/99. It is chosen to make the error equal to a rescaled third order Chebyshev polynomial of the first kind. The coefficients should be precalculated and hardcoded.
Then in the loop, use an iteration which cubes the error.
The Y·E term is new.
If the loop is performed until X agrees with 1/D on its leading P bits, then the number of iterations will be no more than
which is the number of times 99 must be cubed to get to 2^{P+1}. Then
is the quotient to P bits.
Using higher degree polynomials in either the initialization or the iteration results in a degradation of performance because the extra multiplications required would be better spent on doing more iterations.
Goldschmidt division
Goldschmidt (after Robert Elliott Goldschmidt)^{[5]} division uses an iterative process of repeatedly multiplying both the dividend and divisor by a common factor F_{i}, chosen such that the divisor converges to 1. This causes the dividend to converge to the sought quotient Q:
The steps for Goldschmidt division are:
 Generate an estimate for the multiplication factor F_{i} .
 Multiply the dividend and divisor by F_{i} .
 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 F_{i} is based on D:
Multiplying the dividend and divisor by the factor yields:
After a sufficient number k of iterations .
The Goldschmidt method is used in AMD Athlon CPUs and later models.^{[6]}^{[7]}
Binomial theorem
The Goldschmidt method can be used with factors that allow simplifications by the binomial theorem. Assume N/D has been scaled by a power of two such that . We choose and . This yields
 .
After steps , the denominator can be rounded to with a relative error
which is maximum at when , thus providing a minimum precision of binary digits.
Largeinteger 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. The result is that the computational complexity of the division is of the same order (up to a multiplicative constant) as that of the multiplication. Examples include reduction to multiplication by Newton's method as described above,^{[8]} as well as the slightly faster Barrett reduction and Montgomery reduction algorithms.^{[9]}^{[verification needed]} 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 floatingpoint 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 the division step would reduce 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, (N·X)/Y gives exactly the same result as N/D in integer arithmetic even when (X/Y) is not exactly equal to 1/D, but "close enough" that the error introduced by the approximation is in the bits that are discarded by the shift operation.^{[10]}^{[11]}^{[12]}
As a concrete fixedpoint arithmetic example, for 32bit unsigned integers, division by 3 can be replaced with a multiply by 2863311531 ⁄ 2^{33}, a multiplication by 2863311531 (hexadecimal 0xAAAAAAAB) followed by a 33 right bit shift. The value of 2863311531 is calculated as 2^{33} ⁄ 3, then rounded up.
Likewise, division by 10 can be expressed as a multiplication by 3435973837 (0xCCCCCCCD) followed by division by 2^{35} (or 35 right bit shift).
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.^{[13]} Of particular interest is division by 10, for which the exact quotient is obtained, with remainder if required.^{[14]}
Rounding error
This section needs expansion. You can help by adding to it. (September 2012)

Roundoff error can be introduced by division operations due to limited precision.
See also
References
 ^ SRT Division: Architectures, Models, and Implementations; Stanford University; 1998.
 ^ SRT Division Algorithsm as Dynamical Systems; Society for Industrial and Applied Mathematics; 2005.
 ^ "Statistical Analysis of Floating Point Flaw". Intel Corporation. 1994. Retrieved 22 October 2013.
 ^ An Analysis of Division Algorithms and Implementations; Stanford University; 1995.
 ^ Goldschmidt, Robert E. (1964). Applications of Division by Convergence (PDF) (Thesis). M.Sc. dissertation. M.I.T. OCLC 34136725.
 ^ Oberman, Stuart F. (1999). "Floating Point Division and Square Root Algorithms and Implementation in the AMDK7 Microprocessor" (PDF). Proceedings of the IEEE Symposium on Computer Arithmetic: 106–115.
 ^ Soderquist, Peter; Leeser, Miriam (July–August 1997). "Division and Square Root: Choosing the Right Implementation" (PDF). IEEE Micro. 17 (4): 56–66.
 ^ Hasselström, Karl (2003). Fast Division of Large Integers: A Comparison of Algorithms (PDF) (M.Sc. in Computer Science thesis). Royal Institute of Technology. Archived from the original (PDF) on 8 July 2017. Retrieved 20170708.
 ^ Barrett, Paul (1987). "Implementing the Rivest Shamir and Adleman public key encryption algorithm on a standard digital signal processor". Proceedings on Advances in cryptologyCRYPTO '86. London, UK: SpringerVerlag. pp. 311–323. ISBN 0387180478.
 ^ Granlund, Torbjörn; Montgomery, Peter L. (June 1994). "Division by Invariant Integers using Multiplication" (PDF). SIGPLAN Notices. ACM. 29 (6): 61–72. doi:10.1145/773473.178249.
 ^ Möller, Niels; Granlund, Torbjörn (February 2011). "Improved Division by Invariant Integers" (PDF). IEEE Transactions on Computers. IEEE. 60 (2): 165–175. doi:10.1109/TC.2010.143.
 ^ ridiculous_fish. "Labor of Division (Episode III): Faster Unsigned Division by Constants". 2011.
 ^ LaBudde, Robert A.; Golovchenko, Nikolai; Newton, James; and Parker, David; Massmind: "Binary Division by a Constant"
 ^ Vowels, R. A. (1992). "Division by 10". Australian Computer Journal. 24 (3): 81–85.
Further reading
 Warren Jr., Henry S. (2013). Hacker's Delight (2 ed.). Addison Wesley  Pearson Education, Inc. ISBN 9780321842688.
External links
 Computer Arithmetic Algorithms JavaScript Simulator – contains simulators for many different division algorithms
 Doras, Cory (19 October 2011). "Labor of Division (Episode III): Faster Unsigned Division by Constants" (PDF). ridiculous_fish. (Extends division by constants.)
 http://www.dauniv.ac.in/downloads/CArch_PPTs/CompArchCh03L07IntegerDivision.pdf
 http://www.seas.ucla.edu/~ingrid/ee213a/lectures/division_presentV2.pdf
 http://www.m1c4a1.wz.cz/docs/goldschmidt.pdf