是的,即使没有欺骗也是可能的:
为了速度而牺牲准确性:sqrt算法是迭代的,可以重新实现更少的迭代。
查找表:可以用于迭代的起点,或与插值结合使用以帮助你一直到达目标。
缓存:你总是对同一组有限的值进行平方根运算吗?如果是,缓存可以很好地工作。我在图形应用程序中发现这很有用,因为相同大小的许多形状正在计算相同的东西,因此结果可以被有用地缓存。
来自未来的11年问候。
考虑到这仍然偶尔会得到投票,我想添加一条关于性能的注释,现在比以前更受内存访问的限制。当优化像这样的东西时,您绝对必须使用逼真的基准测试(理想情况下是整个应用程序)-您应用程序的内存访问模式将对查找表和缓存等解决方案产生巨大影响,仅仅比较您优化版本的“循环”将使您大失所望:还很难将程序时间分配给单个指令,您的分析工具可能会误导您。
在相关说明上,请考虑使用simd /向量化指令计算平方根,例如_mm512_sqrt_ps或类似产品,如果它们适合您的用例。
请参阅英特尔的optimisation reference manual第15.12.3节,其中描述了近似方法,使用向量化指令,这些方法可能非常适合其他体系结构。
1/sqrt(x)
),而不是sqrt(x)
。 - James1/(1/sqrt(x))
可以作为一种妥协。另外,当你提到 LUTs 时,我想知道你是否考虑了在 C++ 中实现 Gal 的精确表格?我很想了解更多信息。 - rytis这里有一个很好的比较表格: http://assemblyrequired.crashworks.org/timing-square-root/
简而言之,SSE2的ssqrts比FPU fsqrt快约两倍,而一种近似+迭代方法比它更快约4倍(总共8倍)。
此外,如果你试图获取单精度平方根,请确保你实际上正在得到它。我听说至少有一个编译器会将浮点参数转换为双精度,调用双精度sqrt,然后再转换回浮点数。
改变你的算法比改变它们的实现可能会带来更多的速度提升:尝试减少对sqrt()
的调用而不是让调用变得更快。 (如果你认为这不可能-你提到的sqrt()
的改进只是用于计算平方根的算法的改进。)
由于sqrt()
使用非常频繁,你标准库的实现在一般情况下可能已经接近最优了。除非你有一个受限的域(例如,如果你需要较少的精度),算法可以采取一些捷径,否则很难想出一个更快的实现。
请注意,由于该函数使用了10%的执行时间,即使你成功地想出了一个只需75%的时间就能完成std::sqrt()
的实现,这仍然只会将你的执行时间降低 2.5%。对于大多数应用程序用户来说,他们甚至不会注意到这一点,除非他们使用手表来测量。
rsqrtss
,ARMv7上的rsqrte
等)。 - Stephen Canon不知道你是否已经解决了这个问题,但我之前读过相关内容,似乎最快的方法是用内联汇编版本替换 sqrt
函数;
你可以在这里看到许多替代方案的介绍。
最好的方法是这个神奇的代码片段:
double inline __declspec (naked) __fastcall sqrt(double n)
{
_asm fld qword ptr [esp+4]
_asm fsqrt
_asm ret 8
}
这个方法比标准的sqrt
函数快大约4.7倍,精度相同。
esp+4
可能不是 n
,并且内联 asm 中有一个 ret
)。 它强制要求您的输入通过存储/重新加载来进行。这又增加了大约5个周期的延迟,加上Haswell上10-23个周期的fsqrt
延迟。sqrtpd
(double)为16个周期的延迟,而sqrtps
(single)为11个周期的延迟。 除非您的编译器真的很差(也许您没有进行优化编译?),否则您的编译器不应选择这种方式。此外,ret 8
?那是一个3uop指令。 - Peter Cordessqrt(int)
和std::sqrt(double)
,这是错误的。实际上,内部的std::sqrt算法速度对于所有参数类型并不相同。顺便说一句,我重新实现了他的测试,只有2或3个变量略微更快,但与std变体相比,精度大大降低了。 - metablaster// LUT for fast sqrt of floats. Table will be consist of 2 parts, half for sqrt(X) and half for sqrt(2X).
const int nBitsForSQRTprecision = 11; // Use only 11 most sagnificant bits from the 23 of float. We can use 15 bits instead. It will produce less error but take more place in a memory.
const int nUnusedBits = 23 - nBitsForSQRTprecision; // Amount of bits we will disregard
const int tableSize = (1 << (nBitsForSQRTprecision+1)); // 2^nBits*2 because we have 2 halves of the table.
static short sqrtTab[tableSize];
static unsigned char is_sqrttab_initialized = FALSE; // Once initialized will be true
// Table of precalculated sqrt() for future fast calculation. Approximates the exact with an error of about 0.5%
// Note: To access the bits of a float in C quickly we must misuse pointers.
// More info in: http://en.wikipedia.org/wiki/Single_precision
void build_fsqrt_table(void){
unsigned short i;
float f;
UINT32 *fi = (UINT32*)&f;
if (is_sqrttab_initialized)
return;
const int halfTableSize = (tableSize>>1);
for (i=0; i < halfTableSize; i++){
*fi = 0;
*fi = (i << nUnusedBits) | (127 << 23); // Build a float with the bit pattern i as mantissa, and an exponent of 0, stored as 127
// Take the square root then strip the first 'nBitsForSQRTprecision' bits of the mantissa into the table
f = sqrtf(f);
sqrtTab[i] = (short)((*fi & 0x7fffff) >> nUnusedBits);
// Repeat the process, this time with an exponent of 1, stored as 128
*fi = 0;
*fi = (i << nUnusedBits) | (128 << 23);
f = sqrtf(f);
sqrtTab[i+halfTableSize] = (short)((*fi & 0x7fffff) >> nUnusedBits);
}
is_sqrttab_initialized = TRUE;
}
// Calculation of a square root. Divide the exponent of float by 2 and sqrt() its mantissa using the precalculated table.
float fast_float_sqrt(float n){
if (n <= 0.f)
return 0.f; // On 0 or negative return 0.
UINT32 *num = (UINT32*)&n;
short e; // Exponent
e = (*num >> 23) - 127; // In 'float' the exponent is stored with 127 added.
*num &= 0x7fffff; // leave only the mantissa
// If the exponent is odd so we have to look it up in the second half of the lookup table, so we set the high bit.
const int halfTableSize = (tableSize>>1);
const int secondHalphTableIdBit = halfTableSize << nUnusedBits;
if (e & 0x01)
*num |= secondHalphTableIdBit;
e >>= 1; // Divide the exponent by two (note that in C the shift operators are sign preserving for signed operands
// Do the table lookup, based on the quaternary mantissa, then reconstruct the result back into a float
*num = ((sqrtTab[*num >> nUnusedBits]) << nUnusedBits) | ((e + 127) << 23);
return n;
}
sqrt
更糟糕,尤其是针对测试sqrt
吞吐量的微基准测试。x86架构具有非常快速的硬件sqrt
,甚至还有更快的rsqrt
(近似倒数平方根),你可以通过牛顿-拉弗逊步骤来进一步优化。 - Peter Cordes
rsqrtss
(近似倒数平方根)会更快。 - Stephen Canonrsqrtps
和牛顿迭代的不错版本,对于单精度float
可达到+/-2ulp。参见 https://dev59.com/H1wY5IYBdhLWcg3w9r7a 了解-mrecip
编译器优化,它应该控制自动使用rsqrt
(但实际上只有 1/x 的rcp
)。 - Peter Cordes