How to divide a 128-bit dividend by a 64-bit divisor, where the dividend's bits are all 1's, and where I only need the 64 LSBs of the quotient?

Question

I need to calculate (2¹²⁸ - 1) / x. The divisor, x, is an unsigned 64-bit number. The dividend is composed of two unsigned 64-bit numbers (high and low), where both numbers are UINT64_MAX. I can only use 64-bit arithmetic and need it to be portable (no use of GNU's __int128, MSCV's _udiv128, assembly, or anything like that). I don't need the high part of the quotient, I only need the lower 64 bits.

How can I do this operation?

Also: x >= 3, x is not a power of 2.

Edit: I created my own solution (answer below). But I welcome any other solution that performs better :)

Temperate

Welcome to SO. You might find reading the site help center useful when it comes to How to Ask, and this question checklist. Code that you've worked on to solve the problem should include a minimal reproducible example, and be included in your question. While homework questions are not off-topic per se, you don't expect somebody to do your homework for you – Kemp

Wikipedia has some binary division algorithms in its article. But simplest might be to treat the problem like the "long division" method taught in schools, with each 64 bit word as a "digit". – Townsfolk

Answer 1

This is what I ended up coding. I'm sure there are much faster alternatives, but at least this is functional.

Based on: https://en.wikipedia.org/wiki/Division_algorithm#Integer_division_(unsigned)_with_remainder. Adapted for for this particular use-case.

// q = (2^128 - 1) / d, where q is the 64 LSBs of the quotient
uint64_t two_pow_128_minus_1_div_d(uint64_t d) {
    uint64_t q = 0, r_hi = 0, r_lo = 0;

    for (int i = 127; i >= 0; --i) {
        r_hi = (r_hi << 1) | (r_lo >> 63);
        r_lo <<= 1;

        r_lo |= 1UL;

        if (r_hi || r_lo >= d) {
            const uint64_t borrow = d > r_lo;
            r_lo -= d;
            r_hi -= borrow;

            if (i < 64)
                q |= 1UL << i;
        }
    }
    return q;
}

Temperate

Answer 2

I am not aware of any optimizations that apply to integer division with a constant dividend. To double check, I tried a test case with an all-ones dividend with Compiler Explorer. Using gcc, icc, and clang, with highest optimization level specified, the generated code showed no optimizations being applied to the division.

It is certainly possible to build high-performance 128-bit division routines, but from personal experience I know that this is quite error prone, and very sophisticated testing is needed to achieve good test coverage including corner cases, as exhaustive test is not possible at this operand size. The effort for design and test easily exceeds what seems reasonable for an answer on Stackoverflow by two decimal orders of magnitude.

An easy way to perform integer division is to use the algorithm we all learned in grade school, only in binary. This makes the decision about the next quotient bit particularly easy: It is 1 when the current partial remainder is greater than, or equal to, the divisor, and 0 otherwise. Using longhand binary division, the only integer operations we need are additions and subtractions.

We can build portable primitives for performing these on operands of any bit length by mimicking the way a processor's machine instructions are used to effect operations on multi-word integers: ADD with carry-out, ADD with carry-in, ADD with carry-in and carry-out; analogous for SUB. In the code below I am using simple C macros for that; certainly more sophisticated approaches are possible.

Since the system I am working on right now does not have support for 128-bit integers, I prototyped and tested this approach for 64-bit integers. The 128-bit version then was an exercise in simple mechanical renaming. On a modern 64-bit processor I would expect this 128-bit division function to execute in roughly 3000 cycles.

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

#define SUBCcc(a,b,cy,t0,t1,t2) \
  (t0=(b)+cy, t1=(a), cy=t0<cy, t2=t1<t0, cy=cy+t2, t1-t0)

#define SUBcc(a,b,cy,t0,t1) \
  (t0=(b), t1=(a), cy=t1<t0, t1-t0)

#define SUBC(a,b,cy,t0,t1) \
  (t0=(b)+cy, t1=(a), t1-t0)

#define ADDCcc(a,b,cy,t0,t1) \
  (t0=(b)+cy, t1=(a), cy=t0<cy, t0=t0+t1, t1=t0<t1, cy=cy+t1, t0=t0)

#define ADDcc(a,b,cy,t0,t1) \
  (t0=(b), t1=(a), t0=t0+t1, cy=t0<t1, t0=t0)

#define ADDC(a,b,cy,t0,t1) \
  (t0=(b)+cy, t1=(a), t0+t1)

typedef struct {
    uint64_t l;
    uint64_t h;
} my_uint128;

my_uint128 bitwise_division_128 (my_uint128 dvnd, my_uint128 dvsr)
{
    my_uint128 quot, rem, tmp;
    uint64_t cy, t0, t1, t2;
    int bits_left = CHAR_BIT * sizeof (my_uint128);
    
    quot.h = dvnd.h;
    quot.l = dvnd.l;
    rem.h = 0;
    rem.l = 0;
    do {
        quot.l = ADDcc  (quot.l, quot.l, cy, t0, t1);
        quot.h = ADDCcc (quot.h, quot.h, cy, t0, t1);
        rem.l  = ADDCcc (rem.l,  rem.l,  cy, t0, t1);
        rem.h  = ADDC   (rem.h,  rem.h,  cy, t0, t1);
        tmp.l  = SUBcc  (rem.l,  dvsr.l, cy, t0, t1);
        tmp.h  = SUBCcc (rem.h,  dvsr.h, cy, t0, t1, t2);
        if (!cy) { // remainder >= divisor
            rem.l = tmp.l;
            rem.h = tmp.h;
            quot.l = quot.l | 1;
        }
        bits_left--;
    } while (bits_left);
    return quot;
}

typedef struct {
    uint32_t l;
    uint32_t h;
} my_uint64;

my_uint64 bitwise_division_64 (my_uint64 dvnd, my_uint64 dvsr)
{
    my_uint64 quot, rem, tmp;
    uint32_t cy, t0, t1, t2;
    int bits_left = CHAR_BIT * sizeof (my_uint64);
    
    quot.h = dvnd.h;
    quot.l = dvnd.l;
    rem.h = 0;
    rem.l = 0;
    do {
        quot.l = ADDcc  (quot.l, quot.l, cy, t0, t1);
        quot.h = ADDCcc (quot.h, quot.h, cy, t0, t1);
        rem.l  = ADDCcc (rem.l,  rem.l,  cy, t0, t1);
        rem.h  = ADDC   (rem.h,  rem.h,  cy, t0, t1);
        tmp.l  = SUBcc  (rem.l,  dvsr.l, cy, t0, t1);
        tmp.h  = SUBCcc (rem.h,  dvsr.h, cy, t0, t1, t2);
        if (!cy) { // remainder >= divisor
            rem.l = tmp.l;
            rem.h = tmp.h;
            quot.l = quot.l | 1;
        }
        bits_left--;
    } while (bits_left);
    return quot;
}

/*
  https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
  From: geo <[email protected]>
  Newsgroups: sci.math,comp.lang.c,comp.lang.fortran
  Subject: 64-bit KISS RNGs
  Date: Sat, 28 Feb 2009 04:30:48 -0800 (PST)

  This 64-bit KISS RNG has three components, each nearly
  good enough to serve alone.    The components are:
  Multiply-With-Carry (MWC), period (2^121+2^63-1)
  Xorshift (XSH), period 2^64-1
  Congruential (CNG), period 2^64
*/
static uint64_t kiss64_x = 1234567890987654321ULL;
static uint64_t kiss64_c = 123456123456123456ULL;
static uint64_t kiss64_y = 362436362436362436ULL;
static uint64_t kiss64_z = 1066149217761810ULL;
static uint64_t kiss64_t;
#define MWC64  (kiss64_t = (kiss64_x << 58) + kiss64_c, \
                kiss64_c = (kiss64_x >> 6), kiss64_x += kiss64_t, \
                kiss64_c += (kiss64_x < kiss64_t), kiss64_x)
#define XSH64  (kiss64_y ^= (kiss64_y << 13), kiss64_y ^= (kiss64_y >> 17), \
                kiss64_y ^= (kiss64_y << 43))
#define CNG64  (kiss64_z = 6906969069ULL * kiss64_z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)

int main (void)
{
    uint64_t a, b, res, ref;
    my_uint64 aa, bb, rr;
    do {
        a = KISS64;
        b = KISS64;
        ref = a / b;

        aa.l = (uint32_t)a;
        aa.h = (uint32_t)(a >> 32);
        bb.l = (uint32_t)b;
        bb.h = (uint32_t)(b >> 32);
        rr = bitwise_division_64 (aa, bb);
        res = (((uint64_t)rr.h) << 32) + rr.l;

        if (ref != res) {
            printf ("a=%016llx b=%016llx res=%016llx ref=%016llx\n", a, b, res, ref);
            return EXIT_FAILURE;
        }
    } while (a);
    return EXIT_SUCCESS;
}

A faster approach than bit-wise computation is to compute the reciprocal of the divisor, multiply by the dividend resulting in a preliminary quotient, then compute the remainder to precisely adjust the quotient. The entire computation can be accomplished in fixed-point arithmetic. However, on modern processors with fast floating-point units it is more convenient to generate a starting approximation for the reciprocal with a double-precision division. A single Halley iteration with cubic convergence then results in a full-precision reciprocal.

The Halley iteration for the reciprocal is very integer multiplication intensive, with a 64x64-bit multiply with 128-bit result (umul64wide() in the code below) being the building block crucial to performance. On modern 64-bit architectures this is typically a single machine instruction executing in a few cycles, however this is not accessible to portable code. Portable code emulating the instruction requires about 15 to 20 instructions depending on architecture and compiler.

The entire 128-bit division should take roughly 300 cycles, or ten times as fast as the simple bit-wise computation. Because the code is fairly complex, it requires a significant amount of testing to ensure correct operation. In the framework below I am using pattern-based and random tests for moderately intensive testing, using the straightforward bit-wise implementation as a reference.

The implementation of udiv128() below assumes that the programming enviornment uses IEEE-754 compliant floating-point arithmetic, that the double type is mapped to IEEE-754's binary64 type, and that division of double operands is correctly rounded.

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <limits.h>

typedef struct {
    uint64_t l;
    uint64_t h;
} my_uint128;

my_uint128 make_my_uint128 (uint64_t h, uint64_t l);
my_uint128 add128 (my_uint128 a, my_uint128 b);
my_uint128 sub128 (my_uint128 a, my_uint128 b);
my_uint128 lsl128 (my_uint128 a, int sh);
my_uint128 lsr128 (my_uint128 a, int sh);
my_uint128 not128 (my_uint128 a);
my_uint128 umul128lo (my_uint128 a, my_uint128 b);
my_uint128 umul128hi (my_uint128 a, my_uint128 b);
double my_uint128_to_double (my_uint128 a);
int lt128 (my_uint128 a, my_uint128 b);
int eq128 (my_uint128 a, my_uint128 b);
uint64_t double_as_uint64 (double a);
double uint64_as_double (uint64_t a);

#define FP64_EXPO_BIAS   (1023)
#define FP64_MANT_BITS   (53)
#define FP64_MANT_IBIT   (0x0010000000000000ULL)
#define FP64_MANT_MASK   (0x000fffffffffffffULL)
#define FP64_INC_EXP_128 (0x0800000000000000ULL)
#define FP64_MANT_ADJ    (2)  // adjustment to ensure underestimate

my_uint128 udiv128 (my_uint128 dividend, my_uint128 divisor)
{
    const my_uint128 zero = make_my_uint128 (0ULL, 0ULL);
    const my_uint128 one  = make_my_uint128 (0ULL, 1ULL);
    const my_uint128 two  = make_my_uint128 (0ULL, 2ULL);
    my_uint128 recip, temp, quo, rem;
    my_uint128 neg_divisor = sub128 (zero, divisor);
    double r;

    /* compute initial approximation for reciprocal; must be underestimate! */
    r = 1.0 / my_uint128_to_double (divisor);
    uint64_t i = double_as_uint64 (r) - FP64_MANT_ADJ + FP64_INC_EXP_128;
    temp = make_my_uint128 (0ULL, (i & FP64_MANT_MASK) | FP64_MANT_IBIT);
    int sh = (i >> (FP64_MANT_BITS-1)) - FP64_EXPO_BIAS - (FP64_MANT_BITS-1);
    recip = (sh < 0) ? lsr128 (temp, -sh) : lsl128 (temp, sh);

    /* perform Halley iteration with cubic convergence to refine reciprocal */
    temp = umul128lo (neg_divisor, recip);
    temp = add128 (umul128hi (temp, temp), temp);
    recip = add128 (umul128hi (recip, temp), recip);

    /* compute preliminary quotient and remainder */
    quo = umul128hi (dividend, recip); 
    rem = sub128 (dividend, umul128lo (divisor, quo));

    /* adjust quotient if too small; quotient off by 2 at most */
    if (! lt128 (rem, divisor)) {
        quo = add128 (quo, lt128 (sub128 (rem, divisor), divisor) ? one : two);
    }

    /* handle division by zero */
    if (eq128 (divisor, zero)) quo = not128 (zero);

    return quo;
}

#define SUBCcc(a,b,cy,t0,t1,t2) \
  (t0=(b)+cy, t1=(a), cy=t0<cy, t2=t1<t0, cy=cy+t2, t1-t0)

#define SUBcc(a,b,cy,t0,t1) \
  (t0=(b), t1=(a), cy=t1<t0, t1-t0)

#define SUBC(a,b,cy,t0,t1) \
  (t0=(b)+cy, t1=(a), t1-t0)

#define ADDCcc(a,b,cy,t0,t1) \
  (t0=(b)+cy, t1=(a), cy=t0<cy, t0=t0+t1, t1=t0<t1, cy=cy+t1, t0=t0)

#define ADDcc(a,b,cy,t0,t1) \
  (t0=(b), t1=(a), t0=t0+t1, cy=t0<t1, t0=t0)

#define ADDC(a,b,cy,t0,t1) \
  (t0=(b)+cy, t1=(a), t0+t1)

uint64_t double_as_uint64 (double a) 
{ 
    uint64_t r; 
    memcpy (&r, &a, sizeof r); 
    return r; 
}

double uint64_as_double (uint64_t a) 
{ 
    double r; 
    memcpy (&r, &a, sizeof r); 
    return r; 
}

my_uint128 add128 (my_uint128 a, my_uint128 b)
{
    uint64_t cy, t0, t1;
    a.l = ADDcc (a.l, b.l, cy, t0, t1);
    a.h = ADDC  (a.h, b.h, cy, t0, t1);
    return a;
}

my_uint128 sub128 (my_uint128 a, my_uint128 b)
{
    uint64_t cy, t0, t1;
    a.l = SUBcc (a.l, b.l, cy, t0, t1);
    a.h = SUBC  (a.h, b.h, cy, t0, t1);
    return a;
}

my_uint128 lsl128 (my_uint128 a, int sh)
{
    if (sh >= 64) {
        a.h = a.l << (sh - 64);
        a.l = 0ULL;
    } else if (sh) {
        a.h = (a.h << sh) + (a.l >> (64 - sh));
        a.l = a.l << sh;
    }
    return a;
}

my_uint128 lsr128 (my_uint128 a, int sh)
{
    if (sh >= 64) {
        a.l = a.h >> (sh - 64);
        a.h = 0ULL;
    } else if (sh) {
        a.l = (a.l >> sh) + (a.h << (64 - sh));
        a.h = a.h >> sh;
    } 
    return a;
}

my_uint128 not128 (my_uint128 a)
{
    a.l = ~a.l;
    a.h = ~a.h;
    return a;
}

int lt128 (my_uint128 a, my_uint128 b)
{
    uint64_t cy, t0, t1, t2;
    a.l = SUBcc  (a.l, b.l, cy, t0, t1);
    a.h = SUBCcc (a.h, b.h, cy, t0, t1, t2);
    return cy;
}

int eq128 (my_uint128 a, my_uint128 b)
{
    return (a.l == b.l) && (a.h == b.h);
}

// derived from Hacker's Delight 2nd ed. figure 8-2
my_uint128 umul64wide (uint64_t u, uint64_t v)
{
    my_uint128 r;
    uint64_t u0, v0, u1, v1, w0, w1, w2, t;
    u0 = (uint32_t)u;  u1 = u >> 32;
    v0 = (uint32_t)v;  v1 = v >> 32;
    w0 = u0 * v0;
    t  = u1 * v0 + (w0 >> 32);
    w1 = (uint32_t)t;
    w2 = t >> 32;
    w1 = u0 * v1 + w1;
    r.h = u1 * v1 + w2 + (w1 >> 32);
    r.l = (w1 << 32) + (uint32_t)w0;
    return r;
}

my_uint128 make_my_uint128 (uint64_t h, uint64_t l)
{
    my_uint128 r;
    r.h = h;
    r.l = l;
    return r;
}

my_uint128 umul128lo (my_uint128 a, my_uint128 b)
{
    my_uint128 r;
    r = umul64wide (a.l, b.l);
    r.h = r.h + a.l * b.h + a.h * b.l;
    return r;
}

my_uint128 umul128hi (my_uint128 a, my_uint128 b)
{
    my_uint128 t0, t1, t2, t3;
    t0 = umul64wide (a.l, b.l);
    t3 = add128 (umul64wide (a.h, b.l), make_my_uint128 (0ULL, t0.h));
    t1 = make_my_uint128 (0ULL, t3.l);
    t2 = make_my_uint128 (0ULL, t3.h);
    t1 = add128 (umul64wide (a.l, b.h), t1);
    return add128 (add128 (umul64wide (a.h, b.h), t2), make_my_uint128 (0ULL, t1.h));
}

double my_uint128_to_double (my_uint128 a)
{
    const int intbits = sizeof (a) * CHAR_BIT;
    const my_uint128 zero = make_my_uint128 (0ULL, 0ULL);
    my_uint128 rnd, i = a;
    uint64_t j;
    int sh = 0;
    double r;

    // normalize integer so MSB is set
    if (lt128 (i, make_my_uint128(0x0000000000000001ULL, 0))) {i = lsl128 (i,64); sh += 64; }
    if (lt128 (i, make_my_uint128(0x0000000100000000ULL, 0))) {i = lsl128 (i,32); sh += 32; }
    if (lt128 (i, make_my_uint128(0x0001000000000000ULL, 0))) {i = lsl128 (i,16); sh += 16; }
    if (lt128 (i, make_my_uint128(0x0100000000000000ULL, 0))) {i = lsl128 (i, 8); sh +=  8; } 
    if (lt128 (i, make_my_uint128(0x1000000000000000ULL, 0))) {i = lsl128 (i, 4); sh +=  4; }
    if (lt128 (i, make_my_uint128(0x4000000000000000ULL, 0))) {i = lsl128 (i, 2); sh +=  2; }
    if (lt128 (i, make_my_uint128(0x8000000000000000ULL, 0))) {i = lsl128 (i, 1); sh +=  1; }
    // form mantissa with explicit integer bit 
    rnd = lsl128 (i, FP64_MANT_BITS);
    i = lsr128 (i, intbits - FP64_MANT_BITS);
    j = i.l;
    // add in exponent, taking into account integer bit of mantissa
    if (! eq128 (a, zero)) {
        j += (uint64_t)(FP64_EXPO_BIAS + (intbits-1) - 1 - sh) << (FP64_MANT_BITS-1);
    }
    // round to nearest or even
    rnd.h = rnd.h | (rnd.l != 0);
    if ((rnd.h > 0x8000000000000000ULL) || 
        ((rnd.h == 0x8000000000000000ULL) && (j & 1))) j++;
    // reinterpret bit pattern as IEEE-754 'binary64'
    r = uint64_as_double (j);
    return r;
}

my_uint128 bitwise_division_128 (my_uint128 dvnd, my_uint128 dvsr)
{
    my_uint128 quot, rem, tmp;
    uint64_t cy, t0, t1, t2;
    int bits_left = CHAR_BIT * sizeof (dvsr);
    
    quot.h = dvnd.h;
    quot.l = dvnd.l;
    rem.h = 0;
    rem.l = 0;
    do {
        quot.l = ADDcc  (quot.l, quot.l, cy, t0, t1);
        quot.h = ADDCcc (quot.h, quot.h, cy, t0, t1);
        rem.l  = ADDCcc (rem.l,  rem.l,  cy, t0, t1);
        rem.h  = ADDC   (rem.h,  rem.h,  cy, t0, t1);
        tmp.l  = SUBcc  (rem.l,  dvsr.l, cy, t0, t1);
        tmp.h  = SUBCcc (rem.h,  dvsr.h, cy, t0, t1, t2);
        if (!cy) { // remainder >= divisor
            rem.l = tmp.l;
            rem.h = tmp.h;
            quot.l = quot.l | 1;
        }
        bits_left--;
    } while (bits_left);
    return quot;
}

/*
  https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
  From: geo <[email protected]>
  Newsgroups: sci.math,comp.lang.c,comp.lang.fortran
  Subject: 64-bit KISS RNGs
  Date: Sat, 28 Feb 2009 04:30:48 -0800 (PST)

  This 64-bit KISS RNG has three components, each nearly
  good enough to serve alone.    The components are:
  Multiply-With-Carry (MWC), period (2^121+2^63-1)
  Xorshift (XSH), period 2^64-1
  Congruential (CNG), period 2^64
*/
static uint64_t kiss64_x = 1234567890987654321ULL;
static uint64_t kiss64_c = 123456123456123456ULL;
static uint64_t kiss64_y = 362436362436362436ULL;
static uint64_t kiss64_z = 1066149217761810ULL;
static uint64_t kiss64_t;
#define MWC64  (kiss64_t = (kiss64_x << 58) + kiss64_c, \
                kiss64_c = (kiss64_x >> 6), kiss64_x += kiss64_t, \
                kiss64_c += (kiss64_x < kiss64_t), kiss64_x)
#define XSH64  (kiss64_y ^= (kiss64_y << 13), kiss64_y ^= (kiss64_y >> 17), \
                kiss64_y ^= (kiss64_y << 43))
#define CNG64  (kiss64_z = 6906969069ULL * kiss64_z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)

my_uint128 v[100000]; /* FIXME: size appropriately */

int main (void)
{
    const my_uint128 zero = make_my_uint128 (0ULL, 0ULL);
    const my_uint128 one = make_my_uint128 (0ULL, 1ULL);
    my_uint128 dividend, divisor, quot, ref;
    int i, j, patterns, idx = 0, nbrBits = sizeof (v[0]) * CHAR_BIT;
    int patterns_done = 0;

    /* pattern class 1: 2**i */
    for (i = 0; i < nbrBits; i++) {
        v [idx] = lsl128 (one, i);
        idx++;
    }
    /* pattern class 2: 2**i-1 */
    for (i = 0; i < nbrBits; i++) {
        v [idx] = sub128 (lsl128 (one, i), one);
        idx++;
    }
    /* pattern class 3: 2**i+1 */
    for (i = 0; i < nbrBits; i++) {
        v [idx] = add128 (lsl128 (one, i), one); 
        idx++;
    }
    /* pattern class 4: 2**i + 2**j */
    for (i = 0; i < nbrBits; i++) {
        for (j = 0; j < nbrBits; j++) {
            v [idx] = add128 (lsl128 (one, i), lsl128 (one, j));
            idx++;
        }
    }
    /* pattern class 5: 2**i - 2**j */
    for (i = 0; i < nbrBits; i++) {
        for (j = 0; j < nbrBits; j++) {
            v [idx] = sub128 (lsl128 (one, i), lsl128 (one, j));
            idx++;
        }
    }
    patterns = idx;
    /* pattern class 6: one's complement of pattern classes 1 through 5 */
    for (i = 0; i < patterns; i++) {
        v [idx] = not128 (v [i]);
        idx++;
    }
    /* pattern class 7: two's complement of pattern classes 1 through 5 */
    for (i = 0; i < patterns; i++) {
        v [idx] = sub128 (zero, v[i]);
        idx++;
    }
    patterns = idx;
    printf ("Starting pattern-based tests. Number of patterns: %d\n", patterns);

    for (long long int k = 0; k < 100000000000LL; k++) {
        if (k < patterns * patterns) {
            dividend = v [k / patterns];
            divisor  = v [k % patterns];
        } else {
            if (!patterns_done) {
                printf ("Starting random tests\n");
                patterns_done = 1;
            }
            dividend.l = KISS64;
            dividend.h = KISS64;
            divisor.h  = KISS64;
            divisor.l  = KISS64;
        }
        /* exclude cases with undefined results: division by zero */
        if (! eq128 (divisor, zero)) {
            quot = udiv128 (dividend, divisor);
            ref = bitwise_division_128 (dividend, divisor);
            if (! eq128 (quot, ref)) {
                printf ("@ (%016llx_%016llx, %016llx_%016llx): quot = %016llx_%016llx  ref=%016llx_%016llx\n", 
                        dividend.h, dividend.l, divisor.h, divisor.l, 
                        quot.h, quot.l, ref.h, ref.l);
                return EXIT_FAILURE;
            }
        }
    }
    printf ("unsigned 128-bit division: tests passed\n");
    return EXIT_SUCCESS;
}

Dessau

Awesome answers, thanks! I'd give you more up-votes if I could :) – Temperate

@Temperate To give you an idea: The answer as it is now represents about 10 hours of work. The design space for integer division is quite large; one could easily dedicate an entire master thesis to it. The two main approaches are digit-by-digit computations (from binary digits to high-radix as suggested in a comment, e.g. Knuth's Algorithm D in TAOCP Vol. 2, or Pope-Stein) and methods based on functional iteration, which additionally involve the problem of generating a starting approximation. Here I mostly side-stepped the latter can of worms by using double division as the starting point. – Dessau