muldiv 512x512

\gdef\a{\mathtt{a}} \gdef\b{\mathtt{b}}

Goal. Compute \mathtt{div512}(\a_0, \a_1, \b_0, \b_1) = \floor{\frac{\a_0 + \a_1 ⋅ 2^{256}}{\b_0 + \b_1 ⋅ 2^{256}}} where all arguments are integers in the range [0, 2^{256}).

contract Div512 {
    function div512(uint256 a0, uint256 a1, uint256 b0, uint256 b1)
    internal pure returns (uint256 result)
    {
        if(b1 == 0) {
            return div512x256(a0, a1, b0)
        }
        // ...
    }

Here we can use Knuth's algorithm D, see for example the implementation in OpenZKP.

If \mathtt{b}_1 = 0 we can use the previous method, so assume \mathtt{b}_1 ≥ 1. In this case the high limbs bound the answer:

\floor{\frac {\mathtt{a}_1} {\mathtt{b}_1 + 1} } ≤ \floor{\frac {\mathtt{a}_0 + \mathtt{a}_1 ⋅ 2^{256}} {\mathtt{b}_0 + \mathtt{b}_1 ⋅ 2^{256}} } ≤ \floor{\frac {\mathtt{a}_1 + 1} {\mathtt{b}_1} }

So we can start by subtracting

a' = a - \floor{\frac {\mathtt{a}_1} {\mathtt{b}_1 + 1} } ⋅ b

If a' < b we are done.