【算法】整数开平方算法之一——只使用逻辑运算和移位
# 定义**整数开平方**被广泛用于嵌入式的步进电机控制计算、游戏的寻路算法的内存池分配等多个领域。对比常见的**浮点数开平方**，整数开平方不需要依赖浮点数计算单元，且能在无硬件浮点数单元的嵌入式环境里实现高效率的开平方。
在[数论](https://zh.wikipedia.org/wiki/%E6%95%B0%E8%AE%BA)的定义里，对正整数n进行整数开方(`isqrt`)运算得出来的结果m，是小于或者等于n平方根的最大整数。
例如，`isqrt(27) = 5`因为`5 · 5 = 25 ≤ 27`且`6 · 6 = 36 ＞ 27`
# 一些整数开方算法思路
其实主要有两种：
* 牛顿开方法
* 逐二进制位开方法
就目前而言，根据我的实际测试，对于嵌入式环境，使用**逐二进制位开方法**是效率最高也是性能最高的。
## 牛顿开方法
一种既可以对整数进行开方又可以对浮点数求平方根的算法是牛顿开方法。这种算法的思路就是直接把浮点数开方的算法用到整数上，并且可行。
可以参考[这篇帖子](https://www.0xaa55.com/thread-1544-1-1.html)对牛顿开方迭代算法的描述。
牛顿法开方的整数开方的C语言实现如下（依赖`stdint.h`）
// 牛顿法整数开方
uint32_t int_sqrt(uint32_t s)
{
uint32_t x0 = s >> 1;
if (x0) // 参数检查
{
uint32_t x1 = ( x0 + s / x0 ) >> 1;
while ( x1 < x0 )
{
x0 = x1;
x1 = ( x0 + s / x0 ) >> 1;
}
return x0;
}
else
{
return s;
}
}
上述代码即可实现整数开平方根。但是它有一个缺点：**它依赖整数除法**。通常情况下，**整数除法的运算效率比浮点数除法低**，尤其是有符号整数的除法更是低。可以参考Intel文档对整数除法指令的周期数量和浮点数除法指令的周期数量对比来证明这一点。
虽说浮点数也可以使用这个算法，但实际上在对精度要求不高的场合，浮点数计算通常并不一定使用这种方法求平方根。而是使用[快速反平方根](https://en.wikipedia.org/wiki/Fast_inverse_square_root)算法先求出`1 ÷ 平方根`（也就是所谓的反平方根）的值，再用1去除以反平方根来得到平方根，或者根据使用的情况可以根本不需要除，比如做向量单位化运算的时候。
// 快速反平方根（并非整数开方）
float Q_rsqrt(float number)
{
const float x2 = number * 0.5F;
const float threehalfs = 1.5F;
// 使用联合体对浮点数的bit进行二进制整数的运算
union
{
float f;
unsigned long i;
} conv= { .f = number };
// 此处某魔法数字出现
conv.i= 0x5f3759df - ( conv.i >> 1 );
// 下面这行迭代多次可以获得精度
conv.f*= ( threehalfs - ( x2 * conv.f * conv.f ) );
// 返回浮点数的结果
return conv.f;
}
## 逐二进制位开方法
这种算法对每个二进制位进行位运算和加减法，通过一个bit一个bit去试验，来测出平方根的每个bit的值。
逐二进制开方算法假定对数值x求出的平方根结果的位数肯定小于x的位数的一半。因此对于32位整数求平方根，可以只测试其中的16个bit。
int int_sqrt(int x)
{
int v_bit = 15;
int n = 0; // 结果
int b = 0x8000; // 对应的位
if(x <= 1) return x;
while(b)
{
int temp = ((n << 1) + b) << v_bit--;
if(x >= temp)
{
n += b;
x -= temp;
}
b >>= 1;
}
return n;
}
这个算法不需要做整数除法计算，只需要做若干加法和移位，在合适的条件下，编译器可以通过展开有限的循环来实现优化。**运行效率比牛顿开方法高得多**。
# 整数开方的典型应用
整数开方可以被广泛用于嵌入式平台，典型例子：[步进电机加减速的计算](https://www.0xaa55.com/thread-26071-1-1.html)。
**整数开方算法可以应用于定点数开方**。基于NASM汇编语言的[编译时间光线追踪渲染](https://github.com/0xAA55/NASMCompileTimeRayTracing)计算中，对向量做单位化计算，使用了整数开方法。NASM编译器不支持编译期间的浮点数计算（虽然支持浮点数常量的生成），我使用定点数进行小数的计算。
寻路算法的临时内存池分配。寻路算法在有时候需要存储一整条直线上的每一个距离的数据，此时使用勾股定理可以根据线段两点的坐标求出线段的长度，而沟谷定理计算依赖开方算法。
tomzbj 发表于 2024-8-29 10:47
我也贴个我的整数开平方
可以用许多平台都有的内建clz指令, 找前导零的个数, 然后用32减去前导零个数, 再除 ...
又测试了一下，当x的二进制位数n是偶数时，似乎取y的初始值用n/2-1比n/2的效果更好，大部分情况迭代三次就能得到精确值，用n/2则一般需要四次。有空跑个完整测试看看。 我也贴个我的整数开平方
可以用许多平台都有的内建clz指令, 找前导零的个数, 然后用32减去前导零个数, 再除以2, 作为第一步的估计根.
后面只要再进行三次牛顿迭代即可.
实测整个32位范围内可能偶尔会差1, 不过一般场合这个精度也够用了.
肯定比用循环快一些, 去掉了循环判断和跳转的额外开销.
如果对精度要求不高还可以再减少一次迭代.
int int_sqrt(register int x)
{
register int n = __builtin_clz(x);
n = (32 - n) / 2;
register int y = 1 << n;
y = (y + x / y) / 2;
y = (y + x / y) / 2;
y = (y + x / y) / 2;
return y;
}
本帖最后由 系统消息 于 2020-11-25 11:27 编辑
A5大婶，你的int版的开平方支持负数不？结果截断取整，还是向下取整？
我之前也搜集过这方面的算法来研究，发现跟你的不一样，我也发出来共享一下吧：
// 无符号整数快速开方
unsigned sqrtu(unsigned x)
{
if (x == 0)return 0;
unsigned n = 0;
unsigned t = (x >> 30);
x <<= 2;
if (t > 1)
{
n++;
t -= n;
}
for (unsigned i = 15; i > 0; i--)
{
n <<= 1;
t <<= 2;
t += (x >> 30);
unsigned u = n;
u = (u << 1) + 1;
x <<= 2;
if (t >= u)
{
t -= u;
n++;
}
}
return n;
}
// 浮点数快速开方
float sqrtf(float x)
{
union {
int i;
float f;
};
float xhalf = 0.5f * x;
i = 0x5f375a86 - ((int&)x >> 1); // 计算第一个近似根
f *= 1.5f - xhalf * f * f; // 牛顿迭代法
f *= 1.5f - xhalf * f * f; // 再次牛顿迭代（可选，提高精度用）
return xhalf * f + 0.5f / f; // 两种方式求倒数刚好误差在方向相反，最后平均一下就精确了（跟C库相同精度)
}
系统消息 发表于 2020-11-25 11:25
A5大婶，你的int版的开平方支持负数不？结果截断取整，还是向下取整？
我之前也搜集过这方面的算法来研究， ...
负数开方可还行……
那个浮点数的看起来就是快速反开平方，并且迭代了两次。
整数那个看起来像我那个按bit开方的算法 0xAA55 发表于 2020-11-25 17:17
负数开方可还行……
那个浮点数的看起来就是快速反开平方，并且迭代了两次。
整数的差别还是挺明显的，浮点数版的那个十六进制常量不一样，我之前查资料的时候说这个常量精度高一点，迭代两次主要是提高精度（也可以选择少一次），我的求倒数的时候，用了两种方式平均来减小误差，达到了跟C库几乎相同的精度）。 系统消息 发表于 2020-11-25 18:29
整数的差别还是挺明显的，浮点数版的那个十六进制常量不一样，我之前查资料的时候说这个常量精度高一点， ...
你可以参考这里：https://en.wikipedia.org/wiki/Methods_of_computing_square_roots
float sqrt_approx(float z)
{
int val_int = *(int *)&z;
val_int -= 1 << 23;
val_int >>= 1;
val_int += 1 << 29;
return *(float *)&val_int;
} new_starter 发表于 2020-12-7 22:55
你可以参考这里：https://en.wikipedia.org/wiki/Methods_of_computing_square_roots
这个我记得只能用于2的整数次方吧 本帖最后由 tlwh163 于 2024-6-23 17:43 编辑
你们这 左移2 右移30 是不是就是神奇的 循环左移?
麻烦介绍下原理
牛顿法可能有点问题 做除法的时候 要截断小数部分 不能让它四舍五入
x1 = ( x0 + s / x0 ) >> 1;
---------------------------------
x1 = ( x0 + Fix(s / x0) ) >> 1;
也许C语言是截断的 我这里VFB不是... tomzbj 发表于 2024-8-30 09:47
又测试了一下，当x的二进制位数n是偶数时，似乎取y的初始值用n/2-1比n/2的效果更好，大部分情况迭代三次 ...
建议 / 2 改成 >> 1，/ 是影响性能的罪魁祸首。 tlwh163 发表于 2024-6-23 15:09
你们这 左移2 右移30 是不是就是神奇的 循环左移?
麻烦介绍下原理
BASIC语言的 / 是浮点数除法（对应 fdiv 指令），整数除法要用 \（对应 idiv 指令）。不过这类算法中最好是用 >> 1 来优化性能。 YY菌 发表于 2024-9-5 09:37
建议 / 2 改成 >> 1，/ 是影响性能的罪魁祸首。
对现在的编译器, 这种地方基本上不会有区别了... 你可以实测一下看看 tomzbj 发表于 2024-9-5 23:26
对现在的编译器, 这种地方基本上不会有区别了... 你可以实测一下看看
再智能的编译也是不可能没有区别的，/ 2 和 >> 1 的负数结果是不一样的，编译器不可能给你平替，它又不晓得你会不会用到负数。
页:
[1]