快速平方根和平方计算

在一些3D代码中，计算需要极快的速度，甚至不能忍受整型乘除运算，更别说浮点运算了。但是几乎每一个3D实现中都无法避免平方根和平方运算，因为这些代码必须计算距离。平方运算只是一个浮点乘法，不得已可以忍受，而平方根运算却是不可忍受的低效。快速平方根算法应运而生。下面是VC6中的一个实现：

__declspec( naked ) float fast_sqrt( float x ) { __asm { SUB ESP, 4 MOV EAX, [ESP+8] SUB EAX, 0x3F800000 SAR EAX, 1 ADD EAX, 0x3F800000 MOV [ESP], EAX FLD DWORD PTR [ESP] ADD ESP, 4         RET }    }

粗看上去，这段代码极简单，很难相信它能计算平方根，不过事实证明它能。同时，这段代码几乎是我看过的最晦涩的代码，我几乎花了两个小时才最终弄明白它的含义，编写这段代码的程序员毫无疑问是个高手中的高手。下面解释这段代码。

根据float格式（参见5.2.1节），x可以形式化地表示为：

其中f是尾数；n是指数；E是指数偏移；1是隐含位。

float的指数偏移是127（即0x7F），位于24~30位，即0x3F800000，因此代码中SUB EAX, 0x3F800000和ADD EAX, 0x3F800000的作用就是去除和添加E。那么实质计算只有一条指令SAR EAX, 1完成。此时x的形式是：

那么这个形式的平方根是：

利用泰勒级数在1.0处展开 去掉根号，有 

下面将指数分奇偶两种情形讨论SAR EAX,1指令的含义。

（1）指数是偶数时

设n = 2k，平方根的形式是： 

而根据float格式，SAR EAX, 1的效果是将指数和尾数分别除以2，因此有： 

两者结果一样。可见，在指数是偶数时，SAR EAX, 1实际上就是在计算平方根的泰勒级数的一次项。

（2）指数是奇数

设n = 2k+1，平方根的形式是： 

而SAR EAX, 1的效果是： 

之所以如此，是因为指数是奇数时，指数的最后一位左移后进入尾数，相当于在尾数部分添加2k-1。也就是说，此时除了泰勒级数误差之外，还有一个代替误差： 

相对误差约： 

最大约5.62%。最终的误差是这个误差与泰勒级数展开误差的混合。

从上述分析可以看出，相对误差只与尾数构成有关，与取值范围无关，因此下面以[0,1]区间做试验，结果见表13-1。

表13-1  [0,1]平方根值

x	fast_sqrt	sqrt	%
0.0000000	0.0000000	0.0000000
0.1000000	0.3250000	0.3162278	2.7740209
0.2000000	0.4500000	0.4472136	0.6230575
0.3000000	0.5500000	0.5477226	0.4158006
0.4000000	0.6500000	0.6324555	2.7740209
0.5000000	0.7500000	0.7071068	6.0660190
0.6000000	0.8000000	0.7745967	3.2795545
0.7000000	0.8500000	0.8366600	1.5944345
0.8000001	0.9000000	0.8944272	0.6230575
0.9000001	0.9500000	0.9486833	0.1387951

最大误差约6%左右，对于一般的3D计算而言，这个精度可以满足要求。

平方运算的原理与平方根运算类似，不再赘述。代码如下：

__declspec( naked ) float fast_x2( float x ) { __asm { SUB ESP, 4 MOV EAX, [ESP+8] SUB EAX, 0x3F800000 SAL EAX, 1 ADD EAX, 0x3F800000 MOV [ESP], EAX FLD DWORD PTR [ESP] ADD ESP, 4 RET }    }

试验结果见表13-2。最大误差有11%（表中没有出现）。

 

表13-2  [0,1]平方值

x	fast_x2	x2	%
0.1000000	0.0093750	0.0100000	-6.2500029
0.2000000	0.0375000	0.0400000	-6.2500029
0.3000000	0.0875000	0.0900000	-2.7777750
0.4000000	0.1500000	0.1600000	-6.2500029
0.5000000	0.2500000	0.2500000	0.0000000
0.6000000	0.3500000	0.3600000	-2.7777750
0.7000000	0.4500000	0.4900001	-8.1632685
0.8000001	0.6000001	0.6400001	-6.2499930
0.9000001	0.8000002	0.8100002	-1.2345664

另外，平方运算有一个限制，那就是0.0时无法用这种算法，需要剔除。注意，IEEE的0有两种：+0和-0，因此检测代码是：

if( ( *(unsigned int*)&x == 0 ) || ( *(unsigned int*)&x == 0x80000000 ) ) return 0.0f;

不要使用下列代码：

if( x == 0.0f ) return 0.0f;

这会导致编译器使用浮点指令。