Re: Dividing 512 bit number by 128 bit number in C program

vaidehikedlaya@xxxxxxxxx wrote:
Thanks a lot this helps. Could you send me your C++ code. I would like
to look at it.

I decided to post here for benefit of others who would make the same
mistake I made.

Description of the algorithm:

Division is naturally the most complicated of the basic big-integer
operations. The following is a C++ static member function of class
Integer that does big integer division. It takes dividend and divisor
by const reference, and quotient and remainder by reference,
calculating the quotient and the remainder before returning. There is
no return value.

Division by 0 results in an Integer::DIVISION_BY_ZERO exception.
Division by 1 results in remainder being set to 0 and quotient being
set to the dividend.

The function consists mostly of a large outer loop.

In one of the inner loops, the divisor is positioned so that its most
significant bit is just beneath the most significant bit of the
dividend. A bit-wise subtraction of divisor from dividend is performed
from right to left, starting with the low-order bit of the divisor,
borrowing bits if necessary. The positioning of the divisor is
important - it is not enough to simply align the high-order bits and
forget, as it might happen that the bit sequence of dividend involved
in the subtraction yields a number whose magnitude is insufficient to
accommodate a subtraction by the divisor. A small snippet of code
checks for this situation before allowing the subtraction to occur by
scanning the bits of both the dividend and the divisor from left to
right, and if it happens that a dividend bit is found to be 0 while the
corresponding divisor bit is 1, then naturally, subtraction is
impossible. In this case, a simple shift of divisor one bit to the
right will remedy, but if the divisor cannot be shifted because it is
already as far to the right as it can go (the divisor_offset is 0)
then naturally we are done with algorithm.

As we perform the subtraction, we pay careful attention to borrowing,
taking all tri-state cases into consideration. Those familiar with the
subtraction logic of an ALU will find these cases familiar. I decided
to use IF statements to handle the case instead a truth table and
direct computation for clarity, since the algorithm is already to slow
for practical use.

Each time a subtraction occurs, we set a corresponding bit in the
quotient. This bit just so happens to have the same offset as
divsor_offset.

After each subtraction, if we find that there are fewer bits in the
dividend than divisor, we are done.

The function ::adjust_extent of class Integer is simply present to trim
the memory use of the Integer object.

Bit testing for 1 or 0 is done using array index into the buffer of the
Integer object and C/C++ logical operators.

The sign of the quotient is based on the sign of the dividend and the
divisor. The sign of the remainder, for lack of a better scheme, uses
the exact same logic to compute the sign of the quotient.

Again, this function should not be used in commercial code because it
is too slow. The only use would be to instruct students how *not* to
design a division function in C/C++. An efficient division algorithm
can be found in "D. Knuth, The Art of Computer Programming Vo.2,
Section 4.3.1"

-Le Chaud Lapin-


void Integer::divide (const Integer &dividend, const Integer &divisor,
Integer &quotient, Integer &remainder)
{
// Division by zero is not permitted:
if (divisor.is_zero())
throw Integer::DIVISION_BY_ZERO;

if (divisor.is_one())
{
quotient = dividend;
remainder = 0;
return;
}

quotient.nullify();
quotient.adjust_extent(dividend.extent());

// We copy the dividend into the remainder since the remainder will be
computed by successively reducing it, beginning with remainder ==
dividend:
remainder = dividend;

// Determine the initial number of bits in the remainder:
unsigned int remainder_length_bits = remainder.length_bits();

// Determine the initial number of bits in the divisor:
const unsigned int divisor_length_bits = divisor.length_bits();

// Keep trying to reduce remainder by dividing by divisor until there
are no more spots to shift divisor into (divisor_offset == 0):
while (remainder_length_bits >= divisor_length_bits)
{
unsigned int divisor_offset = remainder_length_bits -
divisor_length_bits;

// Check if the substraction of divisor can be performed:
bool subtractable = true;
unsigned int divisor_bit_index = divisor_length_bits - 1;

do // DEFECT: The high-order bits of both remainder and divisor are
always set so we can skip them.
{
register unsigned int remainder_bit_index = divisor_bit_index +
divisor_offset;

if (remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] & (1 <<
(remainder_bit_index % WORD_LENGTH_BITS)))
{
// Remainder bit is 1.
if (!(divisor.buffer[divisor_bit_index / WORD_LENGTH_BITS] & (1 <<
(divisor_bit_index % WORD_LENGTH_BITS))))
// Divisor bit is 0.
break; // ...subtraction can be performed.
// Else, both remainder bit and divisor bit are equal to 1.
}
else
{
// Remainder bit is 0.
if (divisor.buffer[divisor_bit_index / WORD_LENGTH_BITS] & (1 <<
(divisor_bit_index % WORD_LENGTH_BITS)))
{
// Divisor bit is 1.
// If we cannot shift the divisor to the right one position...
if (divisor_offset == 0)
subtractable = false; // ...subtraction cannot be performed.
// Otherwise...
else
--divisor_offset; // ...subtraction can be performed after
reduction of divisor_offset by 1:
break; // In any case, we break because we have remainder_bit ==
0, divisor_bit == 1, which is always an interesting situation.
}
// Else, both remainder bit and divisor bit are equal to 0.
}
}
while (divisor_bit_index--);

// If subtraction cannot be performed, we are finished:
if (!subtractable)
break;

// Perform a subtraction:
bool borrow = false;

for (divisor_bit_index = 0; divisor_bit_index < divisor_length_bits;
++divisor_bit_index)
{
register unsigned int remainder_bit_index = divisor_bit_index +
divisor_offset;

if (remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] & (1 <<
(remainder_bit_index % WORD_LENGTH_BITS)))
{
// Remainder bit is 1.
if (divisor.buffer[divisor_bit_index / WORD_LENGTH_BITS] & (1 <<
(divisor_bit_index % WORD_LENGTH_BITS)))
{
// Divisor bit is 1.
if (!borrow)
// Set remainder bit to 0 subtract (1-1) and satisfy borrow of 0:
remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] &= ~(1
<< (remainder_bit_index % WORD_LENGTH_BITS));
// Else leave remainder bit == 1 because borrow made it 0, but to
subtract divisor it need its own borrow, so (2-1) == 1
}
else
{
// Divisor bit is 0.
// In this case, we only disturb remainder bit if borrow needs to
be satisfied.
if (borrow)
{
// Set remainder bit to 0 to satisfy borrow:
remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] &= ~(1
<< (remainder_bit_index % WORD_LENGTH_BITS));
borrow = false;
}
}
}
else
{
// Remainder bit is 0.
if (divisor.buffer[divisor_bit_index / WORD_LENGTH_BITS] & (1 <<
(divisor_bit_index % WORD_LENGTH_BITS)))
{
// Divisor bit is 1.
if (!borrow)
{
// Set remainder bit to 1 to subtract (2-1). Indicate necessity
to borrow:
remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] |= (1 <<
(remainder_bit_index % WORD_LENGTH_BITS));
borrow = true;
}
// Else remainder bit should stay 0 because continuing borrow gave
us 2, borrow took 1 and divisor took 1.
}
else
{
// Divisor bit is 0.
if (borrow)
// Set remainder bit to 1 because continuing borrow gave us 2,
borrow took 1:
remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] |= (1 <<
(remainder_bit_index % WORD_LENGTH_BITS));
// Else there is nothing to do because all 3 states are 0.
}
}
}

if (borrow)
{
unsigned int remainder_bit_index = divisor_offset +
divisor_length_bits;
remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] &= ~(1 <<
(remainder_bit_index % WORD_LENGTH_BITS));
}

// Set appropriate quotient bit since subtraction was performed:
quotient.buffer[divisor_offset / WORD_LENGTH_BITS] |= (1 <<
(divisor_offset % WORD_LENGTH_BITS));

while (remainder_length_bits >= divisor_length_bits)
{
if (remainder.buffer[(remainder_length_bits - 1)/ WORD_LENGTH_BITS]
& (Integer::Word(1) << ((remainder_length_bits - 1)%
WORD_LENGTH_BITS)))
break;
--remainder_length_bits;
}
}

quotient.adjust_extent();
quotient.sign_ = static_cast<Integer::Sign>(dividend.sign_ *
divisor.sign_);

remainder.adjust_extent();
remainder.sign_ = static_cast<Integer::Sign>(dividend.sign_ *
divisor.sign_);
}

.

Relevant Pages

  • Re: Sign conventions for remainder
    ... Addition, subtraction, and multiplication ... What about if both the divisor and the dividend are negative? ... The remainder is a piece of the dividend; ...
    (sci.math)
  • Re: Sign conventions for remainder
    ... Addition, subtraction, and multiplication ... > the conventions for remainders. ... What about if both the divisor and the dividend are negative? ... the remainder should still be in the range ...
    (sci.math)
  • Re: Sign conventions for remainder
    ... What about if both the divisor and the dividend are negative? ... > If the divisor is positive, the remainder should still be in the range ... options in your package for the several possible sensible behaviours. ...
    (sci.math)
  • Re: How to perform integer div with SIMD?
    ... First, use the approximate reciprocal lookup to generate a set of starting values, then one or two NR iterations (also in SIMD fp mode) to get a more precise value. ... If negative, subtract one from the result, add divisor to remainder. ...
    (comp.lang.asm.x86)
  • Re: fixed point division in 32 bits
    ... i got some clue from a code posted by you. ... uint64_t remainder; ... uint64_t divisor; ... is your input data limited in magnitude? ...
    (comp.dsp)