long-integerbigintegerinteger-divisiontaocp

Bug in the normalization step of Knuth Algorithm D (TAOCP 4.3.1)?


Knuth Algorithm D, during the normalization step (D1), states to set d to (b-1)//v_hi where b is the basis (word size or half-word size) and v_hi is the upper limb of the denominator. Then multiply the denominator by d, giving a multi-limb integer of the same size.

Knuth Algorithm D step D1 (TAOCP 4.3.1)

The problem is that the multiplication by d does NOT give a multi-limb integer of the same size. Sometimes it overflows. For instance, with b as 10 and v as 19, we obtain d as 9 and the new value for v as 171, which clearly overflows 2 limbs.

Am I misunderstanding the text of step D1 or is the prescription to set d to (b-1)//v_hi subtly incorrect? What would be the correct formulation for d?

The second half of step D1 states "On a binary computer it may be preferable to choose d to be a power of 2 instead of using the value suggested here". Unfortunately, I'm working on a platform that doesn't have a clz instruction and emulating it is extremely expensive. It would be faster to use simple truncating division to compute d.

Ultimately, I just need to compute a multiplier for v that gives the condition v_hi >= b//2. Also unfortunately, a lookup table would be too large on this platform.

Perhaps a more recent version of TAOCP fixes this issue, but I'm away from my books right now.


Solution

  • yes I agree with you, there are tricky corner cases where this overflows. I ran into the same bug, working on a 32-bit architecture one counter-example is v = (1, 2) (that is v = 2^32 + 2, v_1 = 1, v_0 = 2), the normalization step D1 computes d = floor((b-1)/v_1) that is 0xffffffff with b = 2^32. But then d * v = (2^32-1) * (2^32 + 2) = 2^64 + 2^32 - 2 = (1, 0, 0xfffffffe) with an overflow.

    edited answer: in the Errata of Volume 2, at Earliest errata for Volume 2 (3rd ed.), 08 January 2011, the following answer is given:

    floor((b-1)/v_{n-1}) -> floor(b/(v_{n-1}+1))