Patchwork [Ada] Fix error in multi-precision division used in ELIMINATED mode

login
register
mail settings
Submitter Arnaud Charlet
Date Nov. 6, 2012, 9:56 a.m.
Message ID <20121106095623.GA6514@adacore.com>
Download mbox | patch
Permalink /patch/197433/
State New
Headers show

Comments

Arnaud Charlet - Nov. 6, 2012, 9:56 a.m.
The ELIMINATED overlow mode uses a multi-precision integer arithmetic
package that depends on algorithm D from Knuth for multi-precision
division. In very rare cases, this algorithm has a bug causing an
internal overflow resulting in an incorrect result. This patch uses
the fixed version of this algorithm (see code of Div_Rem in s-bignum
for details).

The following program:

     1. pragma Overflow_Checks (Eliminated);
     2. with Ada.Text_IO;
     3. procedure BadBigNum is
     4.     type U32 is mod 2**32;
     5.     type U64 is mod 2**64;
     6.     subtype LLI is Long_Long_Integer;
     7.
     8.     procedure Q (Y : LLI) is  --  Y is 2**32 - 1
     9.        Ym1 : LLI := Y - 1;  --  Ym1 is 2**32 - 2
    10.        Z : U64 := U64 ((((Ym1 * 2**32) + Ym1) * 2**32)
    11.                        / (Ym1 * 2**32 + Y));
    12.        --  Z is 2**32 - 1 = 4_294_967_295
    13.     begin
    14.        Ada.Text_IO.Put_Line ("z =" & U64'Image (Z));
    15.     end Q;
    16.
    17.     X : U32 := -1;  --  2**32 - 1
    18.     Y : LLI := LLI (X); --  2**32 - 1
    19. begin
    20.     Q (Y);
    21. end BadBigNum;

should generate output of

 4294967295

Tested on x86_64-pc-linux-gnu, committed on trunk

2012-11-06  Robert Dewar  <dewar@adacore.com>

	* s-bignum.adb (Div_Rem): Fix bug in step D3.
	* uintp.adb (UI_Div_Rem): Add comment on bug in step D3.

Patch

Index: uintp.adb
===================================================================
--- uintp.adb	(revision 193215)
+++ uintp.adb	(working copy)
@@ -1216,6 +1216,15 @@ 
 
                --  [ CALCULATE Q (hat) ] (step D3 in the algorithm)
 
+               --  Note: this version of step D3 is from the original published
+               --  algorithm, which is known to have a bug causing overflows.
+               --  See: http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz.
+               --  In this code we are safe since our representation of double
+               --  length numbers allows an expanded range.
+
+               --  We don't have a proof of this claim, but the only cases we
+               --  have found that show the bug in step D3 work fine here.
+
                Tmp_Int := Dividend (J) * Base + Dividend (J + 1);
 
                --  Initial guess
@@ -1230,7 +1239,7 @@ 
 
                while Divisor_Dig2 * Q_Guess >
                      (Tmp_Int - Q_Guess * Divisor_Dig1) * Base +
-                                                          Dividend (J + 2)
+                        Dividend (J + 2)
                loop
                   Q_Guess := Q_Guess - 1;
                end loop;
Index: s-bignum.adb
===================================================================
--- s-bignum.adb	(revision 193215)
+++ s-bignum.adb	(working copy)
@@ -776,7 +776,9 @@ 
 
          d    : DD;
          j    : Length;
-         qhat : SD;
+         qhat : DD;
+         rhat : DD;
+         temp : DD;
 
       begin
          --  Initialize data of left and right operands
@@ -847,26 +849,37 @@ 
          --  Loop through digits
 
          loop
+            --  Note: In the original printing, step D3 was as follows:
+
             --  D3. [Calculate qhat.] If uj = v1, set qhat to b-l; otherwise
-            --  set qhat to (uj,uj+1)/v1.
+            --  set qhat to (uj,uj+1)/v1. Now test if v2 * qhat is greater than
+            --  (uj*b + uj+1 - qhat*v1)*b + uj+2. If so, decrease qhat by 1 and
+            --  repeat this test
 
-            if u (j) = v1 then
-               qhat := -1;
-            else
-               qhat := SD ((u (j) & u (j + 1)) / DD (v1));
-            end if;
+            --  This had a bug not discovered till 1995, see Vol 2 errata:
+            --  http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz. Under
+            --  rare circumstances the expression in the test could overflow.
+            --  The code below is the fixed version of this step.
 
-            --  D3 (continued). Now test if v2 * qhat is greater than (uj*b +
-            --  uj+1 - qhat*v1)*b + uj+2. If so, decrease qhat by 1 and repeat
-            --  this test, which determines at high speed most of the cases in
-            --  which the trial value qhat is one too large, and it eliminates
-            --  all cases where qhat is two too large.
+            --  D3. [Calculate qhat.] Set qhat to (uj,uj+1)/v1 and rhat to
+            --  to (uj,uj+1) mod v1.
 
-            while DD (v2) * DD (qhat) >
-                   ((u (j) & u (j + 1)) -
-                     DD (qhat) * DD (v1)) * b + DD (u (j + 2))
+            temp := u (j) & u (j + 1);
+            qhat := temp / DD (v1);
+            rhat := temp mod DD (v1);
+
+            --  D3 (continued). Now test if qhat = b or v2*qhat > (rhat,uj+2):
+            --  if so, decrease qhat by 1, increase rhat by v1, and repeat this
+            --  test if rhat < b. [The test on v2 determines at at high speed
+            --  most of the cases in which the trial value qhat is one too
+            --  large, and eliminates all cases where qhat is two too large.]
+
+            while qhat = b
+              or else DD (v2) * qhat > LSD (rhat) & u (j + 2)
             loop
                qhat := qhat - 1;
+               rhat := rhat + DD (v1);
+               exit when rhat >= b;
             end loop;
 
             --  D4. [Multiply and subtract.] Replace (uj,uj+1..uj+n) by
@@ -892,7 +905,7 @@ 
             begin
                Borrow := 0;
                for K in reverse 1 .. n loop
-                  Temp := DD (qhat) * DD (v (K)) + DD (Borrow);
+                  Temp := qhat * DD (v (K)) + DD (Borrow);
                   Borrow := MSD (Temp);
 
                   if LSD (Temp) > u (j + K) then
@@ -908,7 +921,7 @@ 
                --  D5. [Test remainder.] Set qj = qhat. If the result of step
                --  D4 was negative, we will do the add back step (step D6).
 
-               q (j) := qhat;
+               q (j) := LSD (qhat);
 
                if Negative then