您现在的位置:首页 >> 基础算法 >> window基础 >> 内容

Delphi中编写比FDIV指令更快的浮点除法函数

时间:2011/9/3 15:28:29 点击:

  核心提示:最近无聊中一不留神又在网上看到了Quake中的那段著名的Inverse Square Root代码,忽然产生了一个冲动,将那段代码改成了一个计算浮点除法的函数,不幸就得到了一段比两个浮点数直接相除更快...

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

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

先看一下传说中的那个Inverse sqrt函数的相当神奇的源码C/C++:
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;

作者:Idle_ 来源:转载
共有评论 0相关评论
发表我的评论
  • 大名:
  • 内容:
  • 盒子文章(www.2ccc.com) © 2022 版权所有 All Rights Reserved.
  • 沪ICP备05001939号