There is a difficulty in implementing an
algorithm based on direct evaluation of the expressions used in this proof. Although
p 2W, which is proved above, the case p = 2W can occur
(e.g., for d = 2^{W}
- 2 with W 4). When p = 2W, it is difficult to calculate m, because the dividend in (26) does not fit in a 2W-bit word.

However, it can be implemented by the
"incremental division and remainder" technique of algorithm magic. The algorithm is given in Figure 10-2 for W = 32. It passes back an indicator `a`, which tells
whether or not to generate an `add` instruction. (In the case of signed
division, the caller recognizes this by `M` and `d` having
opposite signs.)

`struct mu {unsigned M;牋牋 // Magic number, `

`牋牋牋牋牋int a;牋牋牋牋牋 // "add" indicator, `

`牋牋牋牋牋int s;};牋牋牋牋 // and shift amount. `

` `

`struct mu magicu(unsigned d) {`

`牋牋牋牋牋牋牋牋牋牋牋牋牋 // Must have 1 <= d <= 2**32-1. `

`牋爄nt p; `

`牋爑nsigned nc, delta, q1, r1, q2, r2; `

`牋爏truct mu magu; `

` `

`牋 magu.a = 0;牋牋牋牋牋牋 // Initialize "add" indicator. `

`牋爊c = -1 - (-d)%d;牋牋牋 // Unsigned arithmetic here. `

`牋爌 = 31;牋牋牋牋牋牋牋牋 // Init. p. `

`牋爍1 = 0x80000000/nc;牋牋 // Init. q1 = 2**p/nc. `

`牋爎1 = 0x80000000 - q1*nc;// Init. r1 = rem(2**p, nc). `

`牋爍2 = 0x7FFFFFFF/d;牋牋?// Init. q2 = (2**p - 1)/d. `

`牋爎2 = 0x7FFFFFFF - q2*d; // Init. r2 = rem(2**p - 1, d). `

`牋燿o {`

`牋牋?p = p + 1; `

`牋牋牋if (r1 >= nc - r1) {`

`牋牋牋牋 q1 = 2*q1 + 1;牋牋牋牋牋?// Update q1. `

`牋牋牋牋?/span>r1 = 2*r1 - nc;}牋牋牋牋?// Update r1. `

牋牋牋else {

牋牋牋牋 q1 = 2*q1;

牋牋牋牋爎1 = 2*r1;}

牋牋牋if (r2 + 1 >= d - r2) {

`牋牋牋牋 if (q2 >= 0x7FFFFFFF) magu.a = 1; `

`牋牋牋牋爍2 = 2*q2 + 1;牋牋牋牋牋?// Update q2. `

`牋牋牋牋爎2 = 2*r2 + 1 - d;}牋牋牋 // Update r2. `

`牋牋牋else {`

`牋牋牋牋 if (q2 >= 0x80000000) magu.a = 1; `

`牋牋牋牋?/span>q2 = 2*q2; `

牋牋牋牋爎2 = 2*r2 + 1;}

牋牋牋delta = d - 1 - r2;

牋爙 while (p < 64 &&

牋牋牋牋牋?q1 < delta || (q1 == delta && r1 == 0)));

`牋 magu.M = q2 + 1;牋?牋牋// Magic number `

`牋爉agu.s = p - 32;牋牋牋?// and shift amount to return `

`牋爎eturn magu;牋牋牋牋牋?// (magu.a was set above). `

`} `

Some key points in understanding this algorithm are as follows:

?span style='font:7.0pt "Times New Roman"'> Unsigned overflow can occur at several places and should be ignored.

?span style='font:7.0pt "Times New Roman"'>
n_{c} = 2^{W} - rem(2^{W}, d) - 1
= (2^{W} - 1) - rem(2^{W} - d, d).

?span style='font:7.0pt "Times New Roman"'>
The quotient and remainder of dividing 2^{p} by n_{c}
cannot be updated in the same way as is done in algorithm magic, because here the quantity `2*r1` can
overflow. Hence the algorithm has the test "`if` `(r1` `>=` `nc` `-` `r1)`,"
whereas "`if` `(2*r1` `>=` `nc)`" would be more natural. A similar remark applies to computing
the quotient and remainder of 2^{p} - 1
divided by d.

?span style='font:7.0pt "Times New Roman"'> 0 d d - 1, so d is representable as a 32-bit unsigned integer.

?span style='font:7.0pt "Times New Roman"'>

?span style='font:7.0pt "Times New Roman"'>
The subtraction of 2^{W}
when the magic number `M` exceeds 2^{W} - 1 is not
explicit in the program; it occurs if the computation of `q2` overflows.

?span style='font:7.0pt "Times New Roman"'>
The "add" indicator, `magu.a`, cannot
be set by a straightforward comparison of `M` to 2^{32}, or
of `q2` to 2^{32} - 1, because of overflow. Instead, the program
tests `q2` before overflow can occur. If `q2` ever gets as large as
2^{32} - 1, so that `M` will be greater than or equal to 2^{32}, then `magu.a` is set
equal to 1. If `q2` stays below 2^{32} - 1, then `magu.a` is left at its
initial value of 0.

?span style='font:7.0pt "Times New Roman"'>
Inequality (27) is equivalent to 2^{p}/n_{c}
> d.

?span style='font:7.0pt "Times New Roman"'>
The loop test needs the condition "`p` `<` `64`"
because without it, overflow of `q1` would cause the program to loop too
many times, giving incorrect results.