最快的整数平方根[指令数量最少]

5 回答

• 1

  我不知道如何把它变成一个有效的算法但是当我在80年代研究这个时，出现了一个有趣的模式 . 当舍入平方根时，平方根有两个整数而不是前一个整数（零之后） .

  因此，一个数字（零）的平方根为零，两个的平方根为1（1和2），4的平方根为2（3,4,5和6），依此类推 . 可能不是一个有用的答案，但有趣的是 .

  回复于 2024-04-27T11:15:23+08:00
• -1

  看看here .

  例如，在3（a）处有这种方法，它可以很容易地适应64-> 32位平方根，并且可以简单地转换为汇编程序：

  /* by Jim Ulery */
  static unsigned julery_isqrt(unsigned long val) {
      unsigned long temp, g=0, b = 0x8000, bshft = 15;
      do {
          if (val >= (temp = (((g << 1) + b)<<bshft--))) {
             g += b;
             val -= temp;
          }
      } while (b >>= 1);
      return g;
  }
  

  没有划分，没有乘法，只有位移 . 但是，所花费的时间在某种程度上是不可预测的，特别是如果您使用分支（在ARM RISC条件指令上可行） .

  通常，this page列出了计算平方根的方法 . 如果你碰巧想要生成一个快速的反平方根（即 x**(-0.5) ），或者只是对优化代码的惊人方法感兴趣，请看一下this，this和this .

  回复于 2024-04-27T11:15:23+08:00
• 7

  这与您的相同，但操作次数较少 . （我在你的代码循环中计算了9个操作，包括for循环中的测试和增量 i 以及3个赋值，但是当在ASM中编码时，其中一些可能会消失？下面的代码中有6个操作，如果你将 g*g>n 算作两个（没有任务）） .

  int isqrt(int n) {
    int g = 0x8000;
    int c = 0x8000;
    for (;;) {
      if (g*g > n) {
        g ^= c;
      }
      c >>= 1;
      if (c == 0) {
        return g;
      }
      g |= c;
    }
  }
  

  我明白了here . 如果您展开循环并根据输入中的最高非零位跳转到适当的位置，则可以消除比较 .

  Update

  我一直在考虑使用牛顿方法 . 理论上，每次迭代的精度位数应该加倍 . 这意味着当答案中的正确位数很少时，牛顿的方法比任何其他建议都要糟糕得多;但是，在答案中存在大量正确位的情况下，情况会发生变化 . 考虑到大多数建议似乎每位需要4个周期，这意味着牛顿方法的一次迭代（除了超过4位，除非为移位= 18个周期的加法1的16个周期用于除法1） .

  所以，我的建议是通过建议的方法之一 Build 8位答案（8 * 4 = 32个周期），然后执行牛顿方法的一次迭代（18个周期），将位数加倍到16个 . 这是一个总数50个循环（加上可能额外的4个循环，在应用牛顿方法之前获得9个位，加上可能有2个循环来克服牛顿方法偶尔遇到的过冲） . 这是最多56个周期，据我所知，可以看到任何其他建议 .

  Second Update

  我编写了混合算法的想法 . 牛顿方法本身没有开销;你只需申请并加倍有效数字 . 在应用牛顿方法之前，问题是要有可预测的有效位数 . 为此，我们需要弄清楚答案的最重要部分将出现在哪里 . 使用另一张海报给出的快速DeBruijn序列方法的修改，我可以在估计中在大约12个周期中执行该计算 . 另一方面，知道答案msb的位置可以加速所有方法（平均而非最坏的情况），所以无论如何它似乎都值得 .

  在计算了答案的msb之后，我运行了上面建议的一些算法，然后用一两轮Newton方法完成它 . 我们通过以下计算决定何时运行牛顿方法：根据评论中的计算，答案的一位需要大约8个周期;一轮牛顿方法需要大约18个周期（除法，加法和移位，也许是赋值），所以我们应该只运行牛顿方法，如果我们要从中得到至少三位 . 因此对于6位答案，我们可以运行线性方法3次得到3个有效位，然后运行牛顿方法1次得到另外3个 . 对于15位答案，我们运行线性方法4次得到4位，然后是牛顿的方法两次得到另一个4然后另一个7.依此类推 .

  这些计算取决于通过线性方法确切知道需要多少个周期，而牛顿方法需要多少个周期 . 如果“经济学”发生变化，例如通过发现以线性方式构建比特的更快方式，何时调用牛顿方法的决定将会改变 .

  我展开循环并将决策实现为开关，我希望这将转换为汇编中的快速表查找 . 我并不完全确定在每种情况下我都有最小的循环次数，因此可能需要进一步调整 . 例如，对于s = 10，您可以尝试获得5位然后应用牛顿方法一次而不是两次 .

  在某些情况下，我愿意接受稍微不正确的答案 . 在应用牛顿方法来纠正与 m^2-1 形式的数字一起出现的一个一个错误之后，至少使用两个循环 . 并且一个循环用于在开始时测试输入0，因为算法可以永远不会取零的平方根，您可以消除该测试 . 最后，如果您在答案中只需要8个有效位，则可以删除其中一个牛顿方法计算 .

  #include <inttypes.h>
  #include <stdint.h>
  #include <stdbool.h>
  #include <stdio.h>
  
  uint32_t isqrt1(uint32_t n);
  
  int main() {
    uint32_t n;
    bool it_works = true;
    for (n = 0; n < UINT32_MAX; ++n) {
      uint32_t sr = isqrt1(n);
      if ( sr*sr > n || ( sr < 65535 && (sr+1)*(sr+1) <= n )) {
        it_works = false;
        printf("isqrt(%" PRIu32 ") = %" PRIu32 "\n", n, sr);
      }
    }
    if (it_works) {
      printf("it works\n");
    }
    return 0;
  }
  
  /* table modified to return shift s to move 1 to msb of square root of x */
  /*
  static const uint8_t debruijn32[32] = {
      0, 31, 9, 30, 3,  8, 13, 29,  2,  5,  7, 21, 12, 24, 28, 19,
      1, 10, 4, 14, 6, 22, 25, 20, 11, 15, 23, 26, 16, 27, 17, 18
  };
  */
  
  static const uint8_t debruijn32[32] = {
    15,  0, 11, 0, 14, 11, 9, 1, 14, 13, 12, 5, 9, 3, 1, 6,
    15, 10, 13, 8, 12,  4, 3, 5, 10,  8,  4, 2, 7, 2, 7, 6
  };
  
  /* based on CLZ emulation for non-zero arguments, from
   * http://stackoverflow.com/questions/23856596/counting-leading-zeros-in-a-32-bit-unsigned-integer-with-best-algorithm-in-c-pro
   */
  uint8_t shift_for_msb_of_sqrt(uint32_t x) {
    x |= x >>  1;
    x |= x >>  2;
    x |= x >>  4;
    x |= x >>  8;
    x |= x >> 16;
    x++;
    return debruijn32 [x * 0x076be629 >> 27];
  }
  
  uint32_t isqrt1(uint32_t n) {
    if (n==0) return 0;
  
    uint32_t s = shift_for_msb_of_sqrt(n);
    uint32_t c = 1 << s;
    uint32_t g = c;
  
    switch (s) {
    case 9:
    case 5:
      if (g*g > n) {
        g ^= c;
      }
      c >>= 1;
      g |= c;
    case 15:
    case 14:
    case 13:
    case 8:
    case 7:
    case 4:
      if (g*g > n) {
        g ^= c;
      }
      c >>= 1;
      g |= c;
    case 12:
    case 11:
    case 10:
    case 6:
    case 3:
      if (g*g > n) {
        g ^= c;
      }
      c >>= 1;
      g |= c;
    case 2:
      if (g*g > n) {
        g ^= c;
      }
      c >>= 1;
      g |= c;
    case 1:
      if (g*g > n) {
        g ^= c;
      }
      c >>= 1;
      g |= c;
    case 0:
      if (g*g > n) {
        g ^= c;
      }
    }
  
    /* now apply one or two rounds of Newton's method */
    switch (s) {
    case 15:
    case 14:
    case 13:
    case 12:
    case 11:
    case 10:
      g = (g + n/g) >> 1;
    case 9:
    case 8:
    case 7:
    case 6:
      g = (g + n/g) >> 1;
    }
  
    /* correct potential error at m^2-1 for Newton's method */
    return (g==65536 || g*g>n) ? g-1 : g;
  }
  

  在我的机器上进行轻度测试（这无疑与你的一样），新的 isqrt1 例程平均比我之前提供的 isqrt 例程快40％ .

  回复于 2024-04-27T11:15:23+08:00
• 5

  如果乘法是相同的速度（或快于！）加法和移位，或者如果缺少快速按容量移位的寄存器指令，则以下内容将没有帮助 . 除此以外：

  你在每个循环周期重新计算 temp*temp ，但 temp = res | add ，它与 res + add 相同，因为它们的位不重叠，并且（a）你已经在前一个循环周期计算了 res*res ，并且（b） add 有特殊的结构（它总是只有一点） . 所以，使用 (a+b)^2 = a^2 + 2ab + b^2 ，你已经 a^2 和 b^2 的事实只是向左移动了两位，比单位 b 向左移动了两次，而 2ab 只是 a 左移了1位多于位置在 b 中的单个位，你可以摆脱乘法：

  unsigned short int int_sqrt32(unsigned int x)
  {
      unsigned short int res = 0;
      unsigned int res2 = 0;
      unsigned short int add = 0x8000;   
      unsigned int add2 = 0x80000000;   
      int i;
      for(i = 0; i < 16; i++)
      {
          unsigned int g2 = res2 + (res << i) + add2;
          if (x >= g2)
          {
              res |= add;
              res2 = g2;
          }
          add >>= 1;
          add2 >>= 2;
      }
      return res;
  }
  

  另外我猜想对所有变量使用相同的类型（ unsigned int ）是个更好的主意，因为根据C标准，所有算术都需要在执行算术运算之前将较窄的整数类型提升（转换）为最广泛的类型，如有必要，随后进行后向转换 . （这当然可以通过足够智能的编译器进行优化，但为什么要承担风险呢？）

  回复于 2024-04-27T11:15:23+08:00
• 2

  从注释轨迹来看，似乎RISC处理器仅提供32x32-> 32位乘法和16x16-> 32位乘法 . 32x-32-> 64位加宽乘法，或 MULHI 指令返回64位的高32位没有提供产品 .

  这似乎排除了基于Newton-Raphson迭代的方法，这可能是低效的，因为它们通常需要 MULHI 指令或加宽乘法用于中间定点算术 .

  下面的C99代码使用不同的迭代方法，只需要16x16-> 32位乘法，但稍微线性收敛，需要多达六次迭代 . 此方法需要 CLZ 功能来快速确定迭代的起始猜测 . Asker在评论中指出，所使用的RISC处理器不提供CLZ功能 . 因此需要仿真CLZ，并且由于仿真会增加存储和指令计数，这可能会使这种方法失去竞争力 . 我执行了一个强力搜索来确定具有最小乘数的deBruijn查找表 .

  这种迭代算法提供了非常接近所需结果的原始结果，即 (int)sqrt(x) ，但由于整数算术的截断性质，总是有些偏高 . 为了得到最终结果，结果有条件地减少，直到结果的平方小于或等于原始参数 .

  在代码中使用 volatile 限定符仅用于确定所有命名变量实际上可以分配为16位数据而不会影响功能 . 我不知道这是否有任何优势，但注意到OP在其代码中专门使用了16位变量 .

  请注意，对于大多数处理器，最后的更正步骤不应涉及任何分支 . 产品 y*y 可以通过结转（或借出）从 x 中减去，然后 y 通过减值进位（或借入）进行修正 . 所以每一步都应该是一个序列 MUL ， SUBcc ， SUBC .

  因为循环实现迭代会产生大量开销，所以我选择完全展开循环，但提供两个早期检查 . 手动计算操作我计算最快情况下的46次操作，平均情况下54次操作，最坏情况下60次操作 .

  #include <stdio.h>
  #include <stdlib.h>
  #include <stdint.h>
  #include <math.h>
  
  static const uint8_t clz_tab[32] = {
      31, 22, 30, 21, 18, 10, 29,  2, 20, 17, 15, 13, 9,  6, 28, 1,
      23, 19, 11,  3, 16, 14,  7, 24, 12,  4,  8, 25, 5, 26, 27, 0};
  
  uint8_t clz (uint32_t a)
  {
      a |= a >> 16;
      a |= a >> 8;
      a |= a >> 4;
      a |= a >> 2;
      a |= a >> 1;
      return clz_tab [0x07c4acdd * a >> 27];
  }
  
  /* 16 x 16 -> 32 bit unsigned multiplication; should be single instruction */
  uint32_t umul16w (uint16_t a, uint16_t b)
  {
      return (uint32_t)a * b;
  }
  
  /* Reza Hashemian, "Square Rooting Algorithms for Integer and Floating-Point
     Numbers", IEEE Transactions on Computers, Vol. 39, No. 8, Aug. 1990, p. 1025
  */
  uint16_t isqrt (uint32_t x)
  {
      volatile uint16_t y, z, lsb, mpo, mmo, lz, t;
  
      if (x == 0) return x; // early out, code below can't handle zero
  
      lz = clz (x);         // #leading zeros, 32-lz = #bits of argument
      lsb = lz & 1;
      mpo = 17 - (lz >> 1); // m+1, result has roughly half the #bits of argument
      mmo = mpo - 2;        // m-1
      t = 1 << mmo;         // power of two for two's complement of initial guess
      y = t - (x >> (mpo - lsb)); // initial guess for sqrt
      t = t + t;            // power of two for two's complement of result
      z = y;
  
      y = (umul16w (y, y) >> mpo) + z;
      y = (umul16w (y, y) >> mpo) + z;
      if (x >= 0x40400) {
          y = (umul16w (y, y) >> mpo) + z;
          y = (umul16w (y, y) >> mpo) + z;
          if (x >= 0x1002000) {
              y = (umul16w (y, y) >> mpo) + z;
              y = (umul16w (y, y) >> mpo) + z;
          }
      }
  
      y = t - y; // raw result is 2's complement of iterated solution
      y = y - umul16w (lsb, (umul16w (y, 19195) >> 16)); // mult. by sqrt(0.5) 
  
      if ((int32_t)(x - umul16w (y, y)) < 0) y--; // iteration may overestimate 
      if ((int32_t)(x - umul16w (y, y)) < 0) y--; // result, adjust downward if 
      if ((int32_t)(x - umul16w (y, y)) < 0) y--; // necessary 
  
      return y; // (int)sqrt(x)
  }
  
  int main (void)
  {
      uint32_t x = 0;
      uint16_t res, ref;
  
      do {
          ref = (uint16_t)sqrt((double)x);
          res = isqrt (x);
          if (res != ref) {
              printf ("!!!! x=%08x  res=%08x  ref=%08x\n", x, res, ref);
              return EXIT_FAILURE;
          }
          x++;
      } while (x);
      return EXIT_SUCCESS;
  }
  

  另一种可能性是使用牛顿迭代作为平方根，尽管划分成本很高 . 对于小输入，只需要一次迭代 . 虽然提问者没有说明这一点，但基于 DIV 操作的16个周期的执行时间，我强烈怀疑这实际上是 32/16->16 位除法，需要额外的保护代码以避免溢出，定义为不适合16的商 . 位 . 我已根据此假设为我的代码添加了适当的安全措施 .

  由于Newton迭代每次应用时都会将好位数加倍，因此我们只需要一个低精度的初始猜测，可以根据参数的五个前导位从表中轻松检索 . 为了获取这些，我们首先将参数标准化为2.30定点格式，附加隐式比例因子为232-（lz＆~1），其中 lz 是参数中前导零的数量 . 与前一种方法一样，迭代并不总能提供准确的结果，因此如果初步结果太大，则必须应用修正 . 我计算快速路径的49个周期，慢速路径的70个周期（平均60个周期） .

  static const uint16_t sqrt_tab[32] = 
  { 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
    0x85ff, 0x8cff, 0x94ff, 0x9aff, 0xa1ff, 0xa7ff, 0xadff, 0xb3ff,
    0xb9ff, 0xbeff, 0xc4ff, 0xc9ff, 0xceff, 0xd3ff, 0xd8ff, 0xdcff, 
    0xe1ff, 0xe6ff, 0xeaff, 0xeeff, 0xf3ff, 0xf7ff, 0xfbff, 0xffff
  };
  
  /* 32/16->16 bit division. Note: Will overflow if x[31:16] >= y */
  uint16_t udiv_32_16 (uint32_t x, uint16_t y)
  {
      uint16_t r = x / y;
      return r;
  }
  
  /* table lookup for initial guess followed by division-based Newton iteration*/ 
  uint16_t isqrt (uint32_t x)
  {
      volatile uint16_t q, lz, y, i, xh;
  
      if (x == 0) return x; // early out, code below can't handle zero
  
      // initial guess based on leading 5 bits of argument normalized to 2.30
      lz = clz (x);
      i = ((x << (lz & ~1)) >> 27);
      y = sqrt_tab[i] >> (lz >> 1);
      xh = (x >> 16); // needed for overflow check on division
  
      // first Newton iteration, guard against overflow in division
      q = 0xffff;
      if (xh < y) q = udiv_32_16 (x, y);
      y = (q + y) >> 1;
      if (lz < 10) {
          // second Newton iteration, guard against overflow in division
          q = 0xffff;
          if (xh < y) q = udiv_32_16 (x, y);
          y = (q + y) >> 1;
      }
  
      if (umul16w (y, y) > x) y--; // adjust quotient if too large
  
      return y; // (int)sqrt(x)
  }
  

  回复于 2024-04-27T11:15:23+08:00