10-10 Incorporation into a Compiler (Unsigned)

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.)

Figure 10-2 Computing the magic number for unsigned division.

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. 
   int p; 
   unsigned nc, delta, q1, r1, q2, r2; 
   struct mu magu; 
 
   magu.a = 0;             // Initialize "add" indicator. 
   nc = -1 - (-d)%d;       // Unsigned arithmetic here. 
   p = 31;                 // Init. p. 
   q1 = 0x80000000/nc;     // Init. q1 = 2**p/nc. 
   r1 = 0x80000000 - q1*nc;// Init. r1 = rem(2**p, nc). 
   q2 = 0x7FFFFFFF/d;    ?// Init. q2 = (2**p - 1)/d. 
   r2 = 0x7FFFFFFF - q2*d; // Init. r2 = rem(2**p - 1, d). 
   do {
    ?p = p + 1; 
      if (r1 >= nc - r1) {
         q1 = 2*q1 + 1;          ?// Update q1. 
        ?/span>r1 = 2*r1 - nc;}        ?// Update r1. 
      else {
         q1 = 2*q1; 
         r1 = 2*r1;} 
      if (r2 + 1 >= d - r2) {
         if (q2 >= 0x7FFFFFFF) magu.a = 1; 
         q2 = 2*q2 + 1;          ?// Update q2. 
         r2 = 2*r2 + 1 - d;}       // Update r2. 
      else {
         if (q2 >= 0x80000000) magu.a = 1; 
        ?/span>q2 = 2*q2; 
         r2 = 2*r2 + 1;} 
      delta = d - 1 - r2; 
   } while (p < 64 && 
          ?q1 < delta || (q1 == delta && r1 == 0))); 
 
 magu.M = q2 + 1;  ?    // Magic number 
   magu.s = p - 32;      ?// and shift amount to return 
   return 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³², or of q2 to 2³² - 1, because of overflow. Instead, the program tests q2 before overflow can occur. If q2 ever gets as large as 2³² - 1, so that M will be greater than or equal to 2³², then magu.a is set equal to 1. If q2 stays below 2³² - 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.