Division algorithm

Contents
Long division
Non-restoring division
SRT division
Large-integer methods
See also
Notes
References

Contents

Long division

–

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 pseudo-code, all values are treated as unsigned integers.

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

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

= 0 then\n q(i) := 1 -- Result-bit 1\n else\n q(i) := 0 -- Result-bit 0\n R := R + D -- New partial remainder is (restored) shifted value\n end\nend\n\n-- Where: N = numerator, D = denominator, n = #bits, R = partial remainder, q(i) = bit #i of quotient\n"}}'>R:=N D:=D<<n-- R and D need twice the word width of N and Q fori:=n−1..0do-- 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) ifR>=0then q(i):=1-- Result-bit 1 else q(i):=0-- Result-bit 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

Non-performing 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.

Non-restoring division

edit

Non-restoring 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, which potentially cuts down the numbers of operations by up to half and lets it be executed faster. The basic algorithm for binary (radix 2) non-restoring division of non-negative numbers is:Q=1111¯11¯11¯ 1. Form the positive term:P−M:Q are stored as zeros (0) as is common, then M is trivial: perform a ones' complement (bit by bit complement) on the original ifR<0then Q:=Q−1 R:=R+D-- Needed only if the remainder is of interest. endif

The actual remainder is R >> n. (As with restoring division, the low-order 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

edit

SRT division is a popular method for division in many microprocessor implementations. The algorithm is named after D. W. Sweeney of IBM, James E. Robertson of University of Illinois, and K. D. Tocher of Imperial College London. They all developed the algorithm independently at approximately the same time (published in February 1957, September 1958, and January 1958 respectively).

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 most significant difference is that a redundant representation is used for the quotient. For example, when implementing radix-4 SRT division, each quotient digit is chosen from five possibilities: −2, −1, 0, +1, or +2. Because of this, the choice of a quotient digit need not be perfect; later quotient digits can correct for slight errors. (For example, the quotient digit pairs (0, +2) and (1, −2) are equivalent, since N to find the Q</math>."}},"i":0}}]}'>final quotient 1/D of the divisor Q=NXS.

In order to apply Newton's method to find the reciprocal of X=1/D. The obvious such function is f(X)=(1/X)−D, for which the Newton–Raphson iteration gives

Xi+1=Xi+Xi(1−DXi) and (2−DXi) with double the given precision of Xi and εi=1−DXi, then:

X0, it is convenient to apply a bit-shift to the divisor D to scale it so that 0.5 ≤ D ≤ 1. Applying the same bit-shift to the numerator N ensures the quotient does not change. Once within a bounded range, a simple polynomial approximation can be used to find an initial estimate.

The linear approximation with minimum worst-case absolute error on the interval T0+T1D are determined as follows. The absolute value of the error is F(D) occurs when F(1/2)=F(1)=−F(−T0/(2T1)). The two equations in the two unknowns have a unique solution F(1)=1/17. Using this approximation, the absolute value of the error of the initial value is less than

X:=14033−6411D+25699D2.

It is chosen to make the error equal to a re-scaled third order Chebyshev polynomial of the first kind, and gives an absolute value of the error less than or equal to 1/99. This improvement is equivalent to S iterations will give an answer accurate to

1/17 15.35 1/99 25.52 N and '''times'''"}},"i":0}}]}'>repeat Yi=Xiεi

Xi+1=Xi+Xiεi+Xiεi2=Xi+Xi(1−DXi)+Xi(1−DXi)2=3Xi−3DXi2+D2Xi3,

with the result that the error term

S iterations is

1/17 35.79 5.63 177.99

A quadratic initial estimate plus two cubic iterations provides ample precision for an IEEE double-precision result. It is also possible to use a mixture of quadratic and cubic iterations.

Using at least one quadratic iteration ensures that the error is positive, i.e. the reciprocal is underestimated. This can simplify a following rounding step if an exactly-rounded quotient is required.

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.Fi+1=2−Di.

Multiplying the dividend and divisor by the factor yields:

⁠D=1−x and (x∈[0,12)), the denominator εn=Q′−N′Q′=x2n

which is maximum at 2n binary digits.

Large-integer methods

edit

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, as well as the slightly faster Burnikel-Ziegler division, Barrett reduction and Montgomery reduction algorithms.2/3, then rounded up. Likewise, division by 10 can be expressed as a multiplication by 3435973837 (0xCCCCCCCD) followed by division by 2 (or 35 right bit shift). OEIS provides sequences of the constants for multiplication as A346495 and for the right shift as A346496.

For general x-bit addition/subtraction, one k=x+⌈log2D⌉ and q and remainder 0≤r<|b|.

In floating-point arithmetic, the quotient r as ϵr:

Notes

edit

↑ Despite how "little" problem the optimization causes, this reciprocal optimization is still usually hidden behind a "fast math" flag in modern compilers as it is inexact.