Hacker's Delight, Second Edition. FIGURE 9-3. Divide long unsigned, using fullword division instruction.

divlu.c

unsigned divlu(unsigned u1, unsigned u0, unsigned v,
               unsigned *r)   {
   const unsigned b = 65536;  // Number base (16 bits).
   unsigned un1, un0,         // Norm. dividend LSD's.
            vn1, vn0,         // Norm. divisor digits.
            q1, q0,           // Quotient digits.
            un32, un21, un10, // Divident digit pairs.
            rhat;             // A remainder.
   int s;                     // Shift amount for norm.

   if (u1 >= v) {             // If overflow, set rem.
      if (r != NULL)          // to an impossible value,
         *r = 0xFFFFFFFF;     // and return the largest
      return 0xFFFFFFFF;}     // possible quotient.

   s = nlz(v);                // 0 <= s <= 31.
   v = v << s;                // Normalize divisor.
   vn1 = v >> 16;             // Break divisor up into
   vn0 = v & 0xFFFF;          // two 16-bit digits.

   un32 = (u1 << s) | (u0 >> 32 - s) & (-s >> 31);
   un10 = u0 << s;            // Shift dividend left.

   un1 = un10 >> 16;          // Compute the first
   un0 = un10 & 0xFFFF;       // quotient digit, q1.

   q1 = un32/vn1;
   rhat = un32 - q1*vn1;
again1:
   if (q1 >= b || q1*vn0 > v*rhat + un1) {
      q1 = q1 - 1;
      rhat = rhat + vn1;
      if (rhat < b) goto again1;}

   un21 = un32*b + un1 - q1*v;   // Multiply and subtract

   q0 = un21/vn1;             // Compute the second
   rhat = un21 - q0*vn1;      // quotient digit, q0.
again2:
   if (q0 >= b || q0*vn0 > b*rhat + un0) {
      q0 = q0 - 1;
      rhat = rhat + vn1;
      if (rhat < b) goto again2;}

   if (r != NULL)             // If remainder is wanted,
      *r = (un21*b + un0 - q0*v) >> s;       // return it.
   return q1*b + q0;
}

// vim: sw=3