Solidity port of Hacker's Delight, Second Edition. FIGURE 9-3. Divide long unsigned, using fullword division instruction. Original version in C at https://gist.github.com/frangio/6d0f8e1297d0077449fac86033a43b12

divlu.sol

// SPDX-License-Identifier: MIT
pragma solidity ^0.8.0;

error Overflow();

uint constant b = 2 ** 128;

/// Divides a two word number by one word number.
/// Returns q and r such that u1 << 256 + u0 = q * v + r.
function divlu(uint u1, uint u0, uint v) pure returns (uint q, uint r) { unchecked {
    if (u1 >= v) {
        revert Overflow();
    }

    uint s = nlz(v);
    v = v << s;
    uint vn1 = v >> 128;
    uint vn0 = v & type(uint128).max;

    uint un32 = (u1 << s) | (u0 >> 256 - s) & uint(-int(s) >> 255);
    uint un10 = u0 << s;

    uint un1 = un10 >> 128;
    uint un0 = un10 & type(uint128).max;

    uint q1 = un32/vn1;
    uint rhat = un32 - q1*vn1;

    while (q1 >= b || q1*vn0 > b*rhat + un1) {
        q1 = q1 - 1;
        rhat = rhat + vn1;
        if (rhat >= b) break;
    }

    uint un21 = un32*b + un1 - q1*v;

    uint q0 = un21/vn1;
    rhat = un21 - q0*vn1;

    while (q0 >= b || q0*vn0 > b*rhat + un0) {
        q0 = q0 - 1;
        rhat = rhat + vn1;
        if (rhat >= b) break;
    }

    q = q1*b + q0;
    r = (un21*b + un0 - q0*v) >> s;
}}

function nlz(uint x) pure returns (uint n) {
    uint y;
    n = 256;
    y = x >> 128; if (y != 0) { n = n - 128; x = y; }
    y = x >>  64; if (y != 0) { n = n -  64; x = y; }
    y = x >>  32; if (y != 0) { n = n -  32; x = y; }
    y = x >>  16; if (y != 0) { n = n -  16; x = y; }
    y = x >>   8; if (y != 0) { n = n -   8; x = y; }
    y = x >>   4; if (y != 0) { n = n -   4; x = y; }
    y = x >>   2; if (y != 0) { n = n -   2; x = y; }
    y = x >>   1; if (y != 0) return n - 2;
    return n - x;
}