Division algorithm

Long division

Long division is the standard algorithm used for pen-and-paper division of multi-digit 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 the digit level) at each stage; the multiples then become the digits of the quotient, and the final difference is then the remainder.

When used with a binary radix, this method forms the basis for the (unsigned) integer division with remainder algorithm below. Short division is an abbreviated form of long division suitable for one-digit divisors. Chunking – also known as the partial quotients method or the hangman method – is a less-efficient form of long division which may be easier to understand. By allowing one to subtract more multiples than what one currently has at each stage, a more freeform variant of long division can be developed as well.

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

<span id="mwAc8"></span><span class="kr" id="mwAdA">if</span><span class="w" id="mwAdE"> </span><span class="nv" id="mwAdI">D</span><span class="w" id="mwAdM"> </span><span class="o" id="mwAdQ">=</span><span class="w" id="mwAdU"> </span><span class="mi" id="mwAdY">0</span><span class="w" id="mwAdc"> </span><span class="kr" id="mwAdg">then</span><span class="w" id="mwAdk"> </span><span class="nb" id="mwAdo">error</span><span class="p" id="mwAds">(</span><span class="nv" id="mwAdw">DivisionByZeroException</span><span class="p" id="mwAd0">)</span><span class="w" id="mwAd4"> </span><span class="kr" id="mwAd8">end</span> <span class="nv" id="mwAeA">Q</span><span class="w" id="mwAeE"> </span><span class="p" id="mwAeI">:</span><span class="err" id="mwAeM">=</span><span class="w" id="mwAeQ"> </span><span class="err" id="mwAeU">0</span><span class="w" id="mwAeY"> </span><span class="c1" id="mwAec">-- Initialize quotient and remainder to zero</span> <span class="py" id="mwAeg">R</span><span class="w" id="mwAek"> </span><span class="p" id="mwAeo">:</span><span class="err" id="mwAes">=</span><span class="w" id="mwAew"> </span><span class="err" id="mwAe0">0</span><span class="w" id="mwAe4"> </span> <span class="py" id="mwAe8">for</span><span class="w" id="mwAfA"> </span><span class="nv" id="mwAfE">i</span><span class="w" id="mwAfI"> </span><span class="p" id="mwAfM">:</span><span class="err" id="mwAfQ">=</span><span class="w" id="mwAfU"> </span><span class="py" id="mwAfY">n</span><span class="w" id="mwAfc"> </span><span class="err" id="mwAfg">−</span><span class="w" id="mwAfk"> </span><span class="mi" id="mwAfo">1</span><span class="w" id="mwAfs"> </span><span class="o" id="mwAfw">..</span><span class="w" id="mwAf0"> </span><span class="mi" id="mwAf4">0</span><span class="w" id="mwAf8"> </span><span class="kr" id="mwAgA">do</span><span class="w" id="mwAgE"> </span><span class="c1" id="mwAgI">-- Where n is number of bits in N</span> <span class="w" id="mwAgM"> </span><span class="nv" id="mwAgQ">R</span><span class="w" id="mwAgU"> </span><span class="p" id="mwAgY">:</span><span class="err" id="mwAgc">=</span><span class="w" id="mwAgg"> </span><span class="py" id="mwAgk">R</span><span class="w" id="mwAgo"> </span><span class="o" id="mwAgs"><<</span><span class="w" id="mwAgw"> </span><span class="mi" id="mwAg0">1</span><span class="w" id="mwAg4"> </span><span class="c1" id="mwAg8">-- Left-shift R by 1 bit</span> <span class="w" id="mwAhA"> </span><span class="nf" id="mwAhE">R</span><span class="p" id="mwAhI">(</span><span class="mi" id="mwAhM">0</span><span class="p" id="mwAhQ">)</span><span class="w" id="mwAhU"> </span><span class="p" id="mwAhY">:</span><span class="o" id="mwAhc">=</span><span class="w" id="mwAhg"> </span><span class="nf" id="mwAhk">N</span><span class="p" id="mwAho">(</span><span class="nv" id="mwAhs">i</span><span class="p" id="mwAhw">)</span><span class="w" id="mwAh0"> </span><span class="c1" id="mwAh4">-- Set the least-significant bit of R equal to bit i of the numerator</span> <span class="w" id="mwAh8"> </span><span class="kr" id="mwAiA">if</span><span class="w" id="mwAiE"> </span><span class="nv" id="mwAiI">R</span><span class="w" id="mwAiM"> </span><span class="err" id="mwAiQ">≥</span><span class="w" id="mwAiU"> </span><span class="nv" id="mwAiY">D</span><span class="w" id="mwAic"> </span><span class="kr" id="mwAig">then</span> <span class="w" id="mwAik"> </span><span class="nv" id="mwAio">R</span><span class="w" id="mwAis"> </span><span class="p" id="mwAiw">:</span><span class="err" id="mwAi0">=</span><span class="w" id="mwAi4"> </span><span class="py" id="mwAi8">R</span><span class="w" id="mwAjA"> </span><span class="err" id="mwAjE">−</span><span class="w" id="mwAjI"> </span><span class="nv" id="mwAjM">D</span> <span class="w" id="mwAjQ"> </span><span class="nf" id="mwAjU">Q</span><span class="p" id="mwAjY">(</span><span class="nv" id="mwAjc">i</span><span class="p" id="mwAjg">)</span><span class="w" id="mwAjk"> </span><span class="p" id="mwAjo">:</span><span class="o" id="mwAjs">=</span><span class="w" id="mwAjw"> </span><span class="mi" id="mwAj0">1</span> <span class="w" id="mwAj4"> </span><span class="kr" id="mwAj8">end</span> <span class="kr" id="mwAkA">end</span>

Example

If we take N=1100₂ (12₁₀) and D=100₂ (4₁₀)

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₂ (3₁₀) and R=0.

Slow division methods

Slow division methods are all based on a standard recurrence equation

${\displaystyle R_{j+1}=B\times R_j-q_{n-(j+1)}\times D,}$

where:

• R_j is the j-th partial remainder of the division (R(0) := N(i) step is included)
• 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 least-significant 0 to most significant n−1
• n is number of digits in the quotient
• D is the divisor

Restoring division

Restoring division operates on fixed-point fractional numbers and depends on the assumption 0 < D < N.

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

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

<span id="mwAqI"></span><span class="nv" id="mwAqM">R</span><span class="w" id="mwAqQ"> </span><span class="p" id="mwAqU">:</span><span class="err" id="mwAqY">=</span><span class="w" id="mwAqc"> </span><span class="py" id="mwAqg">N</span> <span class="nv" id="mwAqk">D</span><span class="w" id="mwAqo"> </span><span class="p" id="mwAqs">:</span><span class="err" id="mwAqw">=</span><span class="w" id="mwAq0"> </span><span class="py" id="mwAq4">D</span><span class="w" id="mwAq8"> </span><span class="o" id="mwArA"><<</span><span class="w" id="mwArE"> </span><span class="nv" id="mwArI">n</span><span class="w" id="mwArM"> </span><span class="c1" id="mwArQ">-- R and D need twice the word width of N and Q</span> <span class="kr" id="mwArU">for</span><span class="w" id="mwArY"> </span><span class="nv" id="mwArc">i</span><span class="w" id="mwArg"> </span><span class="p" id="mwArk">:</span><span class="err" id="mwAro">=</span><span class="w" id="mwArs"> </span><span class="py" id="mwArw">n</span><span class="w" id="mwAr0"> </span><span class="err" id="mwAr4">−</span><span class="w" id="mwAr8"> </span><span class="mi" id="mwAsA">1</span><span class="w" id="mwAsE"> </span><span class="o" id="mwAsI">..</span><span class="w" id="mwAsM"> </span><span class="mi" id="mwAsQ">0</span><span class="w" id="mwAsU"> </span><span class="kr" id="mwAsY">do</span><span class="w" id="mwAsc"> </span><span class="c1" id="mwAsg">-- For example 31..0 for 32 bits</span> <span class="w" id="mwAsk"> </span><span class="nv" id="mwAso">R</span><span class="w" id="mwAss"> </span><span class="p" id="mwAsw">:</span><span class="err" id="mwAs0">=</span><span class="w" id="mwAs4"> </span><span class="err" id="mwAs8">2</span><span class="w" id="mwAtA"> </span><span class="err" id="mwAtE">*</span><span class="w" id="mwAtI"> </span><span class="py" id="mwAtM">R</span><span class="w" id="mwAtQ"> </span><span class="err" id="mwAtU">−</span><span class="w" id="mwAtY"> </span><span class="nv" id="mwAtc">D</span><span class="w" id="mwAtg"> </span><span class="c1" id="mwAtk">-- Trial subtraction from shifted value (multiplication by 2 is a shift in binary representation)</span> <span class="w" id="mwAto"> </span><span class="kr" id="mwAts">if</span><span class="w" id="mwAtw"> </span><span class="nv" id="mwAt0">R</span><span class="w" id="mwAt4"> </span><span class="o" id="mwAt8">>=</span><span class="w" id="mwAuA"> </span><span class="mi" id="mwAuE">0</span><span class="w" id="mwAuI"> </span><span class="kr" id="mwAuM">then</span> <span class="w" id="mwAuQ"> </span><span class="nf" id="mwAuU">q</span><span class="p" id="mwAuY">(</span><span class="nv" id="mwAuc">i</span><span class="p" id="mwAug">)</span><span class="w" id="mwAuk"> </span><span class="p" id="mwAuo">:</span><span class="o" id="mwAus">=</span><span class="w" id="mwAuw"> </span><span class="mi" id="mwAu0">1</span><span class="w" id="mwAu4"> </span><span class="c1" id="mwAu8">-- Result-bit 1</span> <span class="w" id="mwAvA"> </span><span class="kr" id="mwAvE">else</span> <span class="w" id="mwAvI"> </span><span class="nf" id="mwAvM">q</span><span class="p" id="mwAvQ">(</span><span class="nv" id="mwAvU">i</span><span class="p" id="mwAvY">)</span><span class="w" id="mwAvc"> </span><span class="p" id="mwAvg">:</span><span class="o" id="mwAvk">=</span><span class="w" id="mwAvo"> </span><span class="mi" id="mwAvs">0</span><span class="w" id="mwAvw"> </span><span class="c1" id="mwAv0">-- Result-bit 0</span> <span class="w" id="mwAv4"> </span><span class="nv" id="mwAv8">R</span><span class="w" id="mwAwA"> </span><span class="p" id="mwAwE">:</span><span class="err" id="mwAwI">=</span><span class="w" id="mwAwM"> </span><span class="py" id="mwAwQ">R</span><span class="w" id="mwAwU"> </span><span class="o" id="mwAwY">+</span><span class="w" id="mwAwc"> </span><span class="nv" id="mwAwg">D</span><span class="w" id="mwAwk"> </span><span class="c1" id="mwAwo">-- New partial remainder is (restored) shifted value</span> <span class="w" id="mwAws"> </span><span class="kr" id="mwAww">end</span> <span class="kr" id="mwAw0">end</span> <span class="c1" id="mwAw4">-- Where: N = numerator, D = denominator, n = #bits, R = partial remainder, q(i) = bit #i of quotient</span>

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

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:

<span id="mwAyE"></span><span class="c1" id="mwAyI">-- Inputs: N (Numerator), D (Denominator)</span> <span class="c1" id="mwAyM">-- n = number of bits</span> <span class="c1" id="mwAyQ">-- R and D are usually stored in registers of width 2n or similar to handle shifts</span> <span class="nv" id="mwAyU">R</span><span class="w" id="mwAyY"> </span><span class="p" id="mwAyc">:</span><span class="err" id="mwAyg">=</span><span class="w" id="mwAyk"> </span><span class="py" id="mwAyo">N</span><span class="w" id="mwAys"> </span><span class="c1" id="mwAyw">-- Initialize Remainder</span> <span class="kr" id="mwAy0">for</span><span class="w" id="mwAy4"> </span><span class="nv" id="mwAy8">i</span><span class="w" id="mwAzA"> </span><span class="o" id="mwAzE">=</span><span class="w" id="mwAzI"> </span><span class="nv" id="mwAzM">n</span><span class="w" id="mwAzQ"> </span><span class="err" id="mwAzU">−</span><span class="w" id="mwAzY"> </span><span class="mi" id="mwAzc">1</span><span class="w" id="mwAzg"> </span><span class="o" id="mwAzk">..</span><span class="w" id="mwAzo"> </span><span class="mi" id="mwAzs">0</span><span class="w" id="mwAzw"> </span><span class="kr" id="mwAz0">do</span><span class="w" id="mwAz4"> </span><span class="c1" id="mwAz8">-- for example 31..0 for 32 bits</span> <span class="w" id="mwA0A"> </span><span class="c1" id="mwA0E">-- Shift Remainder Left (Algebraically: 2 * R)</span> <span class="w" id="mwA0I"> </span><span class="kr" id="mwA0M">if</span><span class="w" id="mwA0Q"> </span><span class="nv" id="mwA0U">R</span><span class="w" id="mwA0Y"> </span><span class="o" id="mwA0c">>=</span><span class="w" id="mwA0g"> </span><span class="mi" id="mwA0k">0</span><span class="w" id="mwA0o"> </span><span class="kr" id="mwA0s">then</span> <span class="w" id="mwA0w"> </span><span class="nv" id="mwA00">R</span><span class="w" id="mwA04"> </span><span class="p" id="mwA08">:</span><span class="err" id="mwA1A">=</span><span class="w" id="mwA1E"> </span><span class="err" id="mwA1I">2</span><span class="w" id="mwA1M"> </span><span class="err" id="mwA1Q">*</span><span class="w" id="mwA1U"> </span><span class="py" id="mwA1Y">R</span><span class="w" id="mwA1c"> </span><span class="o" id="mwA1g">-</span><span class="w" id="mwA1k"> </span><span class="nv" id="mwA1o">D</span><span class="p" id="mwA1s">;</span><span class="w" id="mwA1w"> </span><span class="c1" id="mwA10">-- Subtract D</span> <span class="w" id="mwA14"> </span><span class="nf" id="mwA18">q</span><span class="p" id="mwA2A">(</span><span class="nv" id="mwA2E">i</span><span class="p" id="mwA2I">)</span><span class="w" id="mwA2M"> </span><span class="p" id="mwA2Q">:</span><span class="o" id="mwA2U">=</span><span class="w" id="mwA2Y"> </span><span class="mi" id="mwA2c">1</span><span class="p" id="mwA2g">;</span><span class="w" id="mwA2k"> </span><span class="c1" id="mwA2o">-- Record quotient bit as 1</span> <span class="w" id="mwA2s"> </span><span class="kr" id="mwA2w">else</span> <span class="w" id="mwA20"> </span><span class="nv" id="mwA24">R</span><span class="w" id="mwA28"> </span><span class="p" id="mwA3A">:</span><span class="err" id="mwA3E">=</span><span class="w" id="mwA3I"> </span><span class="err" id="mwA3M">2</span><span class="w" id="mwA3Q"> </span><span class="err" id="mwA3U">*</span><span class="w" id="mwA3Y"> </span><span class="py" id="mwA3c">R</span><span class="w" id="mwA3g"> </span><span class="o" id="mwA3k">+</span><span class="w" id="mwA3o"> </span><span class="nv" id="mwA3s">D</span><span class="p" id="mwA3w">;</span><span class="w" id="mwA30"> </span><span class="c1" id="mwA34">-- Add D (restore)</span> <span class="w" id="mwA38"> </span><span class="nf" id="mwA4A">q</span><span class="p" id="mwA4E">(</span><span class="nv" id="mwA4I">i</span><span class="p" id="mwA4M">)</span><span class="w" id="mwA4Q"> </span><span class="p" id="mwA4U">:</span><span class="o" id="mwA4Y">=</span><span class="w" id="mwA4c"> </span><span class="o" id="mwA4g">-</span><span class="mi" id="mwA4k">1</span><span class="p" id="mwA4o">;</span><span class="w" id="mwA4s"> </span><span class="c1" id="mwA4w">-- Record quotient bit as -1</span> <span class="w" id="mwA40"> </span><span class="kr" id="mwA44">end</span><span class="w" id="mwA48"> </span><span class="kr" id="mwA5A">if</span> <span class="kr" id="mwA5E">end</span><span class="w" id="mwA5I"> </span><span class="kr" id="mwA5M">for</span>

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:

If the −1 digits of ${\displaystyle Q}$ are stored as zeros (0) as is common, then ${\displaystyle P}$ is ${\displaystyle Q}$ and computing ${\displaystyle M}$ is trivial: perform a ones' complement (bit by bit complement) on the original ${\displaystyle Q}$.

<span id="mwA8M"></span><span class="nv" id="mwA8Q">Q</span><span class="w" id="mwA8U"> </span><span class="p" id="mwA8Y">:</span><span class="err" id="mwA8c">=</span><span class="w" id="mwA8g"> </span><span class="py" id="mwA8k">Q</span><span class="w" id="mwA8o"> </span><span class="err" id="mwA8s">−</span><span class="w" id="mwA8w"> </span><span class="nv" id="mwA80">bit</span><span class="p" id="mwA84">.</span><span class="nf" id="mwA88">bnot</span><span class="p" id="mwA9A">(</span><span class="nv" id="mwA9E">Q</span><span class="p" id="mwA9I">)</span><span class="w" id="mwA9M"> </span><span class="c1" id="mwA9Q">-- Appropriate if −1 digits in Q are represented as zeros as is common.</span>

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 non-standard form to standard form:

<span id="mwA9s"></span><span class="kr" id="mwA9w">if</span><span class="w" id="mwA90"> </span><span class="nv" id="mwA94">R</span><span class="w" id="mwA98"> </span><span class="o" id="mwA-A"><</span><span class="w" id="mwA-E"> </span><span class="mi" id="mwA-I">0</span><span class="w" id="mwA-M"> </span><span class="kr" id="mwA-Q">then</span> <span class="w" id="mwA-U"> </span><span class="nv" id="mwA-Y">Q</span><span class="w" id="mwA-c"> </span><span class="p" id="mwA-g">:</span><span class="err" id="mwA-k">=</span><span class="w" id="mwA-o"> </span><span class="py" id="mwA-s">Q</span><span class="w" id="mwA-w"> </span><span class="err" id="mwA-0">−</span><span class="w" id="mwA-4"> </span><span class="mi" id="mwA-8">1</span> <span class="w" id="mwA_A"> </span><span class="nv" id="mwA_E">R</span><span class="w" id="mwA_I"> </span><span class="p" id="mwA_M">:</span><span class="err" id="mwA_Q">=</span><span class="w" id="mwA_U"> </span><span class="py" id="mwA_Y">R</span><span class="w" id="mwA_c"> </span><span class="o" id="mwA_g">+</span><span class="w" id="mwA_k"> </span><span class="nv" id="mwA_o">D</span><span class="w" id="mwA_s"> </span><span class="c1" id="mwA_w">-- Needed only if the remainder is of interest.</span> <span class="kr" id="mwA_0">end</span><span class="w" id="mwA_4"> </span><span class="kr" id="mwA_8">if</span>

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

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 0 × 4 + 2 = 1 × 4 − 2.) This tolerance allows quotient digits to be selected using only a few most-significant bits of the dividend and divisor, rather than requiring a full-width subtraction. This simplification in turn allows a radix higher than 2 to be used.

Like non-restoring division, the final steps are a final full-width subtraction to resolve the last quotient bit, and conversion of the quotient to standard binary form.

The original Intel Pentium processor's infamous floating-point division bug was caused by an incorrectly coded lookup table. Pentium processors used a 2048-cell table, of which 1066 cells were to be populated, and values from five cells were mistakenly omitted.

Fast division methods

Newton–Raphson division

Newton–Raphson uses Newton's method to find the reciprocal of ${\displaystyle D}$ and multiply that reciprocal by ${\displaystyle N}$ to find the final quotient ${\displaystyle Q}$.

The steps of Newton–Raphson division are:

1. Calculate an estimate ${\displaystyle X_0}$ for the reciprocal ${\displaystyle 1/D}$ of the divisor ${\displaystyle D}$.
2. Compute successively more accurate estimates ${\displaystyle X_1,X_2,\dots ,X_S}$ of the reciprocal. This is where one employs the Newton–Raphson method as such.
3. Compute the quotient by multiplying the dividend by the reciprocal of the divisor: ${\displaystyle Q=N{X}_S}$.

In order to apply Newton's method to find the reciprocal of ${\displaystyle D}$, it is necessary to find a function ${\displaystyle f(X)}$ that has a zero at ${\displaystyle X=1/D}$. The obvious such function is ${\displaystyle f(X)=DX-1}$, but the Newton–Raphson iteration for this is unhelpful, since it cannot be computed without already knowing the reciprocal of ${\displaystyle D}$ (moreover it attempts to compute the exact reciprocal in one step, rather than allow for iterative improvements). A function that does work is ${\displaystyle f(X)=(1/X)-D}$, for which the Newton–Raphson iteration gives

$${\displaystyle X_{i+1}=X_i-\frac{f(X_i)}{f^{\prime }(X_i)}=X_i-\frac{1/X_i-D}{-1/X_i^2}=X_i+X_i(1-D{X}_i)=X_i(2-D{X}_i),}$$

which can be calculated from ${\displaystyle X_i}$ using only multiplication and subtraction, or using two fused multiply–adds.

From a computation point of view, the expressions ${\displaystyle X_{i+1}=X_i+X_i(1-D{X}_i)}$ and ${\displaystyle X_{i+1}=X_i(2-D{X}_i)}$ 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 ${\displaystyle X_i}$ and ${\displaystyle (2-D{X}_i)}$ with double the given precision of ${\displaystyle X_i}$(n bits). In contrast, the product between ${\displaystyle X_i}$ and ${\displaystyle (1-D{X}_i)}$ need only be computed with a precision of n bits, because the leading n bits (after the binary point) of ${\displaystyle (1-D{X}_i)}$ are zeros.

If the error is defined as ${\displaystyle {\epsilon }_i=1-D{X}_i}$, then:

${\displaystyle \begin{array}{cc}{\epsilon }_{i+1}& =1-D{X}_{i+1}\\ & =1-D(X_i(2-D{X}_i))\\ & =1-2D{X}_i+D^2{X}_i^2\\ & =(1-D{X}_i{)}^2\\ & ={{\epsilon }_i}^2.\end{array}}$

This squaring of the error at each iteration step – the so-called 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 ${\displaystyle X_0}$ is poorly chosen.

Initial estimate

For the subproblem of choosing an initial estimate ${\displaystyle X_0}$, 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 ${\displaystyle [0.5,1]}$ is:

${\displaystyle X_0=\frac{48}{17}-\frac{32}{17}D.}$

The coefficients of the linear approximation ${\displaystyle T_0+T_{1}D}$ are determined as follows. The absolute value of the error is ${\displaystyle |{\epsilon }_0|=|1-D(T_0+T_{1}D)|}$. The minimum of the maximum absolute value of the error is determined by the Chebyshev equioscillation theorem applied to ${\displaystyle F(D)=1-D(T_0+T_{1}D)}$. The local minimum of ${\displaystyle F(D)}$ occurs when ${\displaystyle F^{\prime }(D)=0}$, which has solution ${\displaystyle D=-T_0/(2T_1)}$. The function at that minimum must be of opposite sign as the function at the endpoints, namely, ${\displaystyle F(1/2)=F(1)=-F(-T_0/(2T_1))}$. The two equations in the two unknowns have a unique solution ${\displaystyle T_0=48/17}$ and ${\displaystyle T_1=-32/17}$, and the maximum error is ${\displaystyle F(1)=1/17}$. Using this approximation, the absolute value of the error of the initial value is less than

${\displaystyle |{\epsilon }_0|\le \frac{1}{17}\approx \mathrm{0.059.}}$

The best quadratic fit to ${\displaystyle 1/D}$ in the interval is

${\displaystyle X:=\frac{140}{33}-\frac{64}{11}D+\frac{256}{99}{D}^2.}$

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 ${\displaystyle {\log }_2(\log 99/\log 17)\approx 0.7}$ Newton–Raphson iterations, at a computational cost of less than one iteration.

It is possible to generate a polynomial fit of degree larger than 2, computing the coefficients using the Remez algorithm. The trade-off 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, from an initial error ${\displaystyle {\epsilon }_0}$, ${\displaystyle S}$ iterations will give an answer accurate to

${\displaystyle P=-2^S{\log }_2{\epsilon }_0-1=2^S{\log }_2(1/{\epsilon }_0)-1}$

binary places. Typical values are:

A quadratic initial estimate plus two iterations is accurate enough for IEEE single precision, but three iterations are marginal for double precision. A linear initial estimate plus four iterations is sufficient for both double 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 ${\displaystyle \lceil {\log }_2\frac{P+1}{{\log }_{2}17}\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.

Cubic iteration

There is an iteration which uses three multiplications to cube the error:

${\displaystyle {\epsilon }_i=1-D{X}_i}$

${\displaystyle Y_i=X_i{\epsilon }_i}$

${\displaystyle X_{i+1}=X_i+Y_i+Y_i{\epsilon }_i.}$

The Y_iε_i term is new.

Expanding out the above, ${\displaystyle X_{i+1}}$ can be written as

${\displaystyle \begin{array}{cc}{X}_{i+1}& =X_i+X_i{\epsilon }_i+X_i{\epsilon }_i^2\\ & =X_i+X_i(1-D{X}_i)+X_i(1-D{X}_i{)}^2\\ & =3X_i-3D{X}_i^2+D^2{X}_i^3,\end{array}}$

with the result that the error term

${\displaystyle \begin{array}{cc}{\epsilon }_{i+1}& =1-D{X}_{i+1}\\ & =1-3D{X}_i+3D^2{X}_i^2-D^3{X}_i^3\\ & =(1-D{X}_i{)}^3\\ & ={\epsilon }_i^3.\end{array}}$

This is 3/2 the computation of the quadratic iteration, but achieves ${\displaystyle \log 3/\log 2\approx 1.585}$ as much convergence, so is slightly more efficient. Put another way, two iterations of this method raise the error to the ninth power at the same computational cost as three quadratic iterations, which only raise the error to the eighth power.

The number of correct bits after ${\displaystyle S}$ iterations is

${\displaystyle P=-3^S{\log }_2{\epsilon }_0-1=3^S{\log }_2(1/{\epsilon }_0)-1}$

binary places. Typical values are:

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.

Goldschmidt division

Goldschmidt division (after Robert Elliott Goldschmidt) 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:

${\displaystyle Q=\frac{N}{D}\frac{F_1}{F_1}\frac{F_2}{F_2}\frac{F_{\dots }}{F_{\dots }}.}$

The steps for Goldschmidt division are:

1. Generate an estimate for the multiplication factor F_i .
2. Multiply the dividend and divisor by F_i .
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 F_i is based on D:

${\displaystyle F_{i+1}=2-D_i.}$

Multiplying the dividend and divisor by the factor yields:

${\displaystyle \frac{N_{i+1}}{D_{i+1}}=\frac{N_i}{D_i}\frac{F_{i+1}}{F_{i+1}}.}$

After a sufficient number k of iterations ${\displaystyle Q=N_k}$.

The Goldschmidt method is used in AMD Athlon CPUs and later models. It is also known as Anderson Earle Goldschmidt Powers (AEGP) algorithm and is implemented by various IBM processors. Although it converges at the same rate as a Newton–Raphson implementation, one advantage of the Goldschmidt method is that the multiplications in the numerator and in the denominator can be done in parallel.

Binomial theorem

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

$${\displaystyle \frac{N}{1-x}=\frac{N\cdot (1+x)}{1-x^2}=\frac{N\cdot (1+x)\cdot (1+x^2)}{1-x^4}=\cdots =Q^{\prime }=\frac{N^{\prime }=N\cdot (1+x)\cdot (1+x^2)\cdot \cdot \cdot (1+x^{2^{(n-1)}})}{D^{\prime }=1-x^{2^n}\approx 1}}$$.

After n steps ${\displaystyle \left(x\in \left[0,\frac{1}{2}\right)\right)}$, the denominator ${\displaystyle 1-x^{2^n}}$ can be rounded to 1 with a relative error

${\displaystyle {\epsilon }_n=\frac{Q^{\prime }-N^{\prime }}{Q^{\prime }}=x^{2^n}}$

which is maximum at ${\displaystyle 2^{-2^n}}$ when ${\displaystyle x=\frac{1}{2}}$, thus providing a minimum precision of ${\displaystyle 2^n}$ binary digits.

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. 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. 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 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 division operation. Y is chosen to be a power of 2 to convert the division into a simple right shift.

As a concrete fixed-point arithmetic example, for 32-bit unsigned integers, division by 3 can be replaced with a multiply by 2863311531/2³³, a widening multiplication by 2863311531 (hexadecimal 0xAAAAAAAB) followed by a 33 right bit shift. The value of 2863311531 is calculated as 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 unsigned integer division where the divisor D is not a power of 2, the following identity converts the division into two x-bit addition/subtraction, one x-bit by x-bit multiplication (where only the upper half of the result is used) and several shifts, after precomputing ${\displaystyle k=x+\lceil {\log }_{2}D\rceil }$ and $a=\lceil 2^k/D\rceil -2^x$:

${\displaystyle \lfloor \frac{N}{D}\rfloor =\lfloor \frac{\lfloor \frac{1}{2}(N-b)\rfloor +b}{2^{k-x-1}}\rfloor ,}$

where ${\displaystyle b=\lfloor Na/2^x\rfloor }$.

Modern compilers commonly perform this integer multiply-and-shift optimization; for a constant only known at run-time, however, the program must implement the optimization itself. libdivide is a library for computing the constants and executing the multiply-shift operation.

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.

Rounding error

When a division operation is performed, the exact quotient ${\displaystyle q}$ and remainder ${\displaystyle r}$ are approximated to fit within the computer's precision limits. The Division Algorithm states:

${\displaystyle [a=bq+r]}$

where ${\displaystyle 0\le r<|b|}$.

In floating-point arithmetic, the quotient ${\displaystyle q}$ is represented as ${\displaystyle \tilde{q}}$ and the remainder ${\displaystyle r}$ as ${\displaystyle \tilde{r}}$, introducing rounding errors ${\displaystyle {ϵ}_q}$ and ${\displaystyle {ϵ}_r}$:

${\displaystyle [\tilde{q}=q+{ϵ}_q][\tilde{r}=r+{ϵ}_r]}$

This rounding causes a small error, which can propagate and accumulate through subsequent calculations. Such errors are particularly pronounced in iterative processes and when subtracting nearly equal values - is told loss of significance. To mitigate these errors, techniques such as the use of guard digits or higher precision (or arbitrary precision) are employed.