我已经看到了产生平方根的浮点位破解，如fast floating point square root 所示，但这种方法适用于浮点数。
是否有类似的方法可以在没有 32 位无符号整数循环的情况下找到整数平方根？我一直在网上搜索一个，但没有看到任何
（我的想法是纯二进制表示没有足够的信息来执行此操作，但由于它被限制为 32 位，因此我会猜不出来）
【问题讨论】：
uint32_t y = sqrt(x); 将通过转换为 double、取平方根并转换回来来完成这项工作。
标签： c unsigned-integer square-root
这是@njuffa 给出的代码的变体，可能会引起人们的兴趣。它本质上是无分支的，并且避免了最初猜测的查找表，尽管它确实做了三个可能昂贵的除法而不是两个。 （但请参阅下面的变体，它用查找表代替第一个除法。）
#include <stdint.h>
uint8_t clz (uint32_t a); // count leading zeros
uint32_t umul_16_16 (uint16_t a, uint16_t b); // 16x16 bit multiply
uint16_t udiv_32_16 (uint32_t x, uint16_t y); // 32/16 bit division
uint16_t my_isqrt (uint32_t x)
{
uint16_t lz, y;
lz = clz(x); // use clz(x | 1) if clz(0) is undefined
x <<= lz & 30u; // 2**30 <= x < 2**32 (or x == 0)
y = 1u + (x >> 30); // 2 <= y <= 4
y = (y << 1) + udiv_32_16(x >> 27, y); // 8 <= y <= 16
y = (y << 3) + udiv_32_16(x >> 21, y); // 128 <= y <= 256
y = (y << 7) + udiv_32_16(x >> 9, y); // 32768 <= y < 65536
return (y - (umul_16_16(y, y) > x)) >> (lz >> 1);
}
一些注意事项：
clz、umul_16_16 和 udiv_32_16 的定义与 @njuffa 的帖子中的完全相同x 阳性结果，但也恰好为x = 0 产生正确结果。clz(0) 是有效的，并返回31 或32（@njuffa 给出的实现也是如此）。某些架构可能会提供clz 指令，clz(0) 会为此提供未定义的行为。在这种情况下，请改用clz(x | 1)。y <= 65536；证明y < 65536，使结果不会溢出uint16_t，需要逐案分析该算法与 Python 的 math.isqrt 函数中使用的算法相同，并在从此处开始的 cmets 中进行了描述：https://github.com/python/cpython/blob/0b58bac3e7877d722bdbd3c38913dba2cb212f13/Modules/mathmodule.c#L1577
算法描述，基于牛顿法，仔细控制误差项：
x 缩放为 4 的幂，以使输入具有零或一个前导零。y 的初始值是x >> 28 平方根的上限或下限
y << 2 是x >> 24 的平方根的近似值。下一行应用牛顿法的一个步骤来获得更好的近似值，这再次可以证明是x >> 24 平方根的上限或下限。y << 4 和 x >> 16。y << 8 和 x。在这第三个也是最后一个牛顿步骤之后，y 是x 平方根的上限或下限。而且，它永远不等于65536（这个证明有点尴尬），所以它仍然适合uint16_t。y 是否大于x 的真实平方根，如果是，则在返回之前进行调整以补偿x 的原始缩放比例。最后，这是变体的变体，它用查找表代替了牛顿的第一步。和以前一样，第一次除法可以作为 16 位乘 16 位运算执行，如果在目标架构中这样更快的话。
#include <stdint.h>
uint8_t clz (uint32_t a); // count leading zeros
uint32_t umul_16_16 (uint16_t a, uint16_t b); // 16x16 bit multiply
uint16_t udiv_32_16 (uint32_t x, uint16_t y); // 32/16 bit division
static const uint8_t sqrt_tab[32] =
{
16, 24, 32, 40, 48, 56, 64, 72, 64, 64, 72, 72, 80, 80, 88, 88,
88, 88, 96, 96, 96, 104, 104, 104, 112, 112, 112, 112, 120, 120, 120, 120,
};
uint16_t my_isqrt (uint32_t x)
{
uint16_t lz, y;
lz = clz(x); // use clz(x | 1) if clz(0) is undefined
x <<= lz & 30u; // 2**30 <= x < 2**32 (or x == 0)
y = sqrt_tab[x >> 27]; // 64 <= y < 128
y = y + udiv_32_16(x >> 18, y); // 128 <= y <= 256
y = (y << 7) + udiv_32_16(x >> 9, y); // 32768 <= y < 65536
return (y - (umul_16_16(y, y) > x)) >> (lz >> 1);
}
【讨论】：
这个答案假设目标平台不支持浮点，或者非常慢的浮点支持（可能通过仿真）。
正如在 cmets 中所指出的，计数前导零 (CLZ) 指令可用于提供通过浮点操作数的指数部分提供的快速 log2 功能。 CLZ 也可以在不通过内部函数提供功能的平台上以合理的效率进行仿真，如下所示。
可以从查找表 (LUT) 中提取对几位有用的初始近似值，就像在浮点情况下一样，它可以通过牛顿迭代进一步细化。对于 32 位整数平方根，一到两次迭代通常就足够了。下面的 ISO-C99 代码显示了工作示例性实施，包括详尽的测试。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
uint8_t clz (uint32_t a); // count leading zeros
uint32_t umul_16_16 (uint16_t a, uint16_t b); // 16x16 bit multiply
uint16_t udiv_32_16 (uint32_t x, uint16_t y); // 32/16 bit division
/* LUT for initial square root approximation */
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
};
/* table lookup for initial guess followed by division-based Newton iteration */
uint16_t my_isqrt (uint32_t x)
{
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; // use for overflow check on divisions
// 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 (umul_16_16 (y, y) > x) y--; // adjust quotient if too large
return y; // (uint16_t)sqrt((double)x)
}
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
};
/* count leading zeros (for non-zero argument); a machine instruction on many architectures */
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];
}
/* 16x16->32 bit unsigned multiply; machine instruction on many architectures */
uint32_t umul_16_16 (uint16_t a, uint16_t b)
{
return (uint32_t)a * b;
}
/* 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;
}
int main (void)
{
uint32_t x;
uint16_t res, ref;
printf ("testing 32-bit integer square root\n");
x = 0;
do {
ref = (uint16_t)sqrt((double)x);
res = my_isqrt (x);
if (res != ref) {
printf ("error: x=%08x res=%08x ref=%08x\n", x, res, ref);
printf ("exhaustive test FAILED\n");
return EXIT_FAILURE;
}
x++;
} while (x);
printf ("exhaustive test PASSED\n");
return EXIT_SUCCESS;
}
【讨论】：
没有。您需要在某处引入日志；快速浮点平方根之所以有效，是因为位表示中的对数。
最快的方法可能是 n -> floor(sqrt(n)) 的查找表。您不会将所有值存储在表中，而只存储平方根更改的值。使用二分查找以 log(n) 时间在表中查找结果。
【讨论】：