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.