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

[阿呆_]

又换工作了， 现在用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, 4

end;

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;

[原来是果果呀]

顶一下

学习优化的思路顶一下

学习优化的思路

[xihii6]

够强，

[moonlight88162348]

大量计算，才能显现出作用的。

[wxzty]

顶一下

[zgbj999]

[Quote=引用 68 楼 cliffbaby 的回复:]
每天回帖即可获得10分可用分！
[/Quote]
每天回帖即可获得10分可用分！

[jiashie]

[Quote=引用 48 楼 gordon3000 的回复:]

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

同问。

[lwglucky]

意义不大。计算科学不是这么研究的。。建议看看数值分析的书籍和仿真方面的理论书籍。

[liumu1209]

高手啊

[liumu1209]

高手啊

[allen1981813]

有必要么？

[linguangfei2007]

路过,路过,学习

[minyi123yimin]

想學這語言呢
謝謝lz啦

[tomwang1978]

现在软件的要求是如何能够正常的实现功能，服务器性能强了n倍

[kakaroot_cl]

牛人 尽管看不懂 还是帮你顶个

[iamyg]

每天回复既有10分可用分

[Rexdy]

三角函数对我来说不是很常用，楼主加油

[Cjinlei]

很实用。。。

[ahejn]

用DELPHI做程序，还对速度要求这么苛刻，貌似不太合适。

[StevenWURGBLink]

^_^，^_^

[YFLK]

快慢是一回事,计算精度是另一回事,你的程序计算精度能满足要求吗?