【问题标题】:Optimizing Fixed-Point Sqrt优化定点 Sqrt
【发布时间】:2023-04-09 21:47:01
【问题描述】:

我做了我认为不错的fixed-point square root algorithm

template<int64_t M, int64_t P>
typename enable_if<M + P == 32, FixedPoint<M, P>>::type sqrt(FixedPoint<M, P> f)
{
    if (f.num == 0)
        return 0;
    //Reduce it to the 1/2 to 2 range (based around FixedPoint<2, 30> to avoid left/right shift branching)
    int64_t num{ f.num }, faux_half{ 1 << 29 };
    ptrdiff_t mag{ 0 };
    while (num < (faux_half)) {
        num <<= 2;
        ++mag;
    }

    int64_t res = (M % 2 == 0 ? SQRT_32_EVEN_LOOKUP : SQRT_32_ODD_LOOKUP)[(num >> (30 - 4)) - (1LL << 3)];
    res >>= M / 2 + mag - 1; //Finish making an excellent guess
    for (int i = 0; i < 2; ++i)
        //                            \    |   /
        //                             \   |  /
        //                              _| V L
        res = (res + (int64_t(f.num) << P) / res) >> 1; //Use Newton's method to improve greatly on guess
        //                               7 A r
        //                              /  |  \
        //                             /   |   \
        //                       The Infamous Time Eater
    return FixedPoint<M, P>(res, true);
}

但是,在分析之后(在发布模式下),我发现这里的除法占用了该算法所花费时间的 83%。我可以通过用乘法代替除法来加速 6 倍,但这是错误的。不幸的是,我发现整数除法比乘法慢得多。有什么办法可以优化吗?

如果需要这张表。

const array<int32_t, 24> SQRT_32_EVEN_LOOKUP = {
     0x2d413ccd, //magic numbers calculated by taking 0.5 + 0.5 * i / 8 from i = 0 to 23, multiplying by 2^30, and converting to hex
     0x30000000,
     0x3298b076,
     0x3510e528,
     0x376cf5d1,
     0x39b05689,
     0x3bddd423,
     0x3df7bd63,
     0x40000000,
     0x41f83d9b,
     0x43e1db33,
     0x45be0cd2,
     0x478dde6e,
     0x49523ae4,
     0x4b0bf165,
     0x4cbbb9d6,
     0x4e623850,
     0x50000000,
     0x5195957c,
     0x532370b9,
     0x54a9fea7,
     0x5629a293,
     0x57a2b749,
     0x59159016
};

SQRT_32_ODD_LOOKUP 就是SQRT_32_EVEN_LOOKUP 除以sqrt(2)

【问题讨论】:

标签: c++ algorithm performance fixed-point sqrt


【解决方案1】:

重新发明轮子,真的,而且不是很好。正确的解决方案是使用 NR 计算 1/sqrt(x),然后乘以一次以得到 x/sqrt(x) - 只需预先检查 x==0

之所以这么好,是因为 y=1/sqrt(x) 的 NR 步骤只是 y = (3y-x*y*y)*y/2。这都是简单的乘法。

【讨论】:

  • 我认为应该是y = (3 - x*y*y)*y/2
  • 感谢您的提示。我最初认为这行不通,因为平方根倒数会在 1 的另一边(对FixedPoint&lt;2, 30&gt;FixedPoint&lt;32, 0&gt; 非常不利)。但是,使用 y = (3 - x*y*y/C^2)*y/2 抵消平方根的逆以使 y 在 2^29 和 2^31 之间有效。它比旧算法准确度稍差,但速度快了 3 倍。
  • @user155698:您可以使用第 4 个 NR 步骤,但仍然快 2 倍。
  • @chmike:如果你在定点上做(3y-x*y*y)*y/2,在定点上做x*(1/sqrt(x)),那么结果就是在定点上sqrt(x)
  • @MSalters 我查找了如何计算 1/sqrt(x),但地震技巧(近似值)使用浮点数。 (3y-x*y*y)*y/2 产生什么?那是 1/sqrt(x) 吗?我很想知道它是否比不使用乘法和除法的my algorithm 更快。我来比较一下。
猜你喜欢
  • 1970-01-01
  • 2016-09-04
  • 1970-01-01
  • 2011-06-07
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2020-08-11
  • 1970-01-01
相关资源
最近更新 更多