clarify estimate through explicit variable containing
authorLuke Kenneth Casson Leighton <lkcl@lkcl.net>
Fri, 22 Apr 2022 09:41:19 +0000 (10:41 +0100)
committerLuke Kenneth Casson Leighton <lkcl@lkcl.net>
Fri, 22 Apr 2022 09:41:26 +0000 (10:41 +0100)
top two digits

openpower/sv/biginteger/divmnu64.c

index 57d007798eac23968e88efccd332e3404be66ab7..06b0deb13a8009f55988838e145d70d8d4b86866 100644 (file)
@@ -82,11 +82,11 @@ int divmnu(unsigned q[], unsigned r[], const unsigned u[], const unsigned v[],
            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;
 
@@ -124,10 +124,12 @@ int divmnu(unsigned q[], unsigned r[], const unsigned u[], const unsigned v[],
 
     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;