10
\$\begingroup\$

I wrote this out of curiosity. It is based on the theory behind Knuth's Algorithm D and is intended to emulate the behavior of the x86 div instruction (though the result is truncated instead of raising an exception for overflow).

I am primarily interested in correctness (it seems to work correctly), and have not yet compared performance with anything simpler.

Two's complement hardware is assumed.

uint64_t div(uint64_t a_lo, uint64_t a_hi, uint64_t b, uint64_t &r)
{
    uint64_t p_lo;
    uint64_t p_hi;
    uint64_t q = 0;

    auto r_hi = a_hi;
    auto r_lo = a_lo;

    int s = 0;
    if(0 == (b >> 63)){

        // Normalize so quotient estimates are
        // no more than 2 in error.

        // Note: If any bits get shifted out of
        // r_hi at this point, the result would
        // overflow.

        s = 63 - bsr(b);
        const auto t = 64 - s;

        b <<= s;
        r_hi = (r_hi << s)|(r_lo >> t);
        r_lo <<= s;
    }

    const auto b_hi = b >> 32;

    /*
    The first full-by-half division places b
    across r_hi and r_lo, making the reduction
    step a little complicated.

    To make this easier, u_hi and u_lo will hold
    a shifted image of the remainder.

    [u_hi||    ][u_lo||    ]
          [r_hi||    ][r_lo||    ]
                [ b  ||    ]
    [p_hi||    ][p_lo||    ]
                  |
                  V
                [q_hi||    ]
    */

    auto q_hat = r_hi / b_hi;

    p_lo = mul(b, q_hat, p_hi);

    const auto u_hi = r_hi >> 32;
    const auto u_lo = (r_hi << 32)|(r_lo >> 32);

    // r -= b*q_hat
    //
    // At most 2 iterations of this...
    while(
        (p_hi > u_hi) ||
        ((p_hi == u_hi) && (p_lo > u_lo))
        )
    {
        if(p_lo < b){
            --p_hi;
        }
        p_lo -= b;
        --q_hat;
    }

    auto w_lo = (p_lo << 32);
    auto w_hi = (p_hi << 32)|(p_lo >> 32);

    if(w_lo > r_lo){
        ++w_hi;
    }

    r_lo -= w_lo;
    r_hi -= w_hi;

    q = q_hat << 32;

    /*
    The lower half of the quotient is easier,
    as b is now aligned with r_lo.

          |r_hi][r_lo||    ]
                [ b  ||    ]
    [p_hi||    ][p_lo||    ]
                        |
                        V
                [q_hi||q_lo]
    */

    q_hat = ((r_hi << 32)|(r_lo >> 32)) / b_hi;

    p_lo = mul(b, q_hat, p_hi);

    // r -= b*q_hat
    //
    // ...and at most 2 iterations of this.
    while(
        (p_hi > r_hi) ||
        ((p_hi == r_hi) && (p_lo > r_lo))
        )
    {
        if(p_lo < b){
            --p_hi;
        }
        p_lo -= b;
        --q_hat;
    }

    r_lo -= p_lo;

    q |= q_hat;

    r = r_lo >> s;

    return q;
}

For convenience, below are the bsr and mul functions used above. The bsr routine comes from this useful site:

int bsr(uint64_t x)
{
    uint64_t y;
    uint64_t r;

    r = (x > 0xFFFFFFFF) << 5; x >>= r;
    y = (x > 0xFFFF    ) << 4; x >>= y; r |= y;
    y = (x > 0xFF      ) << 3; x >>= y; r |= y;
    y = (x > 0xF       ) << 2; x >>= y; r |= y;
    y = (x > 0x3       ) << 1; x >>= y; r |= y;

    return static_cast<int>(r | (x >> 1));
}

uint64_t mul(uint64_t a, uint64_t b, uint64_t &y)
{
    auto a_lo = a & 0x00000000FFFFFFFF;
    auto a_hi = a >> 32;

    auto b_lo = b & 0x00000000FFFFFFFF;
    auto b_hi = b >> 32;

    auto c0 = a_lo * b_lo;
    auto c1 = a_hi * b_lo;
    auto c2 = a_hi * b_hi;

    auto u1 = c1 + (a_lo * b_hi);
    if(u1 < c1){
        c2 += 1LL << 32;
    }

    auto u0 = c0 + (u1 << 32);
    if(u0 < c0){
        ++c2;
    }

    y = c2 + (u1 >> 32);

    return u0;
}
\$\endgroup\$
  • \$\begingroup\$ Is the answer received useful? If so, consider upvoting/selecting it. If not, try leaving a comment to the user or ask for clarification about his points. \$\endgroup\$ – glampert Mar 27 '15 at 19:33
  • \$\begingroup\$ Almost forgot this place existed :| A point has been given! Not sure the Q&A format really applies to code review (not to this one, anyway). If nothing happens in a week, it will be selected as the answer. \$\endgroup\$ – defube Mar 27 '15 at 21:11
  • \$\begingroup\$ I see one problem right away: you're assuming that the quotient will fit in a 64-bit result. What if the dividend is more than 64 bits and the divisor is 1? \$\endgroup\$ – Edward Falk Dec 2 '15 at 1:44
  • \$\begingroup\$ @EdwardFalk The result is undefined. Generally, div is used for the double-by-single word divide required for computing multiple precision quotients. Knuth's explanation is better than what I can fit in a comment, though in short, the arguments are normalized so their "limbs", as provided to the primitive divide, should never produce an undefined result. \$\endgroup\$ – defube Dec 8 '15 at 1:04
3
\$\begingroup\$

The original code seemed too complicated for this so I set out to write a simpler version. Should be fairly portable too, and fast as it is just bit shifts and subtraction.

(It works similarly to how one would accomplish division in assembly language on a processor without hardware division.)

// 128-bit / 64-bit unsigned divide

#include <stdio.h>
#include <stdint.h>

int main(void)
{
  // numerator
  uint64_t a_lo = 1234567;
  uint64_t a_hi = 0;

  // denominator
  uint64_t b = 10;

  // quotient
  uint64_t q = a_lo << 1;

  // remainder
  uint64_t rem = a_hi;

  uint64_t carry = a_lo >> 63;
  uint64_t temp_carry = 0;
  int i;

  for(i = 0; i < 64; i++)
  {
    temp_carry = rem >> 63;
    rem <<= 1;
    rem |= carry;
    carry = temp_carry;

    if(carry == 0)
    {
      if(rem >= b)
      {
        carry = 1;
      }
      else
      {
        temp_carry = q >> 63;
        q <<= 1;
        q |= carry;
        carry = temp_carry;
        continue;
      }
    }

    rem -= b;
    rem -= (1 - carry);
    carry = 1;
    temp_carry = q >> 63;
    q <<= 1;
    q |= carry;
    carry = temp_carry;
  }

  printf("quotient = %llu\n", (long long unsigned int)q);
  printf("remainder = %llu\n", (long long unsigned int)rem);
  return 0;
}
\$\endgroup\$
0
\$\begingroup\$

Only the second solution is correct. The first approach fails for some cases, specially for numerators close to 128 bits limits.

\$\endgroup\$
  • 4
    \$\begingroup\$ Could you expand on this a bit? What cases does the first fail for and why? I feel like there's a good answer in here just waiting to come out. \$\endgroup\$ – RubberDuck May 26 '16 at 16:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.