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 &amp;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 &gt;&gt; 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 &lt;&lt;= s;
        r_hi = (r_hi &lt;&lt; s)|(r_lo &gt;&gt; t);
        r_lo &lt;&lt;= s;
    }

    const auto b_hi = b &gt;&gt; 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 &gt;&gt; 32;
    const auto u_lo = (r_hi &lt;&lt; 32)|(r_lo &gt;&gt; 32);

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

    auto w_lo = (p_lo &lt;&lt; 32);
    auto w_hi = (p_hi &lt;&lt; 32)|(p_lo &gt;&gt; 32);

    if(w_lo &gt; r_lo){
        ++w_hi;
    }

    r_lo -= w_lo;
    r_hi -= w_hi;

    q = q_hat &lt;&lt; 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 &lt;&lt; 32)|(r_lo &gt;&gt; 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 &gt; r_hi) ||
        ((p_hi == r_hi) &amp;&amp; (p_lo &gt; r_lo))
        )
    {
        if(p_lo &lt; b){
            --p_hi;
        }
        p_lo -= b;
        --q_hat;
    }

    r_lo -= p_lo;

    q |= q_hat;

    r = r_lo &gt;&gt; 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 &gt; 0xFFFFFFFF) &lt;&lt; 5; x &gt;&gt;= r;
    y = (x &gt; 0xFFFF    ) &lt;&lt; 4; x &gt;&gt;= y; r |= y;
    y = (x &gt; 0xFF      ) &lt;&lt; 3; x &gt;&gt;= y; r |= y;
    y = (x &gt; 0xF       ) &lt;&lt; 2; x &gt;&gt;= y; r |= y;
    y = (x &gt; 0x3       ) &lt;&lt; 1; x &gt;&gt;= y; r |= y;

    return static_cast&lt;int&gt;(r | (x &gt;&gt; 1));
}

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

    auto b_lo = b &amp; 0x00000000FFFFFFFF;
    auto b_hi = b &gt;&gt; 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 &lt; c1){
        c2 += 1LL &lt;&lt; 32;
    }

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

    y = c2 + (u1 &gt;&gt; 32);

    return u0;
}
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 &lt;stdio.h&gt;
#include &lt;stdint.h&gt;

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 &lt;&lt; 1;

  // remainder
  uint64_t rem = a_hi;

  uint64_t carry = a_lo &gt;&gt; 63;
  uint64_t temp_carry = 0;
  int i;

  for(i = 0; i &lt; 64; i++)
  {
    temp_carry = rem &gt;&gt; 63;
    rem &lt;&lt;= 1;
    rem |= carry;
    carry = temp_carry;

    if(carry == 0)
    {
      if(rem &gt;= b)
      {
        carry = 1;
      }
      else
      {
        temp_carry = q &gt;&gt; 63;
        q &lt;&lt;= 1;
        q |= carry;
        carry = temp_carry;
        continue;
      }
    }

    rem -= b;
    rem -= (1 - carry);
    carry = 1;
    temp_carry = q &gt;&gt; 63;
    q &lt;&lt;= 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;
}