呃，不小心写了个比FDIV指令更快的浮点除法函数

阿呆_ 2010-05-31 11:43:57

加精

最近无聊中一不留神又在网上看到了Quake中的那段著名的Inverse Square Root代码，忽然产生了一个冲动，将那段代码改成了一个计算浮点除法的函数，不幸就得到了一段比两个浮点数直接相除更快的代码。先说说我是怎么改造的，这段Quake代码计算的是某个浮点数的平方根的倒数，如果我们用这段代码对除数进行处理，得到的结果平方一下，再乘以被除数，那么结果当然就是两个数直接相除的值了。看上去很白痴的行为，可是结果却出人意料，绕这么一大圈来得到结果，比直接两个数相除还快。
简单测试了一下，又发现了另一件有趣的事情，在Intel Core Duo cpu上我的函数大约快50%，可是在Intel P4上我的函数竟然快了220%（很狼狈的结果，在两种cpu上我的函数执行时间几乎相同），看来传说中P4有时比P3还慢的谣言并不是瞎说。

下面详细解说一下我的改造：

先看一下传说中的那个Inverse sqrt函数的相当神奇的源码：

float InvSqrt(float x){

   float xhalf = 0.5f * x;

   int i = *(int*)&x; // store floating-point bits in integer

   i = 0x5f3759d5 - (i >> 1); // initial guess for Newton's method

   x = *(float*)&i; // convert new bits into float

   x = x*(1.5f - xhalf*x*x); // One round of Newton's method

   return x;

}

改造第一步，当然是将其翻译成delphi：

function InvSqrt(x: Single): Single;

var

  xhalf: Single;

  i: Integer;

begin

  xhalf := x * 0.5;

  i := PInteger(@xhalf)^;

  i := $5f3759d5 - (i shr 1);

  x := PSingle(@i)^;

  Result := x * (1.5 - xhalf*x*x);

end;

由于本人对执行效率的变态追求(并不是一个好习惯)，我又进一步将其改成了汇编：

var

  OneHalf:       DWORD  = $3FC00000;



function f_invsqrt(v: Single): Single;

asm

         MOV        EAX, DWORD PTR [EBP+8]  // i := PInteger(@x)^;

         SUB        EAX, $800000            

         PUSH       EAX                     // xhalf := x*2

         ADD        EAX, $800000            

         MOV        EDX, $5F3759D5          

         SHR        EAX, 1

         SUB        EDX, EAX                // i := 5f3759d5-(i shr 1)

         PUSH       EDX

         FLD        DWORD PTR [ESP]         // x := PSingle(@i)^

         FLD        ST(0)

         FMUL       ST(0), ST(0)

         FMUL       DWORD PTR [ESP+4]

         FSUBR      DWORD PTR OneHalf

         FMULP                              // result := x * (1.5 - xhalf*x*x);

         ADD        ESP, 8

end;

再将其改成除法，用delphi写的话就是(当然，这段代码绝对比直接除慢得多的多的多):

function f_div(x, y: Single): Single;

begin

  if y<0 then 

  begin

    y := -y;

    x := -x;

  end;

  result :=f_invsqrt(y) * f_invsqrt(y) * x;

end;

最后优化后的代码如下：

function f_div(v1, v2: Single): Single;

asm

         /* 调整v1, v2的值: v2取绝对值并相应调整v1的正负 */

         MOV        EAX, DWORD PTR [EBP+8]  

         MOV        EDX, $80000000

         AND        EDX, EAX                

         AND        EAX, $7FFFFFFF          

         XOR        DWORD PTR [EBP+12], EDX 

         FLD        DWORD PTR [EBP+12]

         /* 计算v2的Inverse sqrt再平方再乘以v1 */

         SUB        EAX, $800000

         PUSH       EAX

         ADD        EAX, $800000

         MOV        EDX, $5F3759D5

         SHR        EAX, 1

         SUB        EDX, EAX

         PUSH       EDX

         FLD        DWORD PTR [ESP]

         FMUL       ST(0), ST(0)

         FMUL       ST(1), ST(0)

         FMUL       DWORD PTR [ESP+4]

         FSUBR      DWORD PTR OneHalf

         FMUL       ST(0), ST(0)

         FMULP

         ADD        ESP, 8

end;