二进制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
   
10.                
    
11.                 and  edx,ecx        //阶码
    
12.                 and  eax,ecx        //尾数
    
13. 
    
14.                 shr  edx,20                //得到阶段码
    
15.                 add  eax,0x100000 //设置位数最高有效位
    
16.                 mov  ecx,1043
    
17.                 sub  ecx,edx
    
18.                 shr  eax,cl                //得到整数
    
19.                
    
20. 
    
21.         }
    
22. }
    
23. 
    
24. void testdouble2int()
    
25. {
    
26.         double f;
    
27.         int n;
    
28. 
    
29.         f=23.4567;
    
30.         n=double2int(&f);
    
31.         printf("%lf=%d\n",f,n);
    
32. 
    
33.         f=23.9845;
    
34.         n=double2int(&f);
    
35.         printf("%lf=%d\n",f,n);
    
36. 
    
37. 
    
38.         f=2.4567;
    
39.         n=double2int(&f);
    
40.         printf("%lf=%d\n",f,n);
    
41.        
    
42.         f=2.9845;
    
43.         n=double2int(&f);
    
44.         printf("%lf=%d\n",f,n);
    
45. 
    
46. 
    
47.         f=65534.4567;
    
48.         n=double2int(&f);
    
49.         printf("%lf=%d\n",f,n);
    
50. 
    
51.         f=65534.9977;
    
52.         n=double2int(&f);
    
53.         printf("%lf=%d\n",f,n);
    
54.        
    
55. 
    
56.         f=1.00;
    
57.         n=double2int(&f);
    
58.         printf("%lf=%d\n",f,n);
    
59. 
    
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};
   
2. 
   
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. }
    
21. 
    
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. }
    
38. 
    
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等完全平方数时，计算结果错误。

1. 
   
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
    
21. 
    
22. fwait
    
23. mov  ecx,dword ptr [esp-4]
    
24. mov  edx,0xfff00000
    
25. mov  eax,0xfffff
    
26. 
    
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 //得到整数
    
34. 
    
35. pop ecx
    
36. ret
    
37. }
    
38. }

复制代码

无心人

老大啊

我已经修正了栈平衡问题
35#的代码是正确的
另外，这两个代码都需要定义个
double b32 [] = {0.0, 4294967296.0};
我早已经偷偷的把
double b32 = 4294967296.0;
改了啊

否则要用到FCMOVcc等复杂指令的