二进制32位整数快速平方根

无心人

你39#代码存在问题 首行应该是or ecx, ecx吧

liangbch

是，不小心写错了，已修改 下面贴出在无心人的基础上，经我改写的3个FPU的版本。

3. double b32[] = {0.0,4294967296.0};
4. //实验证明zero5不能大于0.49999999999636
5. double zero5= 0.49999999999636;;
7. __declspec(naked)
8. DWORD __fastcall iSqrt_FPU1(DWORD n)
9. {
10. __asm
11. {
12. push ecx
14. shr ecx, 31
15. fld qword ptr [b32 + ecx * 8]
16. fild dword ptr [esp]
17. faddp st(1),st
18. fsqrt
19. fsub qword ptr [zero5]
20. fistp dword ptr [esp]
21. pop eax
22. ret
23. }
24. }
26. __declspec(naked)
27. DWORD __fastcall iSqrt_FPU2(DWORD n)
28. {
29. __asm
30. {
31. or ecx,ecx
32. jnz next00
33. mov eax,0
34. ret
35. next00:
36. push ecx
38. shr ecx, 31
39. fld qword ptr [b32 + ecx * 8]
40. fild dword ptr [esp]
41. faddp st(1),st
42. fsqrt
43. fstp qword ptr [esp-8]
45. fwait
46. mov ecx,dword ptr [esp-4]
47. mov edx,0xfff00000
48. mov eax,0xfffff
50. and edx,ecx //阶码
51. and eax,ecx //尾数
52. shr edx,20 //得到阶段码
53. add eax,0x100000 //设置位数最高有效位
54. mov ecx,1043
55. sub ecx,edx
56. shr eax,cl //得到整数
58. pop ecx
59. ret
60. }
61. }
63. __declspec(naked)
64. DWORD __fastcall iSqrt_FPU3(DWORD n)
65. {
66. __asm
67. {
68. push ecx
69. shr ecx, 31
70. fld qword ptr [b32 + ecx * 8]
71. fild dword ptr [esp]
72. faddp st(1),st
73. fsqrt
75. //设置FPU的取整方式 为了直接使用fistp浮点指令
76. FNSTCW [esp-4] // 保存协处理器控制字,用来恢复
77. FNSTCW [esp-6] // 保存协处理器控制字,用来修改
78. FWAIT
79. OR word ptr [esp-6], 0x0F00 // 改为 RC=11 使FPU向零取整
80. FLDCW [esp-6] // 载入协处理器控制字,RC场已经修改
82. fistp dword ptr [esp]
84. //恢复FPU的取整方式
85. FWAIT
86. FLDCW [esp-4]
88. pop eax
89. ret
90. }
91. }

复制代码

[ 本帖最后由 liangbch 于 2009-2-13 10:28 编辑 ]

gxqcn

我也写了个FPU的版本，有兴趣的朋友可把相关函数添加进来一并测试一下：

1. #include < stdio.h >
2. #include < stdlib.h >
3. #include < time.h >
4. #include < float.h >
6. typedef unsigned int UINT32;
7. typedef UINT32 __fastcall func_t( UINT32 );
9. _declspec(naked)
10. UINT32 __fastcall Sqrt_43F( UINT32 n )
11. {
12. __asm
13. {
14. xor eax, eax;
15. push eax;
16. push ecx;
18. fild qword ptr[esp];
19. fsqrt;
20. pop eax;
22. fistp dword ptr[esp];
23. pop eax;
25. ret;
26. }
27. }
29. #define FN(ID) Sqrt_##ID##F /* ex.: Sqrt_43F */
31. void test_fun( void )
32. {
33. UINT32 i, j;
34. UINT32 fID[] = { 43, /* ... */ };
35. func_t * funcTable[] = { FN(43), /* FN(..),... */ };
36. func_t * pfn;
38. #ifdef _DEBUG
39. UINT32 m = 256;
40. #else
41. UINT32 m = 1UL<<28;
42. time_t t;
43. #endif
45. printf( "\n\nn: 0-0x%x\n", m );
46. for ( i=0; i<sizeof(fID)/sizeof(fID[0]); ++i )
47. {
48. pfn = funcTable[i];
49. #ifdef _DEBUG
50. // check result
51. for ( j=0; j<=m; ++j )
52. {
53. printf( "Sqrt_%uF(%u)=%u\n", fID[i], j, pfn(j) );
54. }
55. system( "pause" );
56. #else
57. t = time(NULL);
58. for ( j=m; 0!=j; --j )
59. {
60. pfn(j);
61. }
62. t = time(NULL)-t;
64. printf( "%u#: %u s\n", fID[i], t );
65. #endif
66. }
67. }
69. int main(int argc, char* argv[])
70. {
71. _control87( RC_CHOP | PC_64, MCW_RC | MCW_PC );
72. test_fun();
74. return 0;
75. }

复制代码

无心人

宝宝： 减0.5不是个好办法，或许因为别的函数改变了RC， 就反而得不到正确的数值了 锅老大： 你这算偷懒的 在函数外 控制RC获得的性能提升有限的

liangbch

gxq 的这个版本最简单，但是在程序开头改变了FPU控制字，不知到在运行期间，是不是会影响的一些浮点数学函数。

无心人

2. double b32[] = {0.0,4294967296.0};
4. __declspec(naked)
5. DWORD __fastcall iSqrt_FPU(DWORD n)
6. {
7. __asm
8. {
9. push ecx
10. sub esp, 4
11. fnstcw word ptr [esp]
12. mov edx, dword ptr [esp]
13. or dword ptr [esp], 0x0C00
14. mov eax, ecx
15. and eax, 0x80000000
16. shr eax, 31
17. fld qword ptr [b32 + eax * 8]
18. fild dword ptr [esp + 4]
19. faddp st(1), st
20. fsqrt
21. fldcw word ptr [esp]
22. fistp dword ptr [esp + 4]
23. mov dword ptr [esp], edx
24. fldcw word ptr [esp]
25. add esp, 4
26. pop eax
27. ret
28. }
29. }
31. __declspec(naked)
32. DWORD iSqrt_FPU1(DWORD n)
33. {
34. __asm
35. {
36. sub esp, 4
37. fnstcw word ptr [esp]
38. mov edx, dword ptr [esp]
39. or dword ptr [esp], 0x0C00
40. mov eax, [esp + 8]
41. and eax, 0x80000000
42. shr eax, 31
43. fld qword ptr [b32 + eax * 8]
44. fild dword ptr [esp + 8]
45. faddp st(1), st
46. fsqrt
47. fldcw word ptr [esp]
48. fistp dword ptr [esp + 8]
49. mov dword ptr [esp], edx
50. fldcw word ptr [esp]
51. add esp, 4
52. mov eax, [esp + 4]
53. ret
54. }
55. }
57. __declspec(naked)
58. DWORD __fastcall iSqrt_FPU2(DWORD n)
59. {
60. __asm
61. {
62. push ecx
63. mov eax, ecx
64. and eax, 0x80000000
65. shr eax, 31
66. fld qword ptr [b32 + eax * 8]
67. fild dword ptr [esp]
68. faddp st(1), st
69. fsqrt
70. sub esp, 8
71. fstp qword ptr [esp]
72. mov edx, dword ptr [esp + 4]
73. mov eax, edx
74. and edx,0x7ff00000
75. and eax,0xfffff
76. shr edx, 20
77. or eax, 0x100000
78. xchg ecx, edx
79. sub ecx, 1043
80. neg ecx
81. shr eax, cl
82. xchg edx, ecx
83. add esp, 12
84. or ecx, ecx
85. cmove eax, ecx
86. ret
87. }
88. }

复制代码

比较时间 FPU2 < FPU < FPU1 另外，程序开始就改变RC的 会造成函数结果错误 输入小的时候测试不到

无心人

下面就是寻求执行时间少于60周期的整数算法 考虑乘法总延迟16周期 除法总延迟是50周期 所以 1、不能出现除法 2、乘法最多2个

liangbch

你的FPU2的版本写成 sub edx, 1022 mov ecx, 21 sub ecx, edx 我的则做了优化，直接写成以下形式，实际上是等价的 mov ecx,1043 sub ecx,edx

无心人

我不理解你的做法 所以改了啊

无心人

呵呵 在持续改进FPU2函数，增加了非跳转处理0的代码 另外，FWAIT似乎并不必须吧