[Ada] Correct multi-precision division algorithm used for universal integers

[Ada] Correct multi-precision division algorithm used for universal integers

[Arnaud Charlet]

[-- Attachment #1: Type: text/plain, Size: 625 bytes --]

The version of multi-precision division algorithm used was based on Algorithm D
published in 2nd edition of Donald Knuth's "The Art of Computer
Programming". This algorithm was corrected twice, in 1995 and 2005, to protect
against a possible overflow. Although this problem may not be present in this
code, due to the value of Base used with respect to the underlying machine
integers, it is preferable to use the corrected version of the algorithm.

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

2012-12-05  Yannick Moy  <moy@adacore.com>

	* uintp.adb (UI_Div_Rem): Correct algorithm D to remove potential
	overflow.


[-- Attachment #2: difs --]
[-- Type: text/plain, Size: 2207 bytes --]

Index: uintp.adb
===================================================================
--- uintp.adb	(revision 194188)
+++ uintp.adb	(working copy)
@@ -1165,6 +1165,7 @@
             Divisor_Dig1 : Int;
             Divisor_Dig2 : Int;
             Q_Guess      : Int;
+            R_Guess      : Int;
 
          begin
             --  [ NORMALIZE ] (step D1 in the algorithm). First calculate the
@@ -1218,30 +1219,26 @@
 
                --  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.
+               --  See: http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz
+               --  and http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz.
+               --  The code below is the fixed version of this step.
 
-               --  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
 
-               if Dividend (J) = Divisor_Dig1 then
-                  Q_Guess := Base - 1;
-               else
-                  Q_Guess := Tmp_Int / Divisor_Dig1;
-               end if;
+               Q_Guess := Tmp_Int / Divisor_Dig1;
+               R_Guess := Tmp_Int rem Divisor_Dig1;
 
                --  Refine the guess
 
-               while Divisor_Dig2 * Q_Guess >
-                     (Tmp_Int - Q_Guess * Divisor_Dig1) * Base +
-                        Dividend (J + 2)
+               while Q_Guess >= Base
+                 or else Divisor_Dig2 * Q_Guess >
+                           R_Guess * Base + Dividend (J + 2)
                loop
                   Q_Guess := Q_Guess - 1;
+                  R_Guess := R_Guess + Divisor_Dig1;
+                  exit when R_Guess >= Base;
                end loop;
 
                --  [ MULTIPLY & SUBTRACT ] (step D4). Q_Guess * Divisor is