===================================================================
@@ -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;
===================================================================
@@ -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
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.