【问题标题】:single-word division algorithm单词分割算法
【发布时间】:2012-08-10 12:03:08
【问题描述】:

我是开发嵌入式平台的软件,需要单字除法算法。

问题如下: 给定一个由 32 位字(可以是很多)序列表示的大整数, 我们需要将它除以另一个 32 位字,即计算商(也是大整数) 以及余数(32 位)。

当然,如果我在 x86 上开发这个算法,我可以简单地使用 GNU MP 但是这个库对于嵌入式平台来说太大了。此外,我们的处理器 没有硬件整数除法器(整数除法在软件中执行)。 然而,处理器的 FPU 相当快,所以诀窍是尽可能使用浮点运算。

任何想法如何实现这一点?

【问题讨论】:

  • 简单的恢复除法就够了吗?
  • 是的,在 O(n) 时间内工作的余数正常除法就足够了,不需要渐近快速的方法.. 但我假设二进制(基数 2)恢复除法在单个位上运行,我希望避免,即更好地使用全词除法
  • 但是整个单词除法是从哪里来的——FPU?那么浮点数的有效位应该有 64 位,或者可能使用半字..
  • 我们有 53 位尾数的双精度算术。我假设我们不需要尾数中的整个 64 位来实现除法,因为 64 位整数乘以 1.0/x 浮点数将产生具有适当移位的 32 位结果。当然可能会出现一些舍入错误,但这就是如何做到这一点的诀窍..
  • 我不知道,这对我来说听起来有些不安全。这应该在任何情况下都有效:hackersdelight.org/HDcode/divmnu.c.txt

标签: c++ algorithm embedded


【解决方案1】:

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

有关更详细的背景信息,另请参阅 this detailed article

【讨论】:

  • 感谢这篇文章,我学到了一些东西。虽然如果我理解正确,要执行 64/32->32 基本情况除法,我们必须将 64 位整数乘以定点数 2^32/D。似乎这需要一对 32x32->64 乘法,我只能通过内联汇编器获得?
【解决方案2】:

看看这个:算法将整数 a[0..n-1] 除以单个单词 'c' 使用浮点进行 64x32->32 除法。商 'q' 的四肢只是打印在一个循环中,如果你愿意,你可以将其保存在一个数组中。请注意,您不需要 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;
}

【讨论】:

  • 我在支持 gmp 的 x86 Linux 上运行了您的程序,看起来运行良好。尽管我花了一些时间才意识到其中的逻辑。据我了解,您执行 2 个双精度乘法来实现 64/32->32 位除法,因为尾数舍入会搞砸事情......很酷的把戏,谢谢......而且它在主循环中不使用除法。顺便说一句,你为什么在开始时将 c_norm 向左移动?是为了避免分区溢出吗?
  • 你是对的,我使用 c_norm 来防止溢出:例如将 2^63 除以 3 时,结果会超过 32 位。相反,我将 2^63 除以 3*2^30,然后修复商和余数
【解决方案3】:

我相信查找表和Newton Raphson 逐次逼近是硬件设计人员使用的规范选择(他们通常买不起完整的硬件划分的大门)。您可以在准确性和执行时间之间进行权衡。

【讨论】:

  • @Marco:正如我所说,我们没有硬件整数除法器,但我们确实有相当快的浮点运算,因此不使用它会很愚蠢。使用迭代算法可能并不比仅使用基数 2 恢复除法算法更好
  • Newton Raphson 避免了硬件除法 - 只需要乘法和减法(以及 LUT)。
  • 是的,但是对于 64/32->32 除法,您需要执行多少次迭代?这是单字除法算法的基本情况。此外,在双精度中使用 Newton Raphson 进行 64 位除数,我们再次遇到舍入问题,因为 dp 尾数是 53 位长
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2015-08-30
  • 2022-01-18
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多