# Mostly portable 128 by 64 bit division

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;
}

• 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. Mar 27, 2015 at 19:33
• 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. Mar 27, 2015 at 21:11
• 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? Dec 2, 2015 at 1:44
• @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. Dec 8, 2015 at 1:04

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;
}

• 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. Jun 3, 2020 at 4:59
• 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. Jun 6, 2020 at 16:02

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

• 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. May 26, 2016 at 16:39