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 <de...@adacore.com> * s-bignum.adb (Div_Rem): Fix bug in step D3. * uintp.adb (UI_Div_Rem): Add comment on bug in step D3.
Index: uintp.adb =================================================================== --- uintp.adb (revision 193215) +++ uintp.adb (working copy) @@ -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; Index: s-bignum.adb =================================================================== --- s-bignum.adb (revision 193215) +++ s-bignum.adb (working copy) @@ -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