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 <m...@adacore.com>
* uintp.adb (UI_Div_Rem): Correct algorithm D to remove potential
overflow.
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
