int m, int n)
{
- const unsigned long long b = 4294967296LL; // Number base (2**32).
- unsigned *un, *vn; // Normalized form of u, v.
- unsigned long long qhat; // Estimated quotient digit.
- unsigned long long rhat; // A remainder.
- unsigned long long p; // Product of two digits.
+ const unsigned long long b = 1LL<<32; // Number base (2**32).
+ unsigned *un, *vn; // Normalized form of u, v.
+ unsigned long long qhat; // Estimated quotient digit.
+ unsigned long long rhat; // A remainder.
+ unsigned long long p; // Product of two digits.
long long t, k;
int s, i, j;
for (j = m - n; j >= 0; j--)
{ // Main loop.
- // Compute estimate qhat of q[j].
- qhat = (un[j + n] * b + un[j + n - 1]) / vn[n - 1];
- rhat = (un[j + n] * b + un[j + n - 1]) - qhat * vn[n - 1];
+ // Compute estimate qhat of q[j] from top 2 digits
+ uint64_t dig2 = ((uint64_t)un[j + n] << 32) | un[j + n - 1];
+ qhat = dig2 / vn[n - 1];
+ rhat = dig2 % vn[n - 1];
again:
+ // use 3rd-from-top digit to obtain better accuracy
if (qhat >= b || qhat * vn[n - 2] > b * rhat + un[j + n - 2])
{
qhat = qhat - 1;