解决方案
是,能够即使没有诡计:
1),用于速度牺牲精度:所述SQRT算法是迭代的,与更少的迭代重新实施
2)查找表:要么只是迭代的起点,或者插值相结合,让你一路都有
3)缓存:你总是sqrting有限组值相同?如果是这样,缓存可以很好地工作。我发现在正在计算大量形状相同大小相同的事情图形应用这是有用的,所以结果可以有效地进行高速缓存。
其他提示
这里有一个很大的对比表: http://assemblyrequired.crashworks.org/timing-square-root/
长话短说,SSE2的ssqrts大约为2x比FPU FSQRT更快,以及一个近似+迭代是约4倍比更快(8X整体)。
另外,如果你想采取单精度开方,确保实际上你在说什么。我听说过至少一个编译器,将float参数转换为双的,称之为双精度开方,再转换回浮动。
您很有可能会改变,以获得更多的速度提升您的算法的不是改变他们的实施的:尝试调用sqrt()
少,而不是拨打电话更快。 (如果你认为这是不可能的 - 你提到的sqrt()
的改进仅仅是:算法的用于计算平方根的改进。)
由于它的使用非常普遍,很可能是你的标准库的实现sqrt()
的是接近最优的一般情况。除非你有一个受限制的领域(例如,如果你需要更少的精度),其中该算法可以采取一些快捷键,这是不太可能有人想出了这是更快的实现。
请注意,由于该功能使用你的执行时间的10%,即使你设法拿出一个实现,只需要的std::sqrt()
的75%的时间,这将仍然只是把你的执行时间下降了<强> 2.5%即可。对于大多数应用程序的用户甚至不会注意到这一点,但如果他们使用的手表来衡量。
如何精确你需要你的sqrt
是?你可以得到合理的近似速度非常快:看到的Quake3的优秀逆灵感平方根功能(注意,代码是GPL的,因此可能不希望将其直接集成)。
不知道如果你固定这一点,但我读过有关它之前,它似乎是最快的事情要做更换 sqrt
功能用嵌入式组件版本;
你可以看到的描述负荷的替代品 在这里,.
最好的是这个片段中的魔法:
double inline __declspec (naked) __fastcall sqrt(double n)
{
_asm fld qword ptr [esp+4]
_asm fsqrt
_asm ret 8
}
这是关于4.7个比标准 sqrt
电话用相同的精确度。
下面是与一个查找仅8KB的表的快速方法。的错误是结果的〜0.5%。您可以轻松地扩大该表,从而减少错误。约5倍的运行速度比正常速度更快的sqrt()
// 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;
}