千分求汇编优化:UInt96x96To192(...)

gxqcn 2007-04-30 01:03:24
// Test.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <stdlib.h>
#if 1
#define _DUMP
#endif

typedef unsigned int UINT32;

/*
L2 L1 L0
x) R2 R1 R0 max
---------------------------- -----------------------------------------------------
L0*R0 FFFFFFFE 00000001
L1*R0 FFFFFFFE FFFFFFFF 00000001
+) L2*R0 FFFFFFFE FFFFFFFF FFFFFFFF 00000001
----------------------------
L0*R1 FFFFFFFF FFFFFFFE 00000000 00000001
L1*R1 00000001 FFFFFFFD FFFFFFFF 00000000 00000001
+) L2*R1 FFFFFFFF FFFFFFFE FFFFFFFF 00000000 00000001
----------------------------
L0*R2 00000001 00000000 FFFFFFFD 00000000 00000000 00000001
L1*R2 00000001 FFFFFFFE FFFFFFFE 00000000 00000000 00000001
+) L2*R2 FFFFFFFF FFFFFFFF FFFFFFFE 00000000 00000000 00000001
---------------------------- -----------------------------------------------------
*/
void UInt96x96To192( UINT32 * pH, UINT32 * pL, const UINT32 * pR )
{
// H[2]:H[1]:H[0] : L[2]:L[1]:L[0] <-- L[2]:L[1]:L[0] * R[2]:R[1]:R[0]
// 即:L、R 分别为 96bits 无符号数,将其相乘后,
// 结果的低 96bits 存在 L 中,高 96bits 存在 H 中。
__asm
{
// 请用 SSE2 或更适合的指令集,以充分利用流水线,加速该段算法的执行,help!
//...
}
}

int main(int argc, char* argv[])
{
// printf("Hello World!\n");

UINT32 u32H[3];
UINT32 u32L[3] = { 0xFFFFFFFF, 0xFFFFFFFE, 0xFFFFFFFD };
const UINT32 u32R[3] = { 0xFFFFFFFD, 0xFFFFFFFE, 0xFFFFFFFF };

#ifdef _DUMP

printf( "%08X %08X %08X * %08X %08X %08X",
u32L[2], u32L[1], u32L[0], u32R[2], u32R[1], u32R[0] );
#endif


UInt96x96To192( u32H, u32L, u32R );


#ifdef _DUMP

printf( "\r\n\t= %08X %08X %08X : %08X %08X %08X\r\n\r\n",
u32H[2], u32H[1], u32H[0], u32L[2], u32L[1], u32L[0] );
#endif


system( "pause" );
return 0;
}

/********************* end of Test.cpp ****************************/



我有一个核心模块:UInt96x96To192() 需要反复调用,所以对效率要求非常高,
希望能通过调用当前最先进的汇编指令,以使其可在尽可能少的 clock 周期内完成。


可惜我的汇编水平还很差,所以求助本版高手,望能圆满解决。
(如达到预期效果,将不惜千分相赠;来着讨论亦均有分)
...全文
3887 点赞 收藏 131
写回复
131 条回复
切换为时间正序
当前发帖距今超过3年,不再开放新的回复
发表回复
thewangj 2011-03-21
这里真是高手如云啊
佩服!学习!
回复
szy41 2008-09-24
mark
回复
yaos 2008-03-31
而且似乎128乘按比例计算已超过96乘了
哈哈
回复
gxqcn 2008-03-25
这个帖子简直可以做为学习 SSE2 汇编等的案例教程:)

现在有个更复杂的擂台赛:UInt128x128To256(...),现在已有非常美妙的参赛作品,
感兴趣者请访问:x86上128位二进制乘法最快速算法征解 - 编程擂台 - 数学研发论坛 (bbs.emath.ac.cn)
回复
无心人 2008-03-03
怎么没人考虑下MMX 寄存器下的运算啊?
慢不慢写个测试下
回复
ollydbg23 2007-12-27
你的这个做的非常好,奇怪,这么久了居然没有人顶一下,
我支持!谢谢!
回复
gxqcn 2007-11-14
(为了与各位网友更好的沟通交流,我新申请了一个收费空间和顶级域名,欢迎大家访问)

中文名称:数学研发网

网址链接:http://www.emath.ac.cn/

Logo图片:

网站简介:学术性数学网站;知识与趣味相交融;提供原创数学工具软件。以等幂和、幻方等数学问题为主要研讨对象,站内有许多原创作品,部分尚为国内外研究的空白。本站集知识性、趣味性于一体;从中您将充分感受到“数学美”的魅力!

关 键 词:数学,研发,等幂和,幻方,数论,组合数学,趣味数学,软件 | mathematics,rd,magic square,software,HugeCalc

备 注:HugeCalc 是一款高精度算法库(同时支持 MBCS + UNICODE 版),适合于大规模科学计算,尤其适用于数论、密码学等领域研究,其核心算法耗费作者十余年的心血。具有占用资源少、效率高、使用便捷、易二次开发、可移植性强、可扩展性好等特点。关键文件 HugeCalc.dll 虽然很小,却提供了公共函数接口 709 个(标准C++接口 473 个;标准C接口 236 个),且其计算速度完全可与大型专业数学工具软件媲美!(在双核上测试,精确计算 40,000,000!,HugeCalc 比最新的 Mathematica V6.01 快了 58%!)
回复
JTZY 2007-08-19
好消息,HugeCalc V7.0.0.0 今天正式发布了!

在潜心研究数月后,终于取得了本质性的突破!
V7.0 版可自动侦测用户 CPU 的型号,并据此自动调整算法及相应参数,使在兼顾老式机器的前提下,
可充分发挥现代及未来 CPU 的功效(如采用 SSE2 指令集、多核并行等)。
在双核上测试,精确计算 40,000,000!,HugeCalc 比最新的 Mathematica V6.01 快了 40%!

大家可下载 HugeCalc 体验:http://www.freewebs.com/maths/software.htm#02


本帖早该结帖了,只是为了等到这一激动的时刻。。。
大家可以在下面这个帖子中继续交流讨论(“千万阶乘擂台”):
http://community.csdn.net/Expert/TopicView.asp?id=5650085

是该结帖了,由于 liangbch(宝宝) 和 DelphiGuy() 之前已各赠奉了 500 分,
所以在结帖时给分将尽量向其他朋友倾斜,还请上述二位勿怪。
回复
vlient 2007-08-19
好久都没看到这么精彩的论道,可惜对汇编优化没有过多的研究,路过,Mark
回复
豆腐 2007-08-02
好帖子,欣赏了!
回复
housisong 2007-07-13
呵呵 SSE4快出来了
回复
JTZY 2007-07-13
今天查阅了 Intel SSE4 Programming Reference,惊见新增了求 MIN/MAX 的指令,尤其是迫切需要的 PMAXUD/PMINUD 指令,虽然仍没有更适合的 unsigned dword integers 直接比较大小的指令,但仍可进一步优化本主题所讨论的函数。

突然,又发现了一条绝妙的指令:PCMGTQ(quadwords 比较大小),与 PADDQ 配合使用也许会有更简洁的方案!

可惜,我没有这么先进的机器。。。
回复
housisong 2007-07-12
更正一下:“把MOVQ改成MOVNTQ就可以了”有误,MMX寄存器才使用MOVNTQ
回复
housisong 2007-07-12
这个帖子好啊:) 好久没有看到这么精彩的讨论帖子了

这里的函数还没有考虑手工内存预读和禁止写缓存;
这两个优化对于大数乘法应该有用(我没有测试过);对于这里的简单测试缓存优化几乎没有
用,因为这里的测试方式CPU默认缓存策略很适合,但在操作大量内存的大数乘法中,应该
考略缓存策略;

具体实现方法:
预读: 在真正使用到源数据之前几百个CPU周期之前(时间由内存和CPU的速度差异
决定)手工读一次其内存数据,也就是将间隔64字节(当前CPU的缓存行大小)内存数据
其中的一个加载一次;
实现方式: a: 手工预读伪代码 (假设源数据在esi寄存器)
for (i;i<BestBufSize;i+=64) //预读数据总大小应该小于CPU的缓存
mov eax,[esi+i]
开始处理esi到esi+BestBufSize的源数据
b:使用专用CPU预读指令prefetch*替换mov eax,[esi+i]
这类指令prefetchnta、prefetcht0、prefetcht1、prefetcht2等作
用是让CPU默认预读、或者指定预读数据到一级缓存或二级缓存...
a的方案比较通用,不需要CPU支持专用的预读指令,而且在我以前的P4核心赛扬
上a方案比b方案还快;
b方案对后来的CPU更友好一些,在我的AMD64X2上b比a快;

禁止写缓冲:当需要写入内存的数据比较多,或不需要立即再次访问写的数据的时候;可
以使用禁止写缓冲优化策略,减少缓冲区污染;
实现方式:使用CPU提供的流式写内存指令MOVNT*
这些指令包括:MOVNTQ、MOVNTI(SSE3)、MOVNTDQ、MOVNTPD,MOVNTPS,
(MOVNTSD、MOVNTSS在SSE4指令集支持)等等
这些指令针对了不同的寄存器类型或者寄存器里的不同数据类型
(sfence指令配合用于刷新写缓存,一般在处理完一块数据以后使用)

完整的缓存优化还包括 把函数内的局部数据布局到64byte的缓存行,寄存器不够用时可以使用
局部堆栈变量(局部变量)来代替(特别是对于常量)
不建议使用各处分布的全局变量,建议可以集中放到一个结构里,从而使他们占用相同的缓存行。。。

(对于内存地址数值给缓存时造成的hash碰撞(CPU的缓存行路数有限),如果发生也要考略避过(发生的情况比较少),...)

大整数运算中乘法是最核心的部分,而大整数乘法复杂度为O(N*log(N)),缓存优化的效果所能起的效果需要楼主测试一下
写缓冲优化是最容易实现的,把MOVQ改成MOVNTQ就可以了:)


回复
JTZY 2007-07-11
本贴之前有近一个月无新贴,但却并没有及时被关闭——请恕我这次的“食言”。。。



可以说通过本贴,我的汇编功力大增,
以前看到公式首先联想到的是 C 代码,现在则直接浮现为汇编代码了。:)

这期间,我把精力主要集中于改进大数算法核心乘法部分,终于取得了令人欣喜的进展,
将最最核心的顶级大数乘法算法效率提高了近一倍左右!(代价是手工编写了大量汇编代码)


从现有的数款 CPU 及不同的操作系统的测试来看,在计算“1千万的阶乘”时均提速了 50% 以上,
如10000000!: Factorial_old(t0/s) Factorial_new(t1/s) (t0 - t1)/t1 * 100%
机台A 68.929053 43.584557 58.15%
机台B 48.787876 31.227357 56.23%
机台C 65.850629 42.810574 53.82%
机台D 41.436746 25.617276 61.75%
其中,
机台A —— Windows XP SP2,Intel Pentium 4 CPU 2.93GHz,512MB DDR - 133MHz
机台B —— Windows XP SP2,AMD Athlon 64 Processor 3200+,1GB DDR - 200MHz
机台C —— Windows Vista 64,Intel Pentium 4 CPU 3.00GHz,1GB DDR - 200MHz
机台D —— Windows XP SP2,AMD Athlon 64 X2 Dual Core Processor 4800+,2GB DDR2 - 800MHz



如果大家感兴趣,请从下列两个站点之一下载(169 KB,本链接仅在正式发布前有效):
http://hugecalc.ik8.com/download/Factorial.rar
http://www.freewebs.com/maths/download/Factorial.rar

欢迎您将测试结果直接发布出来,或填表(我已制作好了 Excel 表)后反馈给我。。。


如果您发现计算“1千万的阶乘”还有更快的程序,请告诉我,必将千分相赠!
回复
JTZY 2007-05-29
我在 Intel PC 上用 VC6.0 + ICC9.1 编译的程序,拿到 AMD CPU 上测试一切正常,
而在家里的 Intel CPU 上却发生不可思议的“离奇现象”,太奇怪了!
回复
JTZY 2007-05-28
从 zyl910 给出的测试数据来看,MUL 所用时间仅占全部的 1/3 不到,
所以如何快速算出这些积的和,如何处理进位反而是关键。

liangbch 上述方法有点借用前面优化的 ALU 算法的味道,
虽然多了 PMULUDQ 次数,但却减少了数据分拆的次数,所以总指令数并不算太多。

我最后给出的方法是一种更好推广的算法,如果有 k 组 XMM 需进位相加,
无需拆分数据,只需 2*(k-1) 次 PADDQ 相加后,及 10 条指令进行修正即可。

另,我用 VC6.0 + ICC8.1 测试,前天早上的“离奇现象” 没有发生,有点莫名其妙哇?。。。
回复
liangbch 2007-05-28
这个一个返璞归真的算法,之所以这么说,是因为:
1。这个程序中没有用到任何复杂的技巧,所有的paddq只处理低64bit,高64bit从不使用,所有的 PMULUDQ指令只计算2个32bit的乘积(共使用了9条PMULUDQ)。
2。使用的指令种类很少,没有用到 pxor,pand,pcmpeqd, PSUBD,PSRLQ,PSLLQ等。总的特点是简单,多次重复。

这个程序的速度并不慢,到目前为止的所有函数中,仅慢于UInt96x96To192_SSE2_GxQ_070524 和 UInt96x96To192_SSE2_GxQ_070523 这两个程序,其中仅比 UInt96x96To192_SSE2_GxQ_070523 慢2.8%,可视为这两个程序速度相当。

从某些方面来看,这个程序编程思想有点反其道而行之的意味,没有充分发挥SSE2的潜力,竟用了9次乘法,但其表现却相当出色。本程序的核心思想来源与GMP.


下面给出代码:
_declspec(naked) void UINT96x96To192SSE2_lbc_0528(UINT32 *dest, const UINT32 *pL, const UINT32 *pR)
{
__asm
{
MOV EDX,[ESP+12] //EDX:PR
MOV EAX,[ESP+8] //EAX:PL
MOV ECX,[ESP+4] //ECX:dest
MOVDQA XMM6,[EAX] //XMM6:PL[0],PL[1],PL[2],0

////////////////////////////////

MOVD XMM3,[EDX] //XMM3:PR[0]
MOVD XMM4,[EDX+4] //XMM4:PR[1]

PMULUDQ XMM3,XMM6 //PR[0]*PL[0]
PMULUDQ XMM4,XMM6 //PR[1]*PL[0]

MOVD [ECX],XMM3
PSRLDQ XMM3,4
////////////////////////////////

MOVD XMM5,[EDX+8] //XMM5=PR[2]

PADDQ XMM4,XMM3

PMULUDQ XMM5,XMM6 //PR[2]*PL[0]

PSHUFD XMM0,XMM4,11111100b //
PSRLDQ XMM4,4 //XMM1:

PADDQ XMM5,XMM4
PSRLDQ XMM6,4
PSHUFD XMM1,XMM5,11111100b //
PSHUFD XMM2,XMM5,11111101b //

///////////////////////////////

MOVD XMM3,[EDX] //XMM3:PR[0]
MOVD XMM4,[EDX+4] //XMM4:PR[1]

PMULUDQ XMM3,XMM6 //XMM3:PR[0]*PL[1]
PMULUDQ XMM4,XMM6 //XMM4:PR[1]*PL[1]

PADDQ XMM3,XMM0

MOVD XMM5,[EDX+8] //XMM5=PR[2]

MOVD [ECX+4],XMM3
PSRLDQ XMM3,4
////////////////////////////////

PADDQ XMM4,XMM1
PADDQ XMM4,XMM3

PMULUDQ XMM5,XMM6 //XMM5:PR[2]*PL[1]

PSHUFD XMM0,XMM4,11111100b
PSRLDQ XMM4,4

////////////////////////////////
PADDQ XMM5,XMM2
PADDQ XMM5,XMM4

PSRLDQ XMM6,4

PSHUFD XMM1,XMM5,11111100b
PSHUFD XMM2,XMM5,11111101b

/////////////////////////////

MOVD XMM3,[EDX] //XMM3:PR[0]
MOVD XMM4,[EDX+4] //XMM4:PR[1]

PMULUDQ XMM3,XMM6 //XMM3:PR[0]*PL[2]
PMULUDQ XMM4,XMM6 //XMM3:PR[1]*PL[2]

PADDQ XMM3,XMM0
MOVD XMM5,[EDX+8] //XMM5=PR[2]

MOVD [ECX+8],XMM3
PSRLDQ XMM3,4
////////////////////////////////

PMULUDQ XMM5,XMM6 //XMM5: PR[2]*PL[2]

PADDQ XMM4,XMM1
PADDQ XMM4,XMM3

MOVD [ECX+12],XMM4
PSRLDQ XMM4,4
////////////////////////////////

PADDQ XMM5,XMM2
PADDQ XMM5,XMM4
MOVQ [ECX+16],XMM5
RET
}
}
回复
JTZY 2007-05-26
在我的机器还发现一个更离奇的现象:

#if 1 // <-- A
__declspec(align(16)) UINT32 u32L[3];
__declspec(align(16)) UINT32 u32R[3];
/*__declspec(align(16))*/ UINT32 u32P[6];
#else // <-- B
__declspec(align(16)) UINT32 u32LBuf[3];
__declspec(align(16)) UINT32 u32RBuf[3];
/*__declspec(align(16))*/ UINT32 u32PBuf[6];

UINT32 *u32L = (UINT32 *)u32LBuf;
UINT32 *u32R = (UINT32 *)u32RBuf;
UINT32 *u32P = (UINT32 *)u32PBuf;
#endif

“#if 1” 对应的是选项“A”;“#if 0” 对应的是选项“B”。测试结果如下:
A B
UInt96x96To192_ALU_v3 4046 4062
UINT96x96To192SSE2_lbc3 3343 3265
UINT96x96To192SSE2e 2843 3296
UInt96x96To192_SSE2_GxQ_org 2515 4234

测试结果表明:
1、切换 A/B 选项对 ALU 系列版本几无影响;
2、切换 A/B 选项对 SSE2 系列影响很大;
3、切换 A/B 选项对 SSE2 系列的不同版本影响程度不一致,甚至将排名反序。
不知该怪现象在其它网友机器上是否也可复制出来?

讨论:
仅仅是传数组名与传指针的区别,不同版本的自差异咋就这么大呢?是哪些指令这么敏感?
回复
JTZY 2007-05-26
下面这个版本是 UInt96x96To192_SSE2_GxQ_070524 版的原型,未作指令混排,更容易理解。
该算法非常精妙,所以代码可以非常精炼;但原理有点诡异,证明有些复杂;
为方便大家阅读,我已补上相应的注释。

to DelphiGuy():
请将这版与 UInt96x96To192_SSE2_GxQ_070524 对比测试一下,着重如下方面:
1、速度孰快?
2、依赖编译选项吗?
3、buffer 对它影响大吗?(在我的机器上即便开 “UINT32 u32PBuf[40960];” 也无影响。)



/*
V5:V4:V3:V2:V1:V0
-----------------------------
L2*R0:L0*R0 <-- xmm0
L2*R1:L0*R1 <-- xmm1
L2*R2:L0*R2 <-- xmm2
L1*R2:L1*R0 <-- xmm3
L1*R1 <-- xmm4
*/
_declspec(naked) /* 该版算法可能有点诡异哦:) */
void UInt96x96To192_SSE2_GxQ_org( UINT32 *p, const UINT32 *pL, const UINT32 *pR )
{
__asm
{
mov edx, dword ptr[esp + 0Ch] ; pR
mov eax, dword ptr[esp + 08h] ; pL
mov ecx, dword ptr[esp + 04h] ; p

;; when data is 128-bit unaligned, please use 'movdqu' or 'lddqu'(SSE3).
movdqa xmm7, [edx] ; load pR
movdqa xmm6, [eax] ; load pL

pshufd xmm0, xmm7, 00000000b ; xmm0= R0:R0:R0:R0
pshufd xmm1, xmm7, 01010101b ; xmm1= R1:R1:R1:R1
pshufd xmm2, xmm7, 10101010b ; xmm2= R2:R2:R2:R2
pshufd xmm3, xmm6, 01010101b ; xmm3= L1:L1:L1:L1
movq xmm4, xmm1 ; xmm4= 00:00:R1:R1

pmuludq xmm0, xmm6 ; xmm0= L2*R0:L0*R0 *
pmuludq xmm4, xmm3 ; xmm4= 00:00:L1*R1 *
pmuludq xmm1, xmm6 ; xmm1= L2*R1:L0*R1 *
pmuludq xmm2, xmm6 ; xmm2= L2*R2:L0*R2 *
pmuludq xmm3, xmm7 ; xmm3= L1*R2:L1*R0 *


movd dword ptr[ecx], xmm0 ; store V0
psrldq xmm0, 4 ; xmm0 >>= 32

movdqa xmm7, xmm1 ; xmm7 = xmm1 <-- [V4:V3][V2:V1]
paddq xmm7, xmm0 ; + xmm0
paddq xmm7, xmm3 ; + xmm3

psrldq xmm0, 4 ; xmm0 >>= 32 <-- [V5:V4][V3:V2]
psrldq xmm1, 4 ; xmm1 >>= 32
psrldq xmm3, 4 ; xmm3 >>= 32

paddq xmm0, xmm1 ; xmm0 += xmm1 + xmm2 + xmm3 + xmm4
paddq xmm0, xmm2 ;
paddq xmm0, xmm3 ;
paddq xmm0, xmm4 ;

pslldq xmm2, 4 ; xmm2 <<= 32
pslldq xmm4, 4 ; xmm4 <<= 32
paddq xmm2, xmm7 ; xmm2 += xmm4 + xmm7 <-- [V4:V3][V2:V1]
paddq xmm2, xmm4 ; now, xmm2=[**:**][V2:V1]! V2:V1 OK!

movd dword ptr[ecx + 04h], xmm2 ; store V1


psrldq xmm2, 4 ; xmm2 >>= 32
movdqa xmm1, xmm2 ; xmm1 = xmm2 - xmm0
psubd xmm1, xmm0 ;
psllq xmm1, 32 ;
psrlq xmm1, 32 ; xmm1 ^= 00:FF:00:FF
paddq xmm0, xmm1 ; now, xmm0[63-0] bits OK!

psubusb xmm2, xmm0 ; if CF=0: 0 <= xmm0[63-32]-xmm2[63-32] <= 4
psrlq xmm2, 63 ; if CF=1: xmm0[63]==0 && xmm2[63]==1
pslldq xmm2, 8 ;
paddq xmm0, xmm2 ; xmm0[127-64] += ( CF << 64 )


movdqu [ecx + 08h], xmm0 ; store [V5:V4:V3:V2]
ret
}
}
回复
加载更多回复
相关推荐
发帖
汇编语言
创建于2007-08-27

2.1w+

社区成员

汇编语言(Assembly Language)是任何一种用于电子计算机、微处理器、微控制器或其他可编程器件的低级语言,亦称为符号语言。
申请成为版主
帖子事件
创建了帖子
2007-04-30 01:03
社区公告
暂无公告