Calculation of `absNc` involves adding the sign bit of the denominator. It incorrectly used bit 31 of the denominator instead of bit 63.
UT q2;
UT t;
T result_magic;
- int iters = 0;
absDenom = abs(denom);
- t = two_nminus1 + ((unsigned int)denom >> 31);
+ t = two_nminus1 + (UT(denom) >> bits_minus_1);
absNc = t - 1 - (t % absDenom); // absolute value of nc
p = bits_minus_1; // initialize p
q1 = two_nminus1 / absNc; // initialize q1 = 2^p / abs(nc)
do
{
- iters++;
p++;
q1 *= 2; // update q1 = 2^p / abs(nc)
r1 *= 2; // update r1 = rem(2^p / abs(nc))