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

	* sem_prag.adb: Minor reformatting.

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.



git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@193217 138bc75d-0d04-0410-961f-82ee72b054a4

4 files changed, 58 insertions, 30 deletions

diff --git a/gcc/ada/ChangeLog b/gcc/ada/ChangeLog
index 632b6030c20..dd82c971064 100644
--- a/gcc/ada/ChangeLog
+++ b/gcc/ada/ChangeLog
@@ -1,3 +1,12 @@
+2012-11-06  Robert Dewar  <dewar@adacore.com>
+
+	* sem_prag.adb: Minor reformatting.
+
+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.
+
 2012-11-06  Hristian Kirtchev  <kirtchev@adacore.com>
 
 	* exp_prag.adb (Expand_Pragma_Loop_Assertion): Update the comment
diff --git a/gcc/ada/s-bignum.adb b/gcc/ada/s-bignum.adb
index cb7a4f6c1be..84e9c8795d7 100644
--- a/gcc/ada/s-bignum.adb
+++ b/gcc/ada/s-bignum.adb
@@ -776,7 +776,9 @@ package body System.Bignums is
 
          d    : DD;
          j    : Length;
-         qhat : SD;
+         qhat : DD;
+         rhat : DD;
+         temp : DD;
 
       begin
          --  Initialize data of left and right operands
@@ -847,26 +849,37 @@ package body System.Bignums is
          --  Loop through digits
 
          loop
-            --  D3. [Calculate qhat.] If uj = v1, set qhat to b-l; otherwise
-            --  set qhat to (uj,uj+1)/v1.
-
-            if u (j) = v1 then
-               qhat := -1;
-            else
-               qhat := SD ((u (j) & u (j + 1)) / DD (v1));
-            end if;
-
-            --  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.
+            --  Note: In the original printing, step D3 was as follows:
 
-            while DD (v2) * DD (qhat) >
-                   ((u (j) & u (j + 1)) -
-                     DD (qhat) * DD (v1)) * b + DD (u (j + 2))
+            --  D3. [Calculate qhat.] If uj = v1, set qhat to b-l; otherwise
+            --  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
+
+            --  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. [Calculate qhat.] Set qhat to (uj,uj+1)/v1 and rhat to
+            --  to (uj,uj+1) mod v1.
+
+            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 @@ package body System.Bignums is
             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 @@ package body System.Bignums is
                --  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
 
diff --git a/gcc/ada/sem_prag.adb b/gcc/ada/sem_prag.adb
index 325ca0c3bc3..c3f27e14f89 100644
--- a/gcc/ada/sem_prag.adb
+++ b/gcc/ada/sem_prag.adb
@@ -11289,16 +11289,13 @@ package body Sem_Prag is
          --------------------
 
          --  pragma Loop_Assertion
-         --    (   [Invariant =>] boolean_Expression
-         --      | [Invariant =>] boolean_Expression ,
-         --         Variant => TERMINATION_VARIANTS
-         --      |  Variant => TERMINATION_VARIANTS );
-         --
-         --  TERMINATION_VARIANTS ::=
-         --    ( TERMINATION_VARIANT {, TERMINATION_VARIANT} )
-         --
+         --        (  [Invariant =>] boolean_Expression );
+         --     |  ( [[Invariant =>] boolean_Expression ,]
+         --            Variant   =>
+         --              ( TERMINATION_VARIANT {, TERMINATION_VARIANT ) );
+
          --  TERMINATION_VARIANT ::= CHANGE_MODIFIER => discrete_EXPRESSION
-         --
+
          --  CHANGE_MODIFIER ::= Increasing | Decreasing
 
          when Pragma_Loop_Assertion => Loop_Assertion : declare
diff --git a/gcc/ada/uintp.adb b/gcc/ada/uintp.adb
index a98bd9f376b..0761f2df70b 100644
--- a/gcc/ada/uintp.adb
+++ b/gcc/ada/uintp.adb
@@ -1216,6 +1216,15 @@ package body Uintp is
 
                --  [ 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 @@ package body Uintp is
 
                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;