具有近似相同大小的操作数的大整数除法
Big integer division with operands aproximately of the same size
我正在尝试实现一个函数来计算两个大小大致相同的数字之间的除法。我使用 unsigned int 的数组来存储每个数字的值（见下面的 class 结构）。我想为两者 unsigned int 的数量之差不大于 1 的情况实现一个有效的除法函数。我研究了长除法算法（如维基百科中所述），但是它太慢了。有没有处理这种特殊情况的有效算法？目前，我正在处理最多 100 位数字。
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
...
}
谢谢
假设您想要一个整数结果，您应该能够相当快地执行此操作。
用小数开头的小数可能最容易理解。
让我们考虑较小的数字，例如：123456789/54321876。其中一个比另一个多一位，所以最大可能的结果大约是 10。
如果我们将两边相除，那么我们（比如）在较大的数字中留下三位数字，在较小的数字中留下两位数字，我们得到 123/54，这会产生与我们相同的整数结果整个工作都非常精确。
在你的情况下，你可以做几乎相同的事情，除了不是十进制数字，你基本上会在较大的项目中保持三个精确的词，在较小的项目中保持两个词的精确度（实际上，即使这是保守的--2 和 1 个单词分别 可能 就足够了，但我现在懒得证明这一点，如果你这样做，你可能还需要更加小心四舍五入).
所以是的，您认为这是一种特殊情况，您可以根据几乎相同大小的数字来节省工作，这是正确的。
我相信 Knuth 算法 D 的修改版本（计算机编程艺术第 2 卷，第 4.3.1 节）应该可以在随机输入的恒定平均时间和线性时间最坏情况下解决问题。显着的简化源于只需要初始迭代产生的第一个商词。
以下是我对通用 implementation from Hacker's Delight 毫不掩饰地改编的尝试。该实现假定 32 位字和 64 位 intermediates/division，但可以在更广泛的算术可用的情况下自然扩展。
我确实担心该功能实际上未经测试，因此请将其更多地视为想法的萌芽，而不是成熟的实施。老实说，我什至不记得为什么算法有效的微妙之处。
// 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;
}
我正在尝试实现一个函数来计算两个大小大致相同的数字之间的除法。我使用 unsigned int 的数组来存储每个数字的值（见下面的 class 结构）。我想为两者 unsigned int 的数量之差不大于 1 的情况实现一个有效的除法函数。我研究了长除法算法（如维基百科中所述），但是它太慢了。有没有处理这种特殊情况的有效算法？目前，我正在处理最多 100 位数字。
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
...
}
谢谢
假设您想要一个整数结果，您应该能够相当快地执行此操作。
用小数开头的小数可能最容易理解。
让我们考虑较小的数字，例如：123456789/54321876。其中一个比另一个多一位，所以最大可能的结果大约是 10。
如果我们将两边相除，那么我们（比如）在较大的数字中留下三位数字，在较小的数字中留下两位数字，我们得到 123/54，这会产生与我们相同的整数结果整个工作都非常精确。
在你的情况下，你可以做几乎相同的事情，除了不是十进制数字，你基本上会在较大的项目中保持三个精确的词，在较小的项目中保持两个词的精确度（实际上，即使这是保守的--2 和 1 个单词分别 可能 就足够了，但我现在懒得证明这一点，如果你这样做，你可能还需要更加小心四舍五入).
所以是的，您认为这是一种特殊情况，您可以根据几乎相同大小的数字来节省工作，这是正确的。
我相信 Knuth 算法 D 的修改版本（计算机编程艺术第 2 卷，第 4.3.1 节）应该可以在随机输入的恒定平均时间和线性时间最坏情况下解决问题。显着的简化源于只需要初始迭代产生的第一个商词。
以下是我对通用 implementation from Hacker's Delight 毫不掩饰地改编的尝试。该实现假定 32 位字和 64 位 intermediates/division，但可以在更广泛的算术可用的情况下自然扩展。
我确实担心该功能实际上未经测试，因此请将其更多地视为想法的萌芽，而不是成熟的实施。老实说，我什至不记得为什么算法有效的微妙之处。
// 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;
}