有用吗？比Delphi原始函数快近3倍的的sin()和cos()

Question

又换工作了， 现在用c/c++。 以后更少有时间在delphi露面了。 看着一堆过去写的delphi代码， 浪费啊，似乎没什么机会用到了。 干脆挑选一下， 以后会不定期发上来给大家分享， 混个脸熟吧。下面实现的sin和cos的原理相当简单， 事先计算好各个角度对应的sin值存在内存表中， 调用函数即是查表， 直接取出对应的值返回即可。 这样用简单的整数地址运算（只牵涉到加减和位运算--大部分都是一个时钟周期一条指令而且可以并行）代替了缓慢的浮点计算，速度当然快多了（题外话，现代CPU的浮点计算并不慢， 因为并行性能好， 在某些情况下甚至比相应的整数运算（整数模拟的定点数）更快， 不过这只是在“某些情况下”， 再加上浮点和整数之间的转换极其耗时， 所以总的来说浮点数运算还是比不上整数的）。为了更进一步提高速度， 我对sin和cos的参数定义做了一点小小的修改， sin和cos接受的参数不是浮点的弧度而是整数的角度（应该更符合普通人的习惯吧），另一点不那么符合大众习惯的是一个圆周不再是360度而是分成了512度（嘿嘿）。 由于现代CPU的特点， 访问的内存数据最好能够常驻L1或L2 cache（访问L1 cache耗时1-3时钟周期， L2 cache在20周期左右， 而访问主内存一般都在200时钟周期以上，这也是过去很多通过查表来替代某些耗时指令的算法今天来运行的话还没有直接使用那个耗时指令来得快的原因），这就要求我们的内存表尽可能小。下面程序中用到的内存表占用516字节，可以很轻松地加载进L1 cache后还剩大把空间供其它程序使用，这使得我们的内存表更容易驻留在L1 cache不被交换出去。表中保存的是1/4圆周的sin值，其它角度都可以通过换算映射到这张表里， 表中每一项是个整数，记录了包括小数点前1位和小数点后30位（二进制）精度的sin值。 代码是用汇编写的，并进行了我力所能及的最大优化：整个函数没有分支跳转指令（如果用delphi写的话代码中必然存在不止一个的if..else结构，要知道现代CPU对分支预测失败的惩罚是相当严厉的，特别是在CPU流水线越来越长的今天，一次分支预测失败耗时在40-200时钟周期，都够让FSINCOS指令完成一大半了），映射角度到1/4圆周的代码，结果正负调整的代码都是通过一些位操作实现的，这就是将圆周定为512度的原因和好处了。具体代码如下：
  SinLUT: array [0..128] of Cardinal = (  // Degrees for full circle: 0 -- 512 clockwise
        $00000000, $00C90E90, $0192155F, $025B0CAF, $0323ECBE, $03ECADCF, $04B54825, $057DB403, 
        $0645E9AF, $070DE172, $07D59396, $089CF867, $09640837, $0A2ABB59, $0AF10A22, $0BB6ECEF, 
        $0C7C5C1E, $0D415013, $0E05C135, $0EC9A7F3, $0F8CFCBE, $104FB80E, $1111D263, $11D3443F, 
        $1294062F, $135410C3, $14135C94, $14D1E242, $158F9A76, $164C7DDD, $17088531, $17C3A931, 
        $187DE2A7, $19372A64, $19EF7944, $1AA6C82B, $1B5D100A, $1C1249D8, $1CC66E99, $1D79775C, 
        $1E2B5D38, $1EDC1953, $1F8BA4DC, $2039F90F, $20E70F32, $2192E09B, $223D66A8, $22E69AC8, 
        $238E7673, $2434F332, $24DA0A9A, $257DB64C, $261FEFFA, $26C0B162, $275FF452, $27FDB2A7, 
        $2899E64A, $29348937, $29CD9578, $2A650525, $2AFAD269, $2B8EF77D, $2C216EAA, $2CB2324C, 
        $2D413CCD, $2DCE88AA, $2E5A1070, $2EE3CEBE, $2F6BBE45, $2FF1D9C7, $30761C18, $30F8801F, 
        $317900D6, $31F79948, $32744493, $32EEFDEA, $3367C090, $33DE87DE, $34534F41, $34C61236, 
        $3536CC52, $35A5793C, $361214B0, $367C9A7E, $36E5068A, $374B54CE, $37AF8159, $3811884D, 
        $387165E3, $38CF1669, $392A9642, $3983E1E8, $39DAF5E8, $3A2FCEE8, $3A8269A3, $3AD2C2E8, 
        $3B20D79E, $3B6CA4C4, $3BB6276E, $3BFD5CC4, $3C42420A, $3C84D496, $3CC511D9, $3D02F757, 
        $3D3E82AE, $3D77B192, $3DAE81CF, $3DE2F148, $3E14FDF7, $3E44A5EF, $3E71E759, $3E9CC076, 
        $3EC52FA0, $3EEB3347, $3F0EC9F5, $3F2FF24A, $3F4EAAFE, $3F6AF2E3, $3F84C8E2, $3F9C2BFB, 
        $3FB11B48, $3FC395F9, $3FD39B5A, $3FE12ACB, $3FEC43C7, $3FF4E5E0, $3FFB10C1, $3FFEC42D, 
        $40000000);function sin512(angle: Integer): Single; 
asm
         // if angle in upper half circle then
         //   EDX = -1
         //   map angle to 0..255 (lower half circle)
         MOV       EDX, 256
         AND       EDX, EAX
         CMP       EDX, 1
         SBB       EDX, EDX
         XOR       EDX, $FFFFFFFF
         AND       EAX, 255
         // if angle > 128 (right quarter circle) then
         //    map angle to 0..128 (left quarter circle)
         MOV       ECX, 128
         CMP       ECX, EAX
         SBB       ECX, ECX
         XOR       EAX, ECX
         SUB       EAX, ECX
         AND       ECX, 256
         ADD       EAX, ECX
         // lookup table
         MOV       EAX, DWORD PTR [EAX*4+SinLUT]
         // adjust +/- of result
         XOR       EAX, EDX
         SUB       EAX, EDX
         // change to float point result
         PUSH      EAX
         FILD      DWORD PTR [ESP]
         FSTP      DWORD PTR [ESP]
         // adjust decimal, = value shr 30
         SUB       DWORD PTR [ESP], $0F000000
         FLD       DWORD PTR [ESP]
         ADD       ESP, 4end;function cos512(angle: Integer): Single; // cos(x) = size(x+quarter circle)
asm
         ADD       EAX, 128
         JMP       sin512
end;另附计算SinLUT的代码：
type
  TDoubleU = record
    case Integer of
      0: (f: Double);
      1: (ul, uh: Cardinal);
  end;var
  i: Integer;
  y: TDoubleU;
  r: double;
begin
  for i := 0 to 128 do
  begin
    r := i / 256 * PI;
    y.f := sin(r);
    y.f := y.f + ldexp(1, 22);
    SinLUT[i] := y.ul;
  end;
end;

Answer 1

简单测试了一下,结果和楼主说的一样,比D版的快3培左右,不过要进行1000000次计算以后才明显的看出来,所以个人认为绝大多数情况下,D版的已经够用了.不过优化思路值得学习,下面是测试代码:
procedure TForm2.Button1Click(Sender: TObject);
var i,c:integer;
    t:Cardinal;
    r1:Single;
    r2:Extended;
    angle1:Integer;
    angle2:Extended;
begin
  angle1:=30;
  angle2:=30;
  c:=1000000;  Memo1.Clear;  t:=GetTickCount;
  for i := 0 to c do
    r1:=sin512(angle1);
  self.Memo1.Lines.Add(Format('sin512 %d',[GetTickCount-t])) ;  t:=GetTickCount;
  for i := 0 to c do
    r1:=cos512(angle1);
  self.Memo1.Lines.Add(Format('cos512 %d',[GetTickCount-t])) ;  t:=GetTickCount;
  for i := 0 to c do
    r1:=sin(angle1);
  self.Memo1.Lines.Add(Format('sin %d',[GetTickCount-t])) ;  t:=GetTickCount;
  for i := 0 to c do
    r1:=cos(angle1);
  self.Memo1.Lines.Add(Format('cos %d',[GetTickCount-t])) ;
end;测试结果:
1-------------
sin512 15
cos512 16
sin 47
cos 622-------------
sin512 16
cos512 15
sin 47
cos 473--------------
sin512 0
cos512 16
sin 47
cos 62测试环境:
CPU:奔腾 t2130 1.86GHz
内存:1G
windows xp sp3
Delphi2010

Answer 2

本帖最后由 starluck 于 2010-06-26 09:15:02 编辑

  

Answer 3

各位老大，计算结果不对
sin(x)<>sin512(x)

Answer 4

不知道楼主做啥软件的，要重写sin ，cos

Answer 5

咳咳～啰嗦两句疑似题外话:)LZ没有明确你所说的“浮点数”是单精型还是双精型。float才是真正的以浮点方式存储的浮点数，运算速度快，占存储空间也小。
而double的存储方式本质上是64位整型，运算速度比float慢，存储空间多一倍。所以如果不是对精度有特殊要求，建议应多使用float。

Answer 6

sin和cos接受的参数不是浮点的弧度而是整数的角度
==========================================
这个光整数不行吧？要是几点几度咋整呢？
就是用度分秒，秒也是应该浮点吧？

Answer 7

如果能加快效率和money，那就很不错

Answer 8

没用。主要问题是输入精度太差，一个算得上实用的三角函数算法，至少应该能以0.1度为单位计算，而您的sin512的输入精度只是360/512=0.703125度，比如我要计算sin(221.6)（单位均为角度），你实际计算的是sin(221)，误差太大了。另外，您这圆周硬性规定为512度的坐标太怪异，用户不太可能使用这种坐标，那他使用您的函数之前，肯定还要进行转换，把角度转换为您定义的坐标系中的“度”，这需要一次乘法和一次除法，也大大降低了您的算法的速度优势。软件实现三角函数算法，主要是cordic算法和多项式展开算法（诸如高阶泰勒级数展开，carmack在DOOM中的sin16就使用这种方法），查表法虽然快，但是输入精度上的限制很难解决，可以多次迭代来间接计算不能直接查表的值（比如利用sin(A+B)=sinAcosB+cosAsinB来分解输入值再分步计算），但是性能大幅度降低。

Answer 9

想增加输入精度其实很简单，只要相应扩大内存表项并继续细分圆周角度即可，比如定义一个圆周包括4096度，内存表也只需要相应扩大到4k而已，也可以一次性加载进L1 cache（L1 cache有16k)，已经够普通应用的精度了吧，针对密集性调用sin/cos的应用(比如实时图像旋转处理）时CPU cache命中率仍然会高于平均，对性能的增长还是很有帮助的。 
同样对于角度之间的映射也可以通过查表实现，没必要增加乘法和除法。

Answer 10

很强大，不过看不大懂，用sin512 算出的值都是-0.995184719562531，不知道为什么？

Answer 11

carmack在DOOM中的查表法finesine实际使用了10240表项。而且不是说你的L1 data有16K，表有4K，就能保证表在L1 data中，只有单任务系统才可能这么假设（也仅仅是可能）。即便只考虑当前任务，由于一个实际应用软件不可能象简单测试那样疯狂算sin/cos，而不管算出来干什么，它肯定进行一定量（不可能很多）的sin/cos计算后，把结果用于其他的计算过程，比如更新弹道轨迹，复制旋转后的像素等等，而新的计算往往会用到新的数据，结果就是可能把原来cache中的数据冲掉。我估计这就是为什么carmack只在DOS版的DOOM中使用查表法，在Windows版本的DOOM3中用更慢的泰勒级数展开算sin的主要原因（只是我的猜测，精度可能也是一个原因）。至于“同样对于角度之间的映射也可以通过查表实现，没必要增加乘法和除法。”，恐怕也这么简单。
精度高于整数的输入，基本上只能是一个浮点值（否则只能是你强迫用户为了用一个sin/cos函数而改变他的整个坐标系定义，那就本末倒置了），为了查表，需要先乘上一个比例（10或者100之类，取决于要求的精度）再转换成整数，还得做一次mod保证不超出表的范围。