I'm trying to implement a function that computes the division between two numbers approximately of the same size. I use an array of unsigned int's to store the values of each number (see below the class structure). I want to implement an efficient division function for the case when the difference between the number of unsigned int's of both is not bigger than 1. I've looked into the long division algorithm (as described in wikipedia), but it's too slow. Is there any efficient algorithm that deals with this special case? At the moment, I'm dealing with numbers with up to 100 digits.
class BigInt {
public:
short sign; //sign of the number
unsigned int size; //number of unsigned int's needed to store the number
unsigned int* values; //values of the number
...
}
Thanks
I believe a modified version of Knuth's algorithm D (The Art of Computer Programming Vol.2, Section 4.3.1) ought to do the trick in constant average time for random input, with a linear time worst-case. Significant simplifications stem from only requiring the first quotient word produced by the initial iteration.
Below is my attempt at unashamedly adapting the generic implementation from Hacker's Delight. The implementation assumes 32-bit words and 64-bit intermediates/division but may be naturally extended where wider arithmetic is available.
I do fear that the function is virtually untested however so please treat this more as the germ of an idea as opposed to a fully-fledged implementation. To be honest I don't even remember the subtleties of why the algorithm works.
// A helper function for finding the position of the most significant set bit.
// Please plug in your intrinsic of choice here
static unsigned int find_first_bit(uint32_t value) {
# ifdef _MSC_VER
unsigned long index;
(void) _BitScanReverse(&index, value);
return index + 1;
# else
unsigned int count = 0;
for(count = 0; value; ++count)
value >>= 1;
return count;
# endif
}
// Multi-word division in 32-bit big-endian segments, returning the most significant
// word of the quotient.
uint32_t divide(const uint32_t *num, const uint32_t *den, size_t len) {
// Skip past leading zero denominator digits. The quotient must fit in a single
// 32-bit word and so only the preceeding numerator digit needs be examined
uint32_t norm_den;
uint64_t norm_num = 0;
size_t top = 0;
while(norm_den = den[top], !norm_den)
norm_num = num[top++];
// Please pad the input to insure at least three denominator words counting from
// the first non-zero digit
assert(len >= top + 3);
// Divide the first two digits of the numerator by the leading digit of the
// denominator as an initial quotient digit guess, yielding an upper bound
// for the quotient at most two steps above the true value.
// Simultaneously normalize the denominator with the MSB in the 31st bit.
unsigned int norm = find_first_bit(norm_den);
norm_num = norm_num << (64 - norm);
norm_num |= ((uint64_t) num[top + 0] << 32 | num[top + 1]) >> norm;
norm_den = ((uint64_t) norm_den << 32 | den[top + 1]) >> norm;
// We are using a straight 64/64 division where 64/32=32 would suffice. The latter
// is available on e.g. x86-32 but difficult to generate short of assembly code.
uint32_t quot = (uint32_t) (norm_num / norm_den);
// Substitute norm_num - quot * norm_den if your optimizer is too thick-headed to
// efficiently extract the remainder
uint32_t rem = norm_num % norm_den;
// Test the next word of the input, reducing the upper bound to within one step
// of the true quotient. See Knuth for proofs of this reduction and the bounds
// of the first guess
norm_num = ((uint64_t) num[top + 1] << 32 | num[top + 2]) >> norm;
norm_num = (uint64_t) rem << 32 | (uint32_t) norm_num;
norm_den = ((uint64_t) den[top + 1] << 32 | den[top + 2]) >> norm;
if((uint64_t) quot * norm_den > norm_num) {
--quot;
// There is no "add-back" step try to avoid and so there is little point
// in looping to refine the guess further since the bound is sufficiently
// tight already
}
// Compare quotient guess multiplied by the denominator to the numerator
// across whole numbers to account for the final quotient step.
// There is no need to bother with normalization here. Furthermore we can
// compare from the most to the least significant and cut off early when the
// intermediate result becomes large, yielding a constant average running
// time for random input
uint64_t accum = 0;
do {
uint64_t prod = (uint64_t) quot * *den++;
accum = accum << 32 | *num++;
// A negative partial result can never recover, so pick the lower
// quotient. A separate test is required to avoid 65-bit arithmetic
if(accum < prod)
return --quot;
accum -= prod;
// Similarly a partial result which spills into the upper 32-bits can't
// recover either, so go with the upper quotient
if((uint64_t) accum >= 0x100000000)
return quot;
} while(--len);
return quot;
}