d : DD;
j : Length;
- qhat : SD;
+ qhat : DD;
+ rhat : DD;
+ temp : DD;
begin
-- Initialize data of left and right operands
-- Loop through digits
loop
- -- D3. [Calculate qhat.] If uj = v1, set qhat to b-l; otherwise
- -- set qhat to (uj,uj+1)/v1.
-
- if u (j) = v1 then
- qhat := -1;
- else
- qhat := SD ((u (j) & u (j + 1)) / DD (v1));
- end if;
-
- -- 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.
+ -- Note: In the original printing, step D3 was as follows:
- while DD (v2) * DD (qhat) >
- ((u (j) & u (j + 1)) -
- DD (qhat) * DD (v1)) * b + DD (u (j + 2))
+ -- D3. [Calculate qhat.] If uj = v1, set qhat to b-l; otherwise
+ -- 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
+
+ -- 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. [Calculate qhat.] Set qhat to (uj,uj+1)/v1 and rhat to
+ -- to (uj,uj+1) mod v1.
+
+ 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
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
-- 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
--------------------
-- pragma Loop_Assertion
- -- ( [Invariant =>] boolean_Expression
- -- | [Invariant =>] boolean_Expression ,
- -- Variant => TERMINATION_VARIANTS
- -- | Variant => TERMINATION_VARIANTS );
- --
- -- TERMINATION_VARIANTS ::=
- -- ( TERMINATION_VARIANT {, TERMINATION_VARIANT} )
- --
+ -- ( [Invariant =>] boolean_Expression );
+ -- | ( [[Invariant =>] boolean_Expression ,]
+ -- Variant =>
+ -- ( TERMINATION_VARIANT {, TERMINATION_VARIANT ) );
+
-- TERMINATION_VARIANT ::= CHANGE_MODIFIER => discrete_EXPRESSION
- --
+
-- CHANGE_MODIFIER ::= Increasing | Decreasing
when Pragma_Loop_Assertion => Loop_Assertion : declare
-- [ 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
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;