最快的整数平方根[最少量的指令]
我需要快速整数平方根，不涉及任何明确的分割。 目标RISC体系结构可以在一个周期内执行像add，mul，sub，shift这样的操作（很好 - 操作的结果是在第三个周期中写入的 - 确实存在交错），所以使用这些操作并且速度很快的Integer算法将非常赞赏。
这是我现在所拥有的，我认为二分法搜索应该更快，因为下面的循环每次执行16次（不管值是多少）。 我还没有广泛调试它（但很快），所以也许可以在那里提前退出：
unsigned short int int_sqrt32(unsigned int x)
{
unsigned short int res=0;
unsigned short int add= 0x8000;
int i;
for(i=0;i<16;i++)
{
unsigned short int temp=res | add;
unsigned int g2=temp*temp;
if (x>=g2)
{
res=temp;
}
add>>=1;
}
return res;
}
看起来上述[在目标RISC的上下文中]的当前性能成本是5条指令（位集，多路，比较，存储，移位）的循环。 可能没有缓存中的空间可以完全展开（但肯定会是部分展开的主要候选对象[例如，4而不是16的循环]）。 所以，成本是16 * 5 = 80条指令（加上循环开销，如果不展开）。 如果完全交错，则其成本仅为80（最后指令为+2）周期。
我可以在82个周期之后获得其他的sqrt实现吗（只使用add，mul，bitshift，store / cmp）吗？
常问问题：
你为什么不依靠编译器来生成一个好的快速代码？
该平台没有工作的C-> RISC编译器。 我将把当前的参考C代码移植到手写的RISC ASM中。
你是否剖析了代码以确定sqrt实际上是否是一个瓶颈？
不，这是没有必要的。 目标RISC芯片大约为20MHz，因此每一个指令都是很重要的。 使用sqrt的核心循环（计算射手和接收者补丁之间的能量传输形式因子）将在每个渲染帧运行约1,000次（当然，假设它足够快），每秒高达60,000次，整个演示大约为1,000,000次。
您是否尝试过优化算法以删除sqrt ？
是的，我已经做到了。 事实上，我已经摆脱了2 sqrt s和很多部门（删除或更换）。 即使在我的千兆赫笔记本上，我也能看到巨大的性能提升（与参考浮点版本相比）。
什么是应用程序？
它是用于compo演示的实时渐进式细化渲染器。 这个想法是每帧有一个拍摄周期，所以它会明显地收敛，并且在每个渲染帧中看起来更好（例如，每秒高达60次，尽管SW光栅化器不会那么快[但至少它可以运行在与RISC并行的另一个芯片上 - 因此，如果渲染场景需要2-3帧，RISC将并行处理2-3帧辐射度数据。
你为什么不直接在目标ASM中工作？
由于光能传递是一个稍微涉及的算法，我需要Visual Studio的即时编辑和继续调试功能。 我在VS周末做了些什么（将几百次代码更改转换为仅用于整数的浮点运算）将在目标平台上花费6个月的时间，只需打印“调试”。
你为什么不能使用分区？
因为目标RISC比以下任何一个速度慢16倍：mul，add，sub，shift，compare，load / store（仅需要1个周期）。 所以，只有在绝对需要的情况下才会使用它（已经几次了，不幸的是，当不能使用换挡时）。
你可以使用查找表吗？
该引擎已经需要其他LUT，并且从主RAM复制到RISC的小缓存非常昂贵（并且绝对不是每个帧）。 但是，如果它给了我至少100到200％的sqrt提升，我可能会备用128-256字节。
sqrt的值范围是多少？
我设法将它降低到仅仅是32位int（4,294,967,295）
看看这里。
例如，在3（a）有这种方法，它适用于64-32位平方根，也可以简单转换为汇编：
/* by Jim Ulery */
static unsigned julery_isqrt(unsigned long val) {
unsigned long temp, g=0, b = 0x8000, bshft = 15;
do {
if (val >= (temp = (((g << 1) + b)<<bshft--))) {
g += b;
val -= temp;
}
} while (b >>= 1);
return g;
}
没有分割，没有乘法，只有位移。 但是，所花费的时间将会有些不可预测，特别是如果您使用分支（在ARM RISC条件指令上工作）。
一般来说，这个页面列出了计算平方根的方法。 如果你碰巧想产生一个快速倒数平方根（即x**(-0.5) ），或者只是对令人惊叹的优化代码感兴趣，请看看这个，这个和这个。
这与你的一样，但操作较少。 （我在你的代码循环中计算了9个操作数，包括在for循环和3个赋值中的测试和增量i ，但也许其中一些在ASM中编码时消失了？下面的代码中有6个ops，如果你计算g*g>n作为两个（不分配））。
int isqrt(int n) {
int g = 0x8000;
int c = 0x8000;
for (;;) {
if (g*g > n) {
g ^= c;
}
c >>= 1;
if (c == 0) {
return g;
}
g |= c;
}
}
我在这里。 如果展开循环并根据输入中最高的非零位跳转到适当的点，则可以消除比较。
更新
我一直在考虑使用牛顿的方法。 理论上，每次迭代的准确位数应该加倍。 这意味着当答案中的正确位数很少时，牛顿的方法比任何其他建议都差得多; 然而，如果答案中有很多正确的位，情况就会改变。 考虑到大多数建议似乎每个比特需要4个周期，这意味着牛顿方法的一次迭代（除法16个周期+加法+1个移位= 18个周期）是不值得的，除非它给出超过4位。
所以，我的建议是通过建议的方法之一（8 * 4 = 32个周期）建立8位答案，然后执行牛顿方法的一次迭代（18个周期），将位数加倍到16。 50个周期（加上可能需要额外的4个周期才能得到9位，然后再应用牛顿法），再加上2个周期来克服偶尔牛顿方法所经历的过冲。 这是最多56个周期，据我可以看到对手的任何其他建议。
第二次更新
我编码混合算法的想法。 牛顿的方法本身没有开销; 你只需申请并加倍有效位数。 问题是在应用牛顿方法之前有一个可预测的有效数字位数。 为此，我们需要确定答案的最重要部分将出现在哪里。 使用另一张海报给出的快速DeBruijn序列方法的修改，我们可以在我估计的约12个周期内执行该计算。 另一方面，知道答案的MSB的位置加快了所有方法（平均值，而不是最坏情况），所以无论如何它似乎都是值得的。
在计算答案的MSB后，我运行了上面提到的算法的几轮，然后用一两轮牛顿法完成它。 我们通过以下计算来决定何时运行牛顿法：根据评论中的计算，答案的一个位需要大约8个周期; 牛顿方法的一轮需要大约18个周期（除法，加法和移位，也可能是分配），所以如果我们要从中得到至少三个比特的话，我们应该只运行牛顿法。 因此，对于6位答案，我们可以运行线性方法3次得到3个有效位，然后运行牛顿方法1次得到另一个3.对于15位答案，我们运行线性方法4次得到4位，然后牛顿方法两次以获得另一个4然后另一个7.依此类推。
这些计算依赖于确切知道线性方法需要多少个周期才能获得一点，而牛顿方法需要多少个周期。 如果“经济学”发生变化，例如通过发现以线性方式建立位的更快方法，那么何时调用牛顿方法的决定将会改变。
我展开了这些循环，并将这些决定作为开关实现，我希望这些开关可以转换为汇编中的快速表查找。 我并不十分确定每种情况下我都有最少的周期数，所以可能进一步调整。 例如，对于s = 10，您可以尝试获得5位，然后应用一次牛顿方法而不是两次。
我已经彻底测试了算法的正确性。 在某些情况下，如果您愿意接受稍微不正确的答案，则可以进行一些额外的小加速。 在应用牛顿方法纠正与m^2-1形式的数字发生的偏离错误后，至少使用两个周期。 在开始时使用循环测试输入0，因为算法无法处理该输入。 如果你知道你永远不会取零的平方根，你可以消除这个测试。 最后，如果您在答案中只需要8位有效位，您可以放弃牛顿方法计算之一。
#include <inttypes.h>
#include <stdint.h>
#include <stdbool.h>
#include <stdio.h>
uint32_t isqrt1(uint32_t n);
int main() {
uint32_t n;
bool it_works = true;
for (n = 0; n < UINT32_MAX; ++n) {
uint32_t sr = isqrt1(n);
if ( sr*sr > n || ( sr < 65535 && (sr+1)*(sr+1) <= n )) {
it_works = false;
printf("isqrt(%" PRIu32 ") = %" PRIu32 "n", n, sr);
}
}
if (it_works) {
printf("it worksn");
}
return 0;
}
/* table modified to return shift s to move 1 to msb of square root of x */
/*
static const uint8_t debruijn32[32] = {
0, 31, 9, 30, 3, 8, 13, 29, 2, 5, 7, 21, 12, 24, 28, 19,
1, 10, 4, 14, 6, 22, 25, 20, 11, 15, 23, 26, 16, 27, 17, 18
};
*/
static const uint8_t debruijn32[32] = {
15, 0, 11, 0, 14, 11, 9, 1, 14, 13, 12, 5, 9, 3, 1, 6,
15, 10, 13, 8, 12, 4, 3, 5, 10, 8, 4, 2, 7, 2, 7, 6
};
/* based on CLZ emulation for non-zero arguments, from
* http://stackoverflow.com/questions/23856596/counting-leading-zeros-in-a-32-bit-unsigned-integer-with-best-algorithm-in-c-pro
*/
uint8_t shift_for_msb_of_sqrt(uint32_t x) {
x |= x >> 1;
x |= x >> 2;
x |= x >> 4;
x |= x >> 8;
x |= x >> 16;
x++;
return debruijn32 [x * 0x076be629 >> 27];
}
uint32_t isqrt1(uint32_t n) {
if (n==0) return 0;
uint32_t s = shift_for_msb_of_sqrt(n);
uint32_t c = 1 << s;
uint32_t g = c;
switch (s) {
case 9:
case 5:
if (g*g > n) {
g ^= c;
}
c >>= 1;
g |= c;
case 15:
case 14:
case 13:
case 8:
case 7:
case 4:
if (g*g > n) {
g ^= c;
}
c >>= 1;
g |= c;
case 12:
case 11:
case 10:
case 6:
case 3:
if (g*g > n) {
g ^= c;
}
c >>= 1;
g |= c;
case 2:
if (g*g > n) {
g ^= c;
}
c >>= 1;
g |= c;
case 1:
if (g*g > n) {
g ^= c;
}
c >>= 1;
g |= c;
case 0:
if (g*g > n) {
g ^= c;
}
}
/* now apply one or two rounds of Newton's method */
switch (s) {
case 15:
case 14:
case 13:
case 12:
case 11:
case 10:
g = (g + n/g) >> 1;
case 9:
case 8:
case 7:
case 6:
g = (g + n/g) >> 1;
}
/* correct potential error at m^2-1 for Newton's method */
return (g==65536 || g*g>n) ? g-1 : g;
}
在我的机器上进行轻量级测试（与我的机器相比毫无isqrt1 ），新的isqrt1例程平均比我之前的isqrt例程快40％。
如果乘法运算速度相同（或快于！）加法和移位，或者如果您没有快速移位的注册指令，那么以下几点将不会有帮助。 除此以外：
您将在每个循环周期重新计算temp*temp ，但temp = res | add temp = res | add ，它与res + add相同，因为它们的位不重叠，并且（a）您已经在前一个循环中计算了res*res ，并且（b） add具有特殊结构（它总是只有一个位）。 因此，使用(a+b)^2 = a^2 + 2ab + b^2的事实，并且您已经有a^2 ，并且b^2仅仅是一位向左偏移两倍单比特b ，和2ab只是a左移1个比单个位的位置的详细位置b ，你可以摆脱乘法：
unsigned short int int_sqrt32(unsigned int x)
{
unsigned short int res = 0;
unsigned int res2 = 0;
unsigned short int add = 0x8000;
unsigned int add2 = 0x80000000;
int i;
for(i = 0; i < 16; i++)
{
unsigned int g2 = res2 + (res << i) + add2;
if (x >= g2)
{
res |= add;
res2 = g2;
}
add >>= 1;
add2 >>= 2;
}
return res;
}
另外我猜想对所有变量使用相同类型（ unsigned int ）是一个更好的主意，因为根据C标准，所有算术都需要将较窄整数类型提升（转换）为算术运算之前涉及的最宽类型然后在必要时进行后续转换。 （这当然可以通过充分智能的编译器进行优化，但为什么要冒这个风险呢？）
上一篇: Fastest Integer Square Root [least amount of instructions]