11
\$\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;
}
| improve this answer | |
\$\endgroup\$
  • 1
    \$\begingroup\$ The test case for this code has only very small values. I replaced a_hi, a_lo, and b with these values: 0xc3ac542864f652f8, 0xa3ccac1832952398,0x5267b85f8a42fc20. This asks for the division 260094241160506081433645566617850291096 / 5937917354192337952. The program gives quotient 15148923828744002542 and remainder 171728852755111384. These values are wrong. But they would be the right values if a_hi was 0x43ac542864f652f8. Perhaps the highest bit of the numerator is not being handled properly. \$\endgroup\$ – Jeffrey Scofield Jun 3 at 4:59
  • \$\begingroup\$ Thanks for revealing the flaw. I don't have time to sort it out so will let your post stand as a warning for now. \$\endgroup\$ – joe63074 Jun 6 at 16:02
1
\$\begingroup\$

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

| improve this answer | |
\$\endgroup\$
  • 5
    \$\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.