// Compute estimate qhat of q[j] from top 2 digits
uint64_t dig2 = ((uint64_t)un[j + n] << 32) | un[j + n - 1];
if (un[j+n] >= vn[n-1]) {
- // rhat can be bigger than 32-bit when the division overflows,
+ // rhat can be bigger than 32-bit when the division overflows
+ qhat = UINT32_MAX;
rhat = dig2 - (uint64_t)UINT32_MAX * vn[n - 1];
} else {
qhat = dig2 / vn[n - 1]; // 64/32 divide
rhat = dig2 % vn[n - 1]; // 64/32 modulo
}
// use 3rd-from-top digit to obtain better accuracy
+ b = 1UL<<32;
while (rhat < b || qhat * vn[n - 2] > b * rhat + un[j + n - 2]) {
qhat = qhat - 1;
rhat = rhat + vn[n - 1];