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

login
register
mail settings
Submitter Arnaud Charlet
Date Dec. 5, 2012, 10:51 a.m.
Message ID <20121205105111.GA12333@adacore.com>
Download mbox | patch
Permalink /patch/203841/
State New
Headers show

Comments

Arnaud Charlet - Dec. 5, 2012, 10:51 a.m.
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.

Patch

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