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

liangbch

应该是，用浮点指令可以做的更好。 刚刚完成了一个完全用整数指令将浮点数转化为整数的程序，见下，这个函数唯一不足之处是专为0将出错。

1. //*p必须是整数或者0
2. int double2int(double *p)
3. {
4. _asm
5. {
6. mov eax,p
7. mov ecx,dword ptr [eax+4]
8. mov edx,0xfff00000
9. mov eax,0xfffff
11. and edx,ecx //阶码
12. and eax,ecx //尾数
14. shr edx,20 //得到阶段码
15. add eax,0x100000 //设置位数最高有效位
16. mov ecx,1043
17. sub ecx,edx
18. shr eax,cl //得到整数
21. }
22. }
24. void testdouble2int()
25. {
26. double f;
27. int n;
29. f=23.4567;
30. n=double2int(&f);
31. printf("%lf=%d\n",f,n);
33. f=23.9845;
34. n=double2int(&f);
35. printf("%lf=%d\n",f,n);
38. f=2.4567;
39. n=double2int(&f);
40. printf("%lf=%d\n",f,n);
42. f=2.9845;
43. n=double2int(&f);
44. printf("%lf=%d\n",f,n);
47. f=65534.4567;
48. n=double2int(&f);
49. printf("%lf=%d\n",f,n);
51. f=65534.9977;
52. n=double2int(&f);
53. printf("%lf=%d\n",f,n);
56. f=1.00;
57. n=double2int(&f);
58. printf("%lf=%d\n",f,n);
60. f=0;
61. n=double2int(&f);
62. printf("%lf=%d\n",f,n);
63. }

复制代码

HouSisong 写的类似的函数是这样的， 见http://blog.csdn.net/housisong/archive/2007/05/19/1616026.aspx

1. inline long _ftol_ieee(float f)
2. {
3. long a = *(long*)(&f);
4. unsigned long mantissa = (a&((1<<23)-1))|(1<<23); //不支持非规格化浮点数
5. long exponent = ((a&0x7fffffff)>>23);
6. long r = (mantissa<<8) >> (31+127-exponent);
7. long sign = (a>>31);
8. return ((r ^ (sign)) - sign ) &~ ((exponent-127)>>31);
9. }

复制代码

[ 本帖最后由 liangbch 于 2009-2-12 12:30 编辑 ]

无心人

早上想是否32位浮点足够，中午发现，可能存在精度丢失 实际计算下 Prelude> let b = 2.0**32 - 1.0 Prelude> b 4.294967295e9 Prelude> let c = sqrt b Prelude> c 65535.999992370605 Prelude> let d = 255.0 / 256.0 Prelude> d 0.99609375 无法满足要求 ======================== 我也想用整数处理浮点截断 呵呵 看来只能考虑64位浮点了 下午调试好另外三个就考虑下

无心人

上面的帖子，我提到了 用修改RC的方法需要多9条指令 且存在16位计算 你的代码是11条 很可能比fistp快 可以考虑 =========================== 不对 需要额外的两个压栈指令 和存储 即 push 0 push 0 fst qword ptr [esp]

无心人

64位浮点 0的话，指数和尾数都是0 1是指数乃1023，尾数是0 看是否有避免跳转判断的代码

无心人

目前的三个版本最终代码

1. double b32[] = {0.0, 4294967296.0};
3. __declspec(naked)
4. DWORD __fastcall iSqrt_SSE3(DWORD n)
5. {
6. __asm
7. {
8. push ecx
9. mov eax, ecx
10. and eax, 0x80000000
11. shr eax, 31
12. fld qword ptr [b32 + eax * 8]
13. fild dword ptr [esp]
14. faddp st(1), st
15. fsqrt
16. fisttp dword ptr [esp]
17. pop eax
18. ret
19. }
20. }
22. __declspec(naked)
23. DWORD __fastcall iSqrt_SSE2(DWORD n)
24. {
25. __asm
26. {
27. cvtsi2sd xmm0, ecx
28. mov eax, ecx
29. and eax, 0x80000000
30. shr eax, 31
31. movsd xmm1, qword ptr [b32 + eax * 8]
32. addsd xmm0, xmm1
33. sqrtsd xmm0, xmm0
34. cvttsd2si eax, xmm0
35. ret
36. }
37. }
39. __declspec(naked)
40. DWORD __fastcall iSqrt_FPU(DWORD n)
41. {
42. __asm
43. {
44. push ecx
45. sub esp, 4
46. fnstcw word ptr [esp]
47. mov edx, dword ptr [esp]
48. or dword ptr [esp], 0x0C00
49. mov eax, ecx
50. and eax, 0x80000000
51. shr eax, 31
52. fld qword ptr [b32 + eax * 8]
53. fild dword ptr [esp + 4]
54. faddp st(1), st
55. fsqrt
56. fldcw word ptr [esp]
57. fistp dword ptr [esp + 4]
58. mov dword ptr [esp], edx
59. fldcw word ptr [esp]
60. add esp, 4
61. pop eax
62. ret
63. }
64. }

复制代码

无心人

Test FPU: sqrt(1023) = 31 FPU: sqrt(4000000000) = 63245 SSE2: sqrt(1023) = 31 SSE2: sqrt(4000000000) = 63245 SSE3: sqrt(1023) = 31 SSE3: sqrt(4000000000) = 63245 FPU Version: 158930.292 ms SSE2 Version: 177428.864 ms SSE3 Version: 155667.881 ms ======================== 从1到2^32-1循环

无心人

可以看到 FPU版的控制 RC方式仅比fisttp的SSE3指令略微多一点 计算指令周期 FPU: 59.2 clock SSE2: 66.1 clock SSE3: 58.0 clock

liangbch

你的35楼的 FPU 版本在我这里无法运行，当执行到 fld qword ptr [b32 + eax * 8] 时，程序崩溃。 另外，你之前写的FPU版本也有问题， FPU 入栈操作多余出栈，最后当使用fld 时因栈满，到导致st(0)变为一个无效值。 另外，你24楼的程序也有问题，在我这里运行时，当执行到 movsd xmm1, qword ptr [b32 + eax * 4] ，程序崩溃。 下面给出一个FPU的版本，这个程序做到栈平衡，但当被开放数是1,4,9,16,25,49等完全平方数时，计算结果错误。

2. double b32 = 4294967296.0;
3. double zero5= -0.5;
4. __declspec(naked)
5. DWORD __fastcall iSqrt_FPU1(DWORD n)
6. {
7. __asm
8. {
9. push ecx
10. fld qword ptr [b32]
11. fldz
12. cmp ecx, 0x8000000
13. fcmovnb st, st(1)
14. fild dword ptr [esp]
15. faddp st(1),st
16. fsqrt
17. fadd qword ptr [zero5]
18. fistp dword ptr [esp]
19. fstp st
20. pop eax
21. ret
22. }
23. }

复制代码

liangbch

在贴一个FPU的版本，完全不是用 SSE/SSE2指令，可正常工作。

1. __declspec(naked)
2. DWORD __fastcall iSqrt_FPU2(DWORD n)
3. {
4. __asm
5. {
6. or ecx,ecx
7. jnz next00
8. mov eax,0
9. ret
10. next00:
11. push ecx
12. fld qword ptr [b32]
13. fldz
14. cmp ecx, 0x80000000
15. fcmovnb st, st(1)
16. fild dword ptr [esp]
17. faddp st(1),st
18. fsqrt
19. fstp qword ptr [esp-8]
20. fstp st
22. fwait
23. mov ecx,dword ptr [esp-4]
24. mov edx,0xfff00000
25. mov eax,0xfffff
27. and edx,ecx //阶码
28. and eax,ecx //尾数
29. shr edx,20 //得到阶段码
30. add eax,0x100000 //设置位数最高有效位
31. mov ecx,1043
32. sub ecx,edx
33. shr eax,cl //得到整数
35. pop ecx
36. ret
37. }
38. }

复制代码

无心人

老大啊 我已经修正了栈平衡问题 35#的代码是正确的 另外，这两个代码都需要定义个 double b32 [] = {0.0, 4294967296.0}; 我早已经偷偷的把 double b32 = 4294967296.0; 改了啊 否则要用到FCMOVcc等复杂指令的