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

liangbch

我也利用54楼的思想写了一个函数。
   这个函数包含2种做法，当 定义了 LOOP_UNROLL， 使用swich case 技巧做了循环展开，速度应该更快。当LOOP_UNROLL未定义，使用常规的循环算法，楼主可否测试一下这个版本的速度。
  下面给出源代码。

1. 
   
2. 
   
3. typedef unsigned long DWORD;
   
4. typedef unsigned char BYTE;
   
5. 
   
6. inline DWORD log2(DWORD n)
   
7. {
   
8.         _asm
   
9.         {
   
10.                 mov ecx,n
    
11.                 bsr eax,ecx
    
12.         }
    
13. }
    
14. 
    
15. const unsigned char sqrtTab2[]=
    
16. {
    
17.         0,  1,  1,  1,  2,  2,  2,  2,
    
18.         2,  3,  3,  3,  3,  3,  3,  3,
    
19.         4,  4,  4,  4,  4,  4,  4,  4,
    
20.         4,  5,  5,  5,  5,  5,  5,  5,
    
21.         5,  5,  5,  5,  6,  6,  6,  6,
    
22.         6,  6,  6,  6,  6,  6,  6,  6,
    
23.         6,  7,  7,  7,  7,  7,  7,  7,
    
24.         7,  7,  7,  7,  7,  7,  7,  7,
    
25.         8,  8,  8,  8,  8,  8,  8,  8,
    
26.         8,  8,  8,  8,  8,  8,  8,  8,
    
27.         8,  9,  9,  9,  9,  9,  9,  9,
    
28.         9,  9,  9,  9,  9,  9,  9,  9,
    
29.         9,  9,  9,  9, 10, 10, 10, 10,
    
30.         10, 10, 10, 10, 10, 10, 10, 10,
    
31.         10, 10, 10, 10, 10, 10, 10, 10,
    
32.         10, 11, 11, 11, 11, 11, 11, 11,
    
33.         11, 11, 11, 11, 11, 11, 11, 11,
    
34.         11, 11, 11, 11, 11, 11, 11, 11,
    
35.         12, 12, 12, 12, 12, 12, 12, 12,
    
36.         12, 12, 12, 12, 12, 12, 12, 12,
    
37.         12, 12, 12, 12, 12, 12, 12, 12,
    
38.         12, 13, 13, 13, 13, 13, 13, 13,
    
39.         13, 13, 13, 13, 13, 13, 13, 13,
    
40.         13, 13, 13, 13, 13, 13, 13, 13,
    
41.         13, 13, 13, 13, 14, 14, 14, 14,
    
42.         14, 14, 14, 14, 14, 14, 14, 14,
    
43.         14, 14, 14, 14, 14, 14, 14, 14,
    
44.         14, 14, 14, 14, 14, 14, 14, 14,
    
45.         14, 15, 15, 15, 15, 15, 15, 15,
    
46.         15, 15, 15, 15, 15, 15, 15, 15,
    
47.         15, 15, 15, 15, 15, 15, 15, 15,
    
48.         15, 15, 15, 15, 15, 15, 15, 15
    
49. };
    
50. 
    
51. const unsigned char sqrTab2[]=
    
52. {
    
53.         0,        1,        4,        9,        16,        25,        36,        49,
    
54.         64,        81,        100,121,144,169,196,225
    
55. };
    
56. 
    
57. 
    
58. #define LOOP_UNROLL
    
59. extern "C"
    
60. DWORD __fastcall UintSqrt2(DWORD x)
    
61. {
    
62.         DWORD bc,m1,m2,n,t,i;
    
63. 
    
64.         if (x<256)
    
65.                 return sqrtTab2[x];
    
66.        
    
67.         bc=log2(x)/2;
    
68.         m2=(x>>(bc*2-6));
    
69.         n=sqrtTab2[m2];
    
70.         m1=sqrTab2[n];
    
71. 
    
72. #ifdef LOOP_UNROLL
    
73.         switch(bc)
    
74.         {
    
75.         case 15:
    
76.                 m2=  (x>>22);
    
77.                 n <<= 1;        m1<<= 2;
    
78.                 t = m1 + (n<<1)+1;
    
79.                 if ( t<=m2)        { m1=t; n++; }
    
80.         case 14:
    
81.                 m2=  (x>>20);
    
82.                 n <<= 1;        m1<<= 2;
    
83.                 t = m1 + (n<<1)+1;
    
84.                 if ( t<=m2)        { m1=t; n++; }
    
85.         case 13:
    
86.                 m2=  (x>>18);
    
87.                 n <<= 1;        m1<<= 2;
    
88.                 t = m1 + (n<<1)+1;
    
89.                 if ( t<=m2)        { m1=t; n++; }
    
90.         case 12:
    
91.                 m2=  (x>>16);
    
92.                 n <<= 1;        m1<<= 2;
    
93.                 t = m1 + (n<<1)+1;
    
94.                 if ( t<=m2)        { m1=t; n++; }
    
95.         case 11:
    
96.                 m2=  (x>>14);
    
97.                 n <<= 1;        m1<<= 2;
    
98.                 t = m1 + (n<<1)+1;
    
99.                 if ( t<=m2)        { m1=t; n++; }
    
100.         case 10:
     
101.                 m2=  (x>>12);
     
102.                 n <<= 1;        m1<<= 2;
     
103.                 t = m1 + (n<<1)+1;
     
104.                 if ( t<=m2)        { m1=t; n++; }
     
105.         case 9:
     
106.                 m2=  (x>>10);
     
107.                 n <<= 1;        m1<<= 2;
     
108.                 t = m1 + (n<<1)+1;
     
109.                 if ( t<=m2)        { m1=t; n++; }
     
110.         case 8:
     
111.                 m2=  (x>>8);
     
112.                 n <<= 1;        m1<<= 2;
     
113.                 t = m1 + (n<<1)+1;
     
114.                 if ( t<=m2)        { m1=t; n++; }
     
115.         case 7:
     
116.                 m2=  (x>>6);
     
117.                 n <<= 1;        m1<<= 2;
     
118.                 t = m1 + (n<<1)+1;
     
119.                 if ( t<=m2)        { m1=t; n++; }
     
120.         case 6:
     
121.                 m2=  (x>>4);
     
122.                 n <<= 1;        m1<<= 2;
     
123.                 t = m1 + (n<<1)+1;
     
124.                 if ( t<=m2)        { m1=t; n++; }
     
125.         case 5:
     
126.                 m2=  (x>>2);
     
127.                 n <<= 1;        m1<<= 2;
     
128.                 t = m1 + (n<<1)+1;
     
129.                 if ( t<=m2)        { m1=t; n++; }
     
130.         case 4:
     
131.                 m2=  x;
     
132.                 n <<= 1;        m1<<= 2;
     
133.                 t = m1 + (n<<1)+1;
     
134.                 if ( t<=m2)        { m1=t; n++; }
     
135.         case 3:
     
136.         case 2:
     
137.         case 1:
     
138.         case 0:;
     
139.         }
     
140. #else
     
141.         for (i=bc;i>=4;i--)
     
142.         {
     
143.                 m2=(x>>(i*2-8));
     
144.                 n  <<= 1;        m1 <<= 2;
     
145.                 t= m1 + (n<<1)+1;
     
146.                 if ( t<=m2)
     
147.                 {        m1=t; n++;        }
     
148.         }
     
149. #endif
     
150.         return n;
     
151. }

复制代码

无心人

GxQ的代码
好像，并不快啊

无心人

GxQ代码的汇编版本
提高50%左右性能

1. 
   
2. _declspec(naked)
   
3. UINT32 __fastcall iSqrt_GxQ_asm( UINT32 n )
   
4. {
   
5. /*
   
6.     UINT32 t;
   
7.     UINT32 ret = 0;
   
8. 
   
9.     #define SQRT_CORE(x)        \
   
10.         t = ret + (1UL<<(x));   \
    
11.         ret >>= 1;              \
    
12.         if ( n >= t )           \
    
13.         {                       \
    
14.             n -= t;             \
    
15.             ret += (1UL<<(x));  \
    
16.         }
    
17. 
    
18.     SQRT_CORE(30, 28, ..., ..., 0);
    
19. */
    
20. __asm
    
21. {
    
22. push ebx  //t
    
23. push esi
    
24. push edi
    
25. xor eax, eax //ret
    
26. mov ebx, eax  
    
27. mov edx, 1
    
28. shl edx, 30
    
29. //30
    
30. add ebx, edx
    
31. shr eax, 1
    
32.         mov esi, ecx
    
33. sub esi, ebx
    
34. mov edi, eax
    
35. add edi, edx
    
36.         cmp ecx, ebx
    
37.         cmovae ecx, esi
    
38. cmovae eax, edi
    
39.         //28
    
40. shr edx, 2
    
41. mov ebx, eax
    
42. add ebx, edx
    
43. shr eax, 1
    
44.         mov esi, ecx
    
45. sub esi, ebx
    
46. mov edi, eax
    
47. add edi, edx
    
48.         cmp ecx, ebx
    
49.         cmovae ecx, esi
    
50. cmovae eax, edi
    
51.         //26
    
52. shr edx, 2
    
53. mov ebx, eax
    
54. add ebx, edx
    
55. shr eax, 1
    
56.         mov esi, ecx
    
57. sub esi, ebx
    
58. mov edi, eax
    
59. add edi, edx
    
60.         cmp ecx, ebx
    
61.         cmovae ecx, esi
    
62. cmovae eax, edi
    
63.         //24
    
64. shr edx, 2
    
65. mov ebx, eax
    
66. add ebx, edx
    
67. shr eax, 1
    
68.         mov esi, ecx
    
69. sub esi, ebx
    
70. mov edi, eax
    
71. add edi, edx
    
72.         cmp ecx, ebx
    
73.         cmovae ecx, esi
    
74. cmovae eax, edi
    
75.         //22
    
76. shr edx, 2
    
77. mov ebx, eax
    
78. add ebx, edx
    
79. shr eax, 1
    
80.         mov esi, ecx
    
81. sub esi, ebx
    
82. mov edi, eax
    
83. add edi, edx
    
84.         cmp ecx, ebx
    
85.         cmovae ecx, esi
    
86. cmovae eax, edi
    
87.         //20
    
88. shr edx, 2
    
89. mov ebx, eax
    
90. add ebx, edx
    
91. shr eax, 1
    
92.         mov esi, ecx
    
93. sub esi, ebx
    
94. mov edi, eax
    
95. add edi, edx
    
96.         cmp ecx, ebx
    
97.         cmovae ecx, esi
    
98. cmovae eax, edi
    
99.         //18
    
100. shr edx, 2
     
101. mov ebx, eax
     
102. add ebx, edx
     
103. shr eax, 1
     
104.         mov esi, ecx
     
105. sub esi, ebx
     
106. mov edi, eax
     
107. add edi, edx
     
108.         cmp ecx, ebx
     
109.         cmovae ecx, esi
     
110. cmovae eax, edi
     
111.         //16
     
112. shr edx, 2
     
113. mov ebx, eax
     
114. add ebx, edx
     
115. shr eax, 1
     
116.         mov esi, ecx
     
117. sub esi, ebx
     
118. mov edi, eax
     
119. add edi, edx
     
120.         cmp ecx, ebx
     
121.         cmovae ecx, esi
     
122. cmovae eax, edi
     
123.         //14
     
124. shr edx, 2
     
125. mov ebx, eax
     
126. add ebx, edx
     
127. shr eax, 1
     
128.         mov esi, ecx
     
129. sub esi, ebx
     
130. mov edi, eax
     
131. add edi, edx
     
132.         cmp ecx, ebx
     
133.         cmovae ecx, esi
     
134. cmovae eax, edi
     
135.         //12
     
136. shr edx, 2
     
137. mov ebx, eax
     
138. add ebx, edx
     
139. shr eax, 1
     
140.         mov esi, ecx
     
141. sub esi, ebx
     
142. mov edi, eax
     
143. add edi, edx
     
144.         cmp ecx, ebx
     
145.         cmovae ecx, esi
     
146. cmovae eax, edi
     
147.         //10
     
148. shr edx, 2
     
149. mov ebx, eax
     
150. add ebx, edx
     
151. shr eax, 1
     
152.         mov esi, ecx
     
153. sub esi, ebx
     
154. mov edi, eax
     
155. add edi, edx
     
156.         cmp ecx, ebx
     
157.         cmovae ecx, esi
     
158. cmovae eax, edi
     
159.         //8
     
160. shr edx, 2
     
161. mov ebx, eax
     
162. add ebx, edx
     
163. shr eax, 1
     
164.         mov esi, ecx
     
165. sub esi, ebx
     
166. mov edi, eax
     
167. add edi, edx
     
168.         cmp ecx, ebx
     
169.         cmovae ecx, esi
     
170. cmovae eax, edi
     
171.         //6
     
172. shr edx, 2
     
173. mov ebx, eax
     
174. add ebx, edx
     
175. shr eax, 1
     
176.         mov esi, ecx
     
177. sub esi, ebx
     
178. mov edi, eax
     
179. add edi, edx
     
180.         cmp ecx, ebx
     
181.         cmovae ecx, esi
     
182. cmovae eax, edi
     
183.         //4
     
184. shr edx, 2
     
185. mov ebx, eax
     
186. add ebx, edx
     
187. shr eax, 1
     
188.         mov esi, ecx
     
189. sub esi, ebx
     
190. mov edi, eax
     
191. add edi, edx
     
192.         cmp ecx, ebx
     
193.         cmovae ecx, esi
     
194. cmovae eax, edi
     
195.         //2
     
196. shr edx, 2
     
197. mov ebx, eax
     
198. add ebx, edx
     
199. shr eax, 1
     
200.         mov esi, ecx
     
201. sub esi, ebx
     
202. mov edi, eax
     
203. add edi, edx
     
204.         cmp ecx, ebx
     
205.         cmovae ecx, esi
     
206. cmovae eax, edi
     
207.         //0
     
208. shr edx, 2
     
209. mov ebx, eax
     
210. add ebx, edx
     
211. shr eax, 1
     
212.         mov esi, ecx
     
213. sub esi, ebx
     
214. mov edi, eax
     
215. add edi, edx
     
216.         cmp ecx, ebx
     
217.         cmovae ecx, esi
     
218. cmovae eax, edi
     
219. //ret
     
220.         pop edi
     
221. pop esi
     
222. pop ebx
     
223.         ret
     
224. 
     
225. }
     
226. }
     

复制代码

FPU2 Version: 3589.069 ms
iSqrt_16 Version: 26237.479 ms
iSqrt_GxQ Version: 16068.268 ms
iSqrt_GxQ_asm Version: 7884.029 ms

无心人

FPU2 Version: 3596.297 ms
iSqrt_16 Version: 26235.439 ms
iSqrt_GxQ Version: 16177.534 ms
iSqrt_GxQ_asm Version: 7868.024 ms
UintSqrt2 Version: 21584.053 ms
=========================================
急需各位机器的测试结果

我这里Core 2 XEON 1.6似乎还是浮点方式快的多

无心人

我的代码参考46
iSqrt_16在55
iSqrt_GxQ在56
iSqrtGxQ_在63
UintSqrt2在62

无心人

测试方法
1、从1循环到2^32-1
2、测试2100000000到2200000000

liangbch

这个主题中，出现了各种不同版本开平方版本，我想对这些版本全部测一遍。

1.首先由作者提出需要测试的程序，在哪一楼，那个函数。
2.某个人负责将这些程序收集起来，做成一个测试程序，打包上传。
3.任何人就可以下载编译并测试了，测试完毕提交测试环境和测试结果。

2点说明：
1. 测试程序应包括正确性测试和性能测试。
2  我的3个整数版本中，运行速度和被开放数的大小有关。 所以我建议，在做性能测试时，增加一个case,测试n=1到65535时的速度。

我先开个头。
我总共编写了6个有效的程序，其中包含一个汇编版，我需要测试的程序有
1.  11楼的 C 整数版
2.  42楼的3个浮点汇编版
3.  61楼的C 整数非牛顿迭代版。

无心人

我需要测试46， 55， 63

这个测试很耗费时间

我的机器测试46三个版本是160秒左右
55版本是1200秒左右
63版本是330秒
56大概是670秒

liangbch

代码整理中，请稍候。。。

liangbch

全部函数和测试程序源码包。