【算法】整数开平方算法之一——只使用逻辑运算和移位
# 定义**整数开平方**被广泛用于嵌入式的步进电机控制计算、游戏的寻路算法的内存池分配等多个领域。对比常见的**浮点数开平方**，整数开平方不需要依赖浮点数计算单元，且能在无硬件浮点数单元的嵌入式环境里实现高效率的开平方。
在[数论](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编译器不支持编译期间的浮点数计算（虽然支持浮点数常量的生成），我使用定点数进行小数的计算。
寻路算法的临时内存池分配。寻路算法在有时候需要存储一整条直线上的每一个距离的数据，此时使用勾股定理可以根据线段两点的坐标求出线段的长度，而沟谷定理计算依赖开方算法。
本帖最后由 系统消息 于 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的整数次方吧
页:
[1]