单词分割算法

Question

单词分割算法

4

我开发嵌入式平台软件，需要一个单词分割算法。问题如下：给定由一系列32位字表示的大整数（可能很多），我们需要将其除以另一个32位字，即计算商（也是大整数）和余数（32位）。当然，如果我在x86上开发此算法，我可以简单地使用GNU MP但这个库对于嵌入式平台来说太大了。此外，我们的处理器没有硬件整数除法器（整数除法在软件中执行）。然而，处理器具有相当快的FPU，因此诀窍是尽可能使用浮点运算。有任何实现方法吗？

- user1545642

简单的恢复除法是否足够好？ - harold

是的，只需要在O(n)时间内工作的普通除法与余数即可，不需要渐进快速的方法。但我假设二进制（基数2）恢复除法是在单个位上操作的，而我希望避免这种情况，即最好使用整个字的除法。 - user1545642

但是整个单词分割从哪里来的呢 - FPU？浮点数的有效位应该有64位，或者使用半字。 - harold

2

我不知道，这对我来说听起来有些不安全。无论如何，这应该可以工作：http://www.hackersdelight.org/HDcode/divmnu.c.txt - harold

也许我会看一下GNU MP中是如何实现的，但我认为他们在那里不使用浮点数。 - user1545642

显示剩余3条评论

3个回答

1

听起来像是一个经典的优化。不要除以 D，而是乘以 0x100000000/D，然后再除以 0x100000000。后者只是一个字移位，也就是微不足道的。计算乘数有点难，但并不是很难。

此外，还可以参考这篇详细的文章，了解更详细的背景信息。

- MSalters

谢谢这篇文章，我学到了一些东西。虽然如果我理解正确的话，要执行64/32->32基本情况下的除法，我们需要将一个64位整数乘以一个固定点数2^32/D。看起来这需要一对32x32->64乘法，而我只能通过内联汇编获得？ - user1545642

1

我相信查找表和牛顿-拉夫逊连续逼近法是硬件设计师的经典选择（通常无法承担完整硬件除法所需的门）。您可以选择在精度和执行时间之间进行权衡。

- marko

@Marco：就像我说的，我们没有硬件整数除法器，但我们确实拥有相当快的浮点运算，因此不利用它是愚蠢的。使用迭代算法可能并不比使用基数2恢复除法算法更好。 - user1545642

牛顿-拉弗森方法避免了硬件除法 - 只需要乘法和减法（以及查找表）即可。 - marko

是的，但是对于64/32->32除法，您需要执行多少次迭代？这是单词除法算法的基本情况。此外，对于64位被除数，使用双精度的牛顿-拉弗森方法，我们再次遇到四舍五入的问题，因为dp尾数为53位。 - user1545642

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user1545642 · Accepted Answer

看看这个：该算法使用浮点数进行64x32->32除法，将整数a [0..n-1]除以单个字'c'。商'q'的limbs只是在循环中打印，如果您喜欢，可以将它们保存在数组中。请注意，您不需要GMP来运行算法-我只是用它来比较结果。

#include <gmp.h>

// divides a multi-precision integer a[0..n-1] by a single word c
void div_by_limb(const unsigned *a, unsigned n, unsigned c) {

  typedef unsigned long long uint64;
  unsigned c_norm = c, sh = 0;
  while((c_norm & 0xC0000000) == 0) { // make sure the 2 MSB are set
     c_norm <<= 1; sh++;
  }
  // precompute the inverse of 'c'
   double inv1 = 1.0 / (double)c_norm, inv2 = 1.0 / (double)c;
   unsigned i, r = 0;

   printf("\nquotient: "); // quotient is printed in a loop
   for(i = n - 1; (int)i >= 0; i--) { // start from the most significant digit
      unsigned u1 = r, u0 = a[i];
      union {
        struct { unsigned u0, u1; };
        uint64 x;
      } s = {u0, u1}; // treat [u1, u0] as 64-bit int
      // divide a 2-word number [u1, u0] by 'c_norm' using floating-point
      unsigned q = floor((double)s.x * inv1), q2;
      r = u0 - q * c_norm;
      // divide again: this time by 'c'
      q2 = floor((double)r * inv2);

      q = (q << sh) + q2; // reconstruct the quotient
      printf("%x", q);
  }

  r %= c; // adjust the residue after normalization
  printf("; residue: %x\n", r);
}

int main() {
  mpz_t z, quo, rem;
  mpz_init(z); // this is a dividend
  mpz_set_str(z, "9999999999999999999999999999999", 10);
  unsigned div = 9; // this is a divisor
  div_by_limb((unsigned *)z->_mp_d, mpz_size(z), div);
  mpz_init(quo); mpz_init(rem);
  mpz_tdiv_qr_ui(quo, rem, z, div); // divide 'z' by 'div'
  gmp_printf("compare: Quo: %Zx; Rem %Zx\n", quo, rem);
  mpz_clear(quo);
  mpz_clear(rem);
  mpz_clear(z);
  return 1;
}