3D *_*der 30 c algorithm math sqrt
我需要快速的整数平方根,不涉及任何明确的除法.目标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个指令的循环(bitset,mul,compare,store,shift).缓存中可能没有完全展开的空间(但这将是部分展开的主要候选者[例如,4循环而不是16],当然).因此,成本是16*5 = 80指令(加上循环开销,如果没有展开).如果完全交错,则仅花费80(对于最后指令为+2)周期.
我可以在82个周期内获得一些其他sqrt实现(仅使用add,mul,bitshift,store/cmp)吗？
常问问题:
为什么不依靠编译器来生成一个好的快速代码？
该平台没有可用的C-> RISC编译器.我将把当前的参考C代码移植到手写的RISC ASM中.
您是否对代码进行了分析,以确定它add是否真的成为瓶颈？
不,没有必要.目标RISC芯片大约为20 MHz,因此每条指令都很重要.核心环(计算射击器和接收器贴片之间的能量转移形状因子),mul使用它,将在每个渲染帧运行~1000次(当然,假设它足够快),每秒高达60,000,整个演示大约有1,000,000次.
您是否尝试过优化算法以删除sub？
是的,我已经这样做了.事实上,我已经摆脱了2 shift秒和许多分裂(移除或替换为移位).即使在我的千兆赫笔记本上,我也可以看到巨大的性能提升(与参考浮动版相比).
申请是什么？
它是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个操作,包括ifor循环中的测试和增量以及3个赋值,但是当在ASM中编码时,其中一些可能会消失？下面的代码中有6个操作,如果算作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用于加法+ 1用于移位= 18次循环)是不值得的,除非它给出超过4位.
所以,我的建议是通过一种建议的方法(8*4 = 32个周期)建立8位的答案,然后执行牛顿方法的一次迭代(18个周期),将位数加倍到16个.这是一个总数50个循环(加上可能额外的4个循环,在应用牛顿方法之前获得9个位,加上可能有2个循环来克服牛顿方法偶尔遇到的过冲).这是最多56个周期,据我所知,可以看到任何其他建议.
第二次更新
我编写了混合算法的想法.牛顿方法本身没有开销; 你只需申请并加倍有效数字.在应用牛顿方法之前,问题是要有可预测的有效位数.为此,我们需要弄清楚答案中最重要的部分出现在哪里.使用另一张海报给出的快速DeBruijn序列方法的修改,我可以在估计中在大约12个周期中执行该计算.另一方面,知道答案的msb的位置可以加速所有方法(平均而非最坏的情况),所以无论如何它似乎都值得.
在计算了答案的msb之后,我运行了上面建议的一些算法,然后用一两轮Newton方法完成它.我们通过以下计算决定何时运行牛顿方法:根据评论中的计算,答案的一位需要大约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 works\n");
}
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程序平均比isqrt我给出的上一个程序快40%.
如果乘法与加法和移位的速度相同（或更快！），或者如果您缺少寄存器中包含的快速移位指令，那么以下内容将没有帮助。否则：
您temp*temp在每个循环周期上重新计算，但是，这与它们的位不重叠temp = res | add相同，并且 (a) 您已经在前一个循环周期上进行了计算，并且 (b)具有特殊的结构（它总是只是一点）。因此，利用和 您已经拥有 和的事实， 和只是向左移动了一位 的两倍，并且比 中的单个位的位置多左移了 1 个位置，你可以摆脱乘法：res + addres*resadd(a+b)^2 = a^2 + 2ab + b^2a^2b^2b2abab
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 标准，所有算术都需要在执行算术运算之前将较窄的整数类型提升（转换）为涉及的最宽的类型，如果需要的话，随后进行后续的反向转换。（这当然可以通过足够智能的编译器进行优化，但为什么要冒这个风险呢？）