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

无心人 · 发表于 2009-2-11 17:43:27

__declspec(naked)
DWORD __fastcall iSqrt_SSE41(DWORD n)
{
__asm
{
mov eax, ecx
and eax, 0x80000000
movd xmm0, eax
psllq xmm0, 32
mov eax, 0
cvtsi2sd xmm1, eax
movsd xmm2, qword ptr [b32]
blendvpd xmm1, xmm2, xmm0
cvtsi2sd xmm0, ecx
addsd xmm0, xmm1
sqrtsd xmm0, xmm0
cvttsd2si eax, xmm0
ret
}
}

复制代码

SSE2的无跳转版本不如带跳转的
只好写SSE4.1的代码了
可惜，我查我服务器，只支持到SSSE3
哎
谁可以帮测试下代码正确性？
================================
由于改进了SSE2版本
所以此版本作废

liangbch · 发表于 2009-2-11 19:06:18

fisttp 是SSE3指令，在我的VC6.0下不能通过编译，即使通过编译也不能运行。我的电脑是PIV2.6.
以下的内容转自：http://blog.csdn.net/housisong/archive/2007/05/19/1616026.aspx

SSE3增加了一条FPU取整指令fisttp,和fistp指令功能几乎相同(我的电脑上经过测试速度也相同)，但默认向0取整，和RC场设置无关，所以使用fisttp的代码就可以不管RC场了，有利于简化代码和优化性能

在上文中，作者还是先了使用 fistp 指令的浮点转化整数的版本，无心人可否尝试一下。看看速度如何，另外，这样的版本也可以在我的PIV电脑上测试。

无心人 · 发表于 2009-2-11 19:20:03

呵呵

我没注意，我还以为是P6以后都支持呢

另外，下午的程序写糊涂了

并没这么复杂

无心人 · 发表于 2009-2-11 19:23:33

double b32[] = {0.0, 4294967296.0};
__declspec(naked)
DWORD __fastcall iSqrt_SSE2(DWORD n)
{
__asm
{
mov eax, ecx
and eax, 0x80000000
shr eax, 31
movsd xmm1, qword ptr [b32 + eax * 8]
cvtsi2sd xmm0, ecx
addsd xmm0, xmm1
sqrtsd xmm0, xmm0
cvttsd2si eax, xmm0
ret
}
}

复制代码

这么做就能避免跳转了
至于fisttp，我想控制RC可能有点复杂，呵呵

liangbch · 发表于 2009-2-11 19:40:08

测试了3个版本：UintSqrt，iSqrt_FPU 和 sqrt_ref, 在我的电脑上运行2^28(从1到2&28-1)次,其速度如下

  iSqrt_FPU: 4375 ms
  UintSqrt:  8406 ms
  sqrt_ref: 16328 ms

说明:
1. UintSqrt为我编写的使用bsr的c语言版本,略有修改
2. sqrt_ref的代码如下:
DWORD sqrt_ref(DWORD n)
{
DWORD r=(DWORD)(sqrt(double(n)));
return r;
}

3.iSqrt_FPU 函数修改了一条指令,将FISTTP 改为FISTP

  综上所述,
1. 系统提供的函数最慢.
2. iSqrt_FPU 目前是个错误的版本,将其修正,使其在PIV上正确运行,需要增加几条指令,将会使得速度变慢.
3. 如果UintSqrt改为汇编优化,速度应可以提高,但提上的幅度应该不会太大.在2^28以内是否会超过修正后的iSqrt_FPU仍是个未知数.

另外，限于我的VC环境，你的使用 SSE2指令的版本在我的电脑上没有编译，也就没有办法测试了。

无心人 · 发表于 2009-2-11 19:41:12

double b32[] = {0.0, 4294967296.0};
__declspec(naked)
DWORD __fastcall iSqrt_FPU_alpha(DWORD n)
{
__asm
{
push ecx
mov eax, ecx
and eax, 0x80000000
fld qword ptr [b32 + eax * 8]
fild dword ptr [esp]
faddp st(1), st
fsqrt
push 0
fnstcw word ptr [esp]
mov edx, dword ptr [esp]
or dword ptr [esp], 0x0C00
fldcw word ptr [esp]
fistp dword ptr [esp + 4]
mov dword ptr [esp], edx
fldcw word ptr [esp]
add esp, 4
pop eax
ret
}
}

复制代码

不用fisttp
多了9条指令
因为无法测试，特取名...FPU_alpha
明天有时间测试

不用__fastcall似乎能节约下ecx和一条指令
================================
似乎不用__fastcall也无法避免堆栈上占用两个双字
所以，还是用__fastcall

liangbch · 发表于 2009-2-11 23:29:53

终于将C版改为汇编语言版本了，测试发现，速度仅仅提高约10%，计算2^28次平方根，用时从 7902.3ms 提高到7192.7 ms。

函数原型：
extern "C" DWORD __fastcall UintSqrt(DWORD n);

汇编源代码：

.486P
.387
OPTION DOTNAME
_TEXT SEGMENT DWORD PUBLIC FLAT 'CODE'
ALIGN 004H
_TEXT ENDS
_DATA SEGMENT DWORD PUBLIC FLAT 'DATA'
ALIGN 004H
_DATA ENDS
_BSS SEGMENT DWORD PUBLIC FLAT 'BSS'
ALIGN 004H
_BSS ENDS
_RDATA SEGMENT DWORD PUBLIC FLAT 'DATA'
ALIGN 004H
_RDATA ENDS
_TEXT1 SEGMENT DWORD PUBLIC FLAT 'CODE'
ALIGN 004H
_TEXT1 ENDS
_DATA1 SEGMENT DWORD PUBLIC FLAT 'DATA'
ALIGN 004H
_DATA1 ENDS
ASSUME CS:FLAT,DS:FLAT,SS:FLAT
_DATA SEGMENT DWORD PUBLIC FLAT 'DATA'
ALIGN 4
_sqrtTab DB 00H
DB 01H
DB 01H
DB 01H
DB 02H
DB 02H
DB 02H
DB 02H
DB 02H
DB 03H
DB 03H
DB 03H
DB 03H
DB 03H
DB 03H
DB 03H
DB 04H
DB 04H
DB 04H
DB 04H
DB 04H
DB 04H
DB 04H
DB 04H
DB 04H
DB 05H
DB 05H
DB 05H
DB 05H
DB 05H
DB 05H
DB 05H
DB 05H
DB 05H
DB 05H
DB 05H
DB 06H
DB 06H
DB 06H
DB 06H
DB 06H
DB 06H
DB 06H
DB 06H
DB 06H
DB 06H
DB 06H
DB 06H
DB 06H
DB 07H
DB 07H
DB 07H
DB 07H
DB 07H
DB 07H
DB 07H
DB 07H
DB 07H
DB 07H
DB 07H
DB 07H
DB 07H
DB 07H
DB 07H
DB 080H
DB 080H
DB 081H
DB 082H
DB 083H
DB 084H
DB 085H
DB 086H
DB 087H
DB 088H
DB 089H
DB 08aH
DB 08bH
DB 08cH
DB 08dH
DB 08eH
DB 08fH
DB 090H
DB 090H
DB 091H
DB 092H
DB 093H
DB 094H
DB 095H
DB 096H
DB 096H
DB 097H
DB 098H
DB 099H
DB 09aH
DB 09bH
DB 09bH
DB 09cH
DB 09dH
DB 09eH
DB 09fH
DB 0a0H
DB 0a0H
DB 0a1H
DB 0a2H
DB 0a3H
DB 0a3H
DB 0a4H
DB 0a5H
DB 0a6H
DB 0a7H
DB 0a7H
DB 0a8H
DB 0a9H
DB 0aaH
DB 0aaH
DB 0abH
DB 0acH
DB 0adH
DB 0adH
DB 0aeH
DB 0afH
DB 0b0H
DB 0b0H
DB 0b1H
DB 0b2H
DB 0b2H
DB 0b3H
DB 0b4H
DB 0b5H
DB 0b5H
DB 0b6H
DB 0b7H
DB 0b7H
DB 0b8H
DB 0b9H
DB 0b9H
DB 0baH
DB 0bbH
DB 0bbH
DB 0bcH
DB 0bdH
DB 0bdH
DB 0beH
DB 0bfH
DB 0c0H
DB 0c0H
DB 0c1H
DB 0c1H
DB 0c2H
DB 0c3H
DB 0c3H
DB 0c4H
DB 0c5H
DB 0c5H
DB 0c6H
DB 0c7H
DB 0c7H
DB 0c8H
DB 0c9H
DB 0c9H
DB 0caH
DB 0cbH
DB 0cbH
DB 0ccH
DB 0ccH
DB 0cdH
DB 0ceH
DB 0ceH
DB 0cfH
DB 0d0H
DB 0d0H
DB 0d1H
DB 0d1H
DB 0d2H
DB 0d3H
DB 0d3H
DB 0d4H
DB 0d4H
DB 0d5H
DB 0d6H
DB 0d6H
DB 0d7H
DB 0d7H
DB 0d8H
DB 0d9H
DB 0d9H
DB 0daH
DB 0daH
DB 0dbH
DB 0dbH
DB 0dcH
DB 0ddH
DB 0ddH
DB 0deH
DB 0deH
DB 0dfH
DB 0e0H
DB 0e0H
DB 0e1H
DB 0e1H
DB 0e2H
DB 0e2H
DB 0e3H
DB 0e3H
DB 0e4H
DB 0e5H
DB 0e5H
DB 0e6H
DB 0e6H
DB 0e7H
DB 0e7H
DB 0e8H
DB 0e8H
DB 0e9H
DB 0eaH
DB 0eaH
DB 0ebH
DB 0ebH
DB 0ecH
DB 0ecH
DB 0edH
DB 0edH
DB 0eeH
DB 0eeH
DB 0efH
DB 0f0H
DB 0f0H
DB 0f1H
DB 0f1H
DB 0f2H
DB 0f2H
DB 0f3H
DB 0f3H
DB 0f4H
DB 0f4H
DB 0f5H
DB 0f5H
DB 0f6H
DB 0f6H
DB 0f7H
DB 0f7H
DB 0f8H
DB 0f8H
DB 0f9H
DB 0f9H
DB 0faH
DB 0faH
DB 0fbH
DB 0fbH
DB 0fcH
DB 0fcH
DB 0fdH
DB 0fdH
DB 0feH
DB 0feH
DB 0ffH
_DATA ENDS
_TEXT SEGMENT DWORD PUBLIC FLAT 'CODE'
ALIGN 4
PUBLIC _sqrtTab
PUBLIC @UintSqrt@4
@UintSqrt@4 PROC NEAR ; COMDAT
; _n = ecx
; if (n<64)
xor eax,eax
cmp ecx, 64
jae SHORT L25237
mov al, BYTE PTR _sqrtTab[ecx]
ret 0
L25237:
push ebx
; bc=log2(n)/2;
bsr eax, ecx
shr eax, 1
jmp DWORD PTR L25370[eax*4]
L25242:
; case 3:
; r=(sqrtTab[n]>>4);
mov al, BYTE PTR _sqrtTab[ecx]
shr eax, 4
; r1=r+1;
; if (r1*r1>=n) r++;
lea ebx, DWORD PTR [eax+1]
mov edx, ebx
imul edx, ebx
cmp ecx, edx
sbb eax,-1
pop ebx
ret 0
L25244:
; case 4:
; r=(sqrtTab[n>>2]>>3);
mov ebx, ecx
shr ebx, 2
mov al, BYTE PTR _sqrtTab[ebx]
shr eax, 3
; r1=r+1;
; if (r1*r1>=n) r++;
lea ebx, DWORD PTR [eax+1]
mov edx, ebx
imul edx, ebx
cmp ecx, edx
sbb eax,-1
pop ebx
ret 0
L25246:
; case 5:
; r=(sqrtTab[n>>4]>>2);
mov ebx, ecx
shr ebx, 4
mov al, BYTE PTR _sqrtTab[ebx]
shr eax, 2
; r1=r+1;
; if (r1*r1>=n) r++;
lea ebx, DWORD PTR [eax+1]
mov edx, ebx
imul edx, ebx
cmp ecx, edx
sbb eax,-1
pop ebx
ret 0
L25248:
; case 6:
; r=(sqrtTab[n>>6]>>1);
mov ebx, ecx
shr ebx, 6
mov al, BYTE PTR _sqrtTab[ebx]
shr eax, 1
; r1=r+1;
; if (r1*r1>=n) r++;
lea ebx, DWORD PTR [eax+1]
mov edx, ebx
imul edx, ebx
cmp ecx, edx
sbb eax,-1
pop ebx
ret 0
L25250:
; case 7:
; r=(sqrtTab[n>>8]);
mov ebx, ecx
shr ebx, 8
mov al, BYTE PTR _sqrtTab[ebx]
; r1=r+1;
; if (r1*r1>=n)
lea ebx, DWORD PTR [eax+1]
mov edx, ebx
imul edx, ebx
cmp ecx, edx
sbb eax,-1
pop ebx
ret 0
L25252:
; case 8:
; r=(sqrtTab[n>>10]<<1);
mov eax,ecx
xor ebx,ebx
shr eax, 10
mov bl, BYTE PTR _sqrtTab[eax]
add ebx, ebx ;r=(sqrtTab[n>>10]<<1);
;r= (n/r + r)/2;
xor edx,edx
mov eax,ecx
div ebx
add eax, ebx
shr eax, 1
; if (r*r>n) r--;
mov edx,eax
imul edx,eax
cmp ecx,edx
sbb eax,0
pop ebx
ret 0
L25254:
; case 9:
; r=(sqrtTab[n>>12]<<2);
mov eax,ecx
xor ebx,ebx
shr eax, 12
mov bl, BYTE PTR _sqrtTab[eax]
shl ebx, 2
;r= (n/r + r)/2;
xor edx,edx
mov eax,ecx
div ebx
add eax, ebx
shr eax, 1
; if (r*r>n) r--;
mov edx,eax
imul edx,eax
cmp ecx,edx
sbb eax,0
pop ebx
ret 0
L25256:
; case 10:
; r=(sqrtTab[n>>14]<<3);
mov eax,ecx
xor ebx,ebx
shr eax, 14
mov bl, BYTE PTR _sqrtTab[eax]
shl ebx, 3
;r= (n/r + r)/2;
xor edx,edx
mov eax,ecx
div ebx
add eax, ebx
shr eax, 1
; if (r*r>n) r--;
mov edx,eax
imul edx,eax
cmp ecx,edx
sbb eax,0
pop ebx
ret 0
L25258:
; case 11:
; r=(sqrtTab[n>>16]<<4);
mov eax,ecx
xor ebx,ebx
shr eax, 16
mov bl, BYTE PTR _sqrtTab[eax]
shl ebx, 4
;r= (n/r + r)/2;
xor edx,edx
mov eax,ecx
div ebx
add eax, ebx
shr eax, 1
; if (r*r>n) r--;
mov edx,eax
imul edx,eax
cmp ecx,edx
sbb eax,0
pop ebx
ret 0
L25260:
; case 12:
; r=(sqrtTab[n>>18]<<5);
mov eax,ecx
xor ebx,ebx
shr eax, 18
mov bl, BYTE PTR _sqrtTab[eax]
shl ebx, 5
;r= (n/r + r)/2;
xor edx,edx
mov eax,ecx
div ebx
add eax, ebx
shr eax, 1
; if (r*r>n) r--;
mov edx,eax
imul edx,eax
cmp ecx,edx
sbb eax,0
pop ebx
ret 0
L25262:
; case 13:
; r=(sqrtTab[n>>20]<<6);
mov eax,ecx
xor ebx,ebx
shr eax, 20
mov bl, BYTE PTR _sqrtTab[eax]
shl ebx, 6
;r= (n/r + r)/2;
xor edx,edx
mov eax,ecx
div ebx
add eax, ebx
shr eax, 1
; if (r*r>n) r--;
mov edx,eax
imul edx,eax
cmp ecx,edx
sbb eax,0
pop ebx
ret 0
L25264:
; case 14:
; r=(sqrtTab[n>>22]<<7);
mov eax,ecx
xor ebx,ebx
shr eax, 22
mov bl, BYTE PTR _sqrtTab[eax]
shl ebx, 7
;r= (n/r + r)/2;
xor edx,edx
mov eax,ecx
div ebx
add ebx, eax
shr ebx, 1
;r= (n/r + r)/2;
xor edx,edx
mov eax,ecx
div ebx
add eax, ebx
shr eax, 1
; if (r*r>n) r--;
mov edx,eax
imul edx,eax
cmp ecx,edx
sbb eax,0
pop ebx
ret 0
L25266:
; case 15:
; r=(sqrtTab[n>>24]<<8);
mov eax,ecx
xor ebx,ebx
shr eax, 24
mov bh, BYTE PTR _sqrtTab[eax]
;r= (n/r + r)/2;
xor edx,edx
mov eax,ecx
div ebx
add ebx, eax
shr ebx, 1
;r= (n/r + r)/2;
xor edx,edx
mov eax,ecx
div ebx
add eax, ebx
shr eax, 1
; if (r*r>n) r--;
mov edx,eax
imul edx,eax
cmp ecx,edx
sbb eax,0
pop ebx
ret 0
align 4
L25370:
DD 0
DD 0
DD 0
DD L25242
DD L25244
DD L25246
DD L25248
DD L25250
DD L25252
DD L25254
DD L25256
DD L25258
DD L25260
DD L25262
DD L25264
DD L25266
@UintSqrt@4 ENDP
_TEXT ENDS
END

复制代码

无心人 · 发表于 2009-2-12 09:08:10

你只是把C改成了汇编形式

用汇编要用汇编的思维重新写算法的

liangbch · 发表于 2009-2-12 11:02:28

No,No,No.
  虽然只是修改，但是修改的很多，尽量使用汇编的特性，以case 8为例，主要区别有：
  1. 寄存器的使用和编译器有很大的不同，直接编译出来的版本用到esi寄存器，需要两条保存寄存器指令，两条恢复指令，而我则用了一条寄存器保存和恢复指令，尽量合理调度指令，在同样的执行路径上，指令数也较叫编译器少。
  2. 为了速度快，宁愿写重复的代码，来减少跳转指令的使用。编译出来的版本使用了在case 8路径上需要执行3条条件转移指令，2条jmp跳转指令，而我的版本只用一条条件转移指令和一条jmp指令
  3. "if (r1*r1>=n) r++;" 这样的代码，编译器是使用jbe指令来做的，而我是采用 sbb eax,0 指令来实现的。
  4. 如果没有imul,div指令，我改编后的代码速度提高的就不是10%了。因为imul,div占主要运行时间，所以改用汇编后提速不大。

顺便说一下，我的C版在调用log2时，并没有使用call指令，log2被声明为inline函数，所以汇编提速（10%）全靠减少指令条数和减少跳转指令。

无心人 · 发表于 2009-2-12 11:42:49

呵呵

是否有整数算法不如浮点的结论了？

账号		自动登录	找回密码
密码			欢迎注册

[讨论] 二进制32位整数快速平方根

回复 18# 无心人的帖子

回复 28# 无心人的帖子

[讨论] 二进制32位整数快速平方根

回复 18# 无心人 的帖子

回复 28# 无心人 的帖子

回复 18# 无心人的帖子

回复 28# 无心人的帖子