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.
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))