13.8  快速平方根和平方计算

在一些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的形式是：

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

（1）指数是偶数时

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

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

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

（2）指数是奇数

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

而SAR EAX, 1的效果是：