固定点下的倒数平方根

Question

我正在寻找适用于定点16.16数字的最佳平方根倒数算法。以下是我目前拥有的代码（但基本上它取平方根并除以原始数字，而我想获得不需要除法的平方根倒数）。如果有任何变化，该代码将编译为armv5te。

uint32_t INVSQRT(uint32_t n)
{
    uint64_t op, res, one;
    op = ((uint64_t)n<<16);
    res = 0;
    one = (uint64_t)1 << 46;
    while (one > op) one >>= 2;
    while (one != 0)
    {
        if (op >= res + one)
        {
            op -= (res + one);
            res +=  (one<<1);
        }
        res >>= 1;
        one >>= 2;
    }
    res<<=16;
    res /= n;
    return(res);
}

- Jonathan

---

啰嗦：您可能是指“倒数”的平方根吧？ - Oliver Charlesworth

请参考以下链接：http://zh.wikipedia.org/wiki/Fast_inverse_square_root - Guerrero

正好我就要回复这个。 - Jonathan

@Guerrero，@Jonathan：是的，那里的名称有些误导（事实上，该文章说“快速求平方根的倒数...是一种计算平方根倒数的方法”）。求平方根的倒数只是平方！ - Oliver Charlesworth

Answer 1

关键是将牛顿法应用于问题x - 1/y^2 = 0。因此，给定x，使用迭代方案解决y。

Y_(n+1) = y_n * (3 - x*y_n^2)/2

除以2只是一个位移，或者说最坏的情况下是乘以0.5。这个方案会收敛到y=1/sqrt(x)，完全满足要求，并且没有真正的除法。

唯一的问题是需要一个合理的初始值来计算y。据我所知，对于迭代收敛的估计y有限制。

- user85109

---

你可以通过使用指数搜索来找到一个区间，该区间内的xyy - 1会出现符号变化，然后在该区间上使用割线法，在此之后使用牛顿法作为起点。 - Michael Anderson

但是在代码中会是什么样子呢？你认为它比我已经有的更有效吗？ - Jonathan

这段代码会使每次迭代中正确数字的数量翻倍。因此，如果您正在处理16位数字，则需要大约4次迭代才能收敛，如果您在第一次逼近中从1个正确数字开始。这段代码应该如何编写呢？它看起来非常像我上面写的那一行。 - user85109

Answer 2

ARMv5TE处理器提供了快速的整数乘法器和“计算前导零”指令。它们通常带有适度大小的缓存。基于此，实现高性能最合适的方法似乎是使用表查找进行初始近似，然后使用两个牛顿拉弗森迭代来实现完全精确的结果。我们可以通过将额外的预计算纳入表中来进一步加速这些迭代，这是Cray计算机四十年前使用的一种技术。

下面的函数fxrsqrt()实现了这种方法。它以8位近似值r为开始，用每个表元素存储3r（在32位条目的低10位中）和r的立方（在32位条目的上22位中），而不是存储r。这允许快速计算第一次迭代，如下所示：r₁ = 0.5 * (3 * r - a * r³)。然后按照常规方式计算第二次迭代，r₂ = 0.5 * r₁ * (3 - r₁ * (r₁ * a))。

为了能够准确地执行这些计算，无论输入的大小如何，参数a在计算开始时进行了归一化，实质上将其表示为与2^scal乘以比例因子的2.32定点数。在计算结束时，根据公式1/sqrt(2²ⁿ) = 2^-n对结果进行非标准化。通过将最高位被舍弃的位为1的结果四舍五入，提高了精度，从而几乎所有结果都被正确地舍入。详尽的测试报告如下：results too low: 639 too high: 1454 not correctly rounded: 2093。

该代码使用了两个辅助函数：__clz()确定一个非零32位参数中前导0位的数量。__umulhi()计算两个无符号32位整数的完整64位乘积的32位最高位。这两个函数应通过编译器内置函数或使用一些内联汇编来实现。在下面的代码中，我展示了适用于ARM CPU的可移植实现以及适用于x86平台的内联汇编版本。在ARMv5TE平台上，__clz()应映射到CLZ指令，__umulhi()应映射到UMULL。

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>

#define USE_OWN_INTRINSICS 1

#if USE_OWN_INTRINSICS
__forceinline int __clz (uint32_t a)
{
    int r;
    __asm__ ("bsrl %1,%0\n\t" : "=r"(r): "r"(a));
    return 31 - r;
}

uint32_t __umulhi (uint32_t a, uint32_t b)
{
    uint32_t r;
    __asm__ ("movl %1,%%eax\n\tmull %2\n\tmovl %%edx,%0\n\t"
             : "=r"(r) : "r"(a), "r"(b) : "eax", "edx");
    return r;
}
#else // USE_OWN_INTRINSICS
int __clz (uint32_t a)
{
    uint32_t r = 32;
    if (a >= 0x00010000) { a >>= 16; r -= 16; }
    if (a >= 0x00000100) { a >>=  8; r -=  8; }
    if (a >= 0x00000010) { a >>=  4; r -=  4; }
    if (a >= 0x00000004) { a >>=  2; r -=  2; }
    r -= a - (a & (a >> 1));
    return r;
}

uint32_t __umulhi (uint32_t a, uint32_t b)
{
    return (uint32_t)(((uint64_t)a * b) >> 32);
}
#endif // USE_OWN_INTRINSICS

/*
 * For each sub-interval in [1, 4), use an 8-bit approximation r to reciprocal
 * square root. To speed up subsequent Newton-Raphson iterations, each entry in
 * the table combines two pieces of information: The least-significant 10 bits
 * store 3*r, the most-significant 22 bits store r**3, rounded from 24 down to
 * 22 bits such that accuracy is optimized.
 */
uint32_t rsqrt_tab [96] = 
{
    0xfa0bdefa, 0xee6af6ee, 0xe5effae5, 0xdaf27ad9,
    0xd2eff6d0, 0xc890aec4, 0xc10366bb, 0xb9a71ab2,
    0xb4da2eac, 0xadce7ea3, 0xa6f2b29a, 0xa279a694,
    0x9beb568b, 0x97a5c685, 0x9163027c, 0x8d4fd276,
    0x89501e70, 0x8563da6a, 0x818ac664, 0x7dc4fe5e,
    0x7a122258, 0x7671be52, 0x72e44a4c, 0x6f68fa46,
    0x6db22a43, 0x6a52623d, 0x67041a37, 0x65639634,
    0x622ffe2e, 0x609cba2b, 0x5d837e25, 0x5bfcfe22,
    0x58fd461c, 0x57838619, 0x560e1216, 0x53300a10,
    0x51c72e0d, 0x50621a0a, 0x4da48204, 0x4c4c2e01,
    0x4af789fe, 0x49a689fb, 0x485a11f8, 0x4710f9f5,
    0x45cc2df2, 0x448b4def, 0x421505e9, 0x40df5de6,
    0x3fadc5e3, 0x3e7fe1e0, 0x3d55c9dd, 0x3d55d9dd,
    0x3c2f41da, 0x39edd9d4, 0x39edc1d4, 0x38d281d1,
    0x37bae1ce, 0x36a6c1cb, 0x3595d5c8, 0x3488f1c5,
    0x3488fdc5, 0x337fbdc2, 0x3279ddbf, 0x317749bc,
    0x307831b9, 0x307879b9, 0x2f7d01b6, 0x2e84ddb3,
    0x2d9005b0, 0x2d9015b0, 0x2c9ec1ad, 0x2bb0a1aa,
    0x2bb0f5aa, 0x2ac615a7, 0x29ded1a4, 0x29dec9a4,
    0x28fabda1, 0x2819e99e, 0x2819ed9e, 0x273c3d9b,
    0x273c359b, 0x2661dd98, 0x258ad195, 0x258af195,
    0x24b71192, 0x24b6b192, 0x23e6058f, 0x2318118c,
    0x2318718c, 0x224da189, 0x224dd989, 0x21860d86,
    0x21862586, 0x20c19183, 0x20c1b183, 0x20001580
};

/* This function computes the reciprocal square root of its 16.16 fixed-point 
 * argument. After normalization of the argument if uses the most significant
 * bits of the argument for a table lookup to obtain an initial approximation 
 * accurate to 8 bits. This is followed by two Newton-Raphson iterations with
 * quadratic convergence. Finally, the result is denormalized and some simple
 * rounding is applied to maximize accuracy.
 *
 * To speed up the first NR iteration, for the initial 8-bit approximation r0
 * the lookup table supplies 3*r0 along with r0**3. A first iteration computes
 * a refined estimate r1 = 1.5 * r0 - x * r0**3. The second iteration computes
 * the final result as r2 = 0.5 * r1 * (3 - r1 * (r1 * x)).
 *
 * The accuracy for all arguments in [0x00000001, 0xffffffff] is as follows: 
 * 639 results are too small by one ulp, 1454 results are too big by one ulp.
 * A total of 2093 results deviate from the correctly rounded result.
 */
uint32_t fxrsqrt (uint32_t a)
{
    uint32_t s, r, t, scal;

    /* handle special case of zero input */
    if (a == 0) return ~a;
    /* normalize argument */
    scal = __clz (a) & 0xfffffffe;
    a = a << scal;
    /* initial approximation */
    t = rsqrt_tab [(a >> 25) - 32];
    /* first NR iteration */
    r = (t << 22) - __umulhi (t, a);
    /* second NR iteration */
    s = __umulhi (r, a);
    s = 0x30000000 - __umulhi (r, s);
    r = __umulhi (r, s);
    /* denormalize and round result */
    r = ((r >> (18 - (scal >> 1))) + 1) >> 1;
    return r;
}

/* reference implementation, 16.16 reciprocal square root of non-zero argment */
uint32_t ref_fxrsqrt (uint32_t a)
{
    double arg = a / 65536.0;
    double rsq = sqrt (1.0 / arg);
    uint32_t r = (uint32_t)(rsq * 65536.0 + 0.5);
    return r;
}

int main (void)
{
    uint32_t arg = 0x00000001;
    uint32_t res, ref;
    uint32_t err, lo = 0, hi = 0;

    do {
        res = fxrsqrt (arg);
        ref = ref_fxrsqrt (arg);

        err = 0;
        if (res < ref) {
            err = ref - res;
            lo++;
        }
        if (res > ref) {
            err = res - ref;
            hi++;
        }
        if (err > 1) {
            printf ("!!!! arg=%08x  res=%08x  ref=%08x\n", arg, res, ref);
            return EXIT_FAILURE;
        }
        arg++;
    } while (arg);
    printf ("results too low: %u  too high: %u  not correctly rounded: %u\n", 
            lo, hi, lo + hi);
    return EXIT_SUCCESS;
}

- njuffa

Answer 3

我有一个方案，称之为“32位定点数的快速逆平方根”。没有表格，没有参考资料，直接用好的猜测。如果您想要，请跳到下面的源代码，但要注意几件事情。

1. (x * y)>>16可以替换为任何您想要的定点乘法方案。
2. 这不需要64位[长字]，我只是为了演示的便利而使用它。长字用于防止乘法溢出。定点数学库将具有更好处理此问题的定点乘法函数。
3. 初始猜测非常好，因此您在第一次尝试中得到相对精确的结果。
4. 代码比演示所需的要冗长。
5. 小于65536（＜1）和大于32767<<16的值不能使用。
6. 如果您的硬件具有除法功能，则通常不比使用平方根表和除法更快。如果没有，这将避免除法。

int fxisqrt(int input){

    if(input <= 65536){
        return 1;
    }

    long xSR = input>>1;
    long pushRight = input;
    long msb = 0;
    long shoffset = 0;
    long yIsqr = 0;
    long ysqr = 0;
    long fctrl = 0;
    long subthreehalf = 0;

    while(pushRight >= 65536){
        pushRight >>=1;
        msb++;
    }

    shoffset = (16 - ((msb)>>1));
    yIsqr = 1<<shoffset;
    //y = (y * (98304 - ( ( (x>>1) * ((y * y)>>16 ) )>>16 ) ) )>>16;   x2
    //Incantation 1
    ysqr = (yIsqr * yIsqr)>>16;
    fctrl = (xSR * ysqr)>>16;
    subthreehalf = 98304 - fctrl;
    yIsqr = (yIsqr * subthreehalf)>>16;
    //Incantation 2 - Increases precision greatly, but may not be neccessary
    ysqr = (yIsqr * yIsqr)>>16;
    fctrl = (xSR * ysqr)>>16;
    subthreehalf = 98304 - fctrl;
    yIsqr = (yIsqr * subthreehalf)>>16;
    return yIsqr;
}

- ponut64