大整数加法----（10进制，速度快）

只是呼吸 发表于 2016-12-11 00:24:16

本帖最后由只是呼吸于 2016-12-11 02:26 编辑

在大整数的运算程序中，无符号加法应该算最容易写的了。但要把加法的运算程序尽可能运行快一点，仍然要动一番脑筋。
   这个加法程序是我在本论坛学习的过程中，参照了管理员liangbch的两篇文章写成的。

            http://bbs.emath.ac.cn/thread-216-11-1.html 在这篇文章的第102楼，学习了汇编语言的接口技术。

            http://bbs.emath.ac.cn/thread-521-3-1.html    在这篇文章，学习了把单精度除法转化为乘法的技术。

            下面是我写的大整数加法程序：
/*
*声明：本函数的接口技术、单精度乘法代替单精度除法技术，学习自管理员liangbch的两篇帖子：http://bbs.emath.ac.cn/thread-216-11-1.html第102楼和 http://bbs.emath.ac.cn/thread-521-3-1.html
*       liangbch依法享有著作权和解释权
*\
*dst: 目的数组，存放和。
*src1:第一个加数，（操作数）
*src2:第二个加数，（操作数）
*size1; 第一个加数的长度。
*size2: 第二个加数的长度。
*       所有数组均采用10^8进制，数组较大的下标对应10进制数的高位，数组较小的下标对应10进制数的低位。
*/

_declspec(naked)////这个程序能处理任意长度的两个无符号大整数之和，速度快。
void big_add_ALU(unsigned long*dst,unsigned long *src1,unsigned long *src2,int size1,int size2)
{
#undefBASE10000_MASK
#define BASE10000_MASK BASE
   _asm
   {
push esi
   push edi
   push ebx
   push ebp

   mov       edi, dword ptr
mov       ebp, dword ptr
xor       eax,eax
mov       ebx, BASE10000_MASK
movd    mm4,eax
   mov       esi, dword ptr
mov       eax, 0xabcc7712
movd    mm7,ebx
movd    mm5,eax

   xor       ebx,ebx
mov       ecx, dword ptr
   mov       eax, dword ptr
cmp       eax, ecx
jl       fffd

jmp       jjj
fffd:
mov       ecx,eax //< 则跳转
jjj:

loop02:

movd mm1, dword ptr
movd mm2, dword ptr

paddq mm1, mm2
paddq mm1, mm4

movq mm2, mm1 // mm2, 总和
pmuludq mm1,mm5
psrlq mm1,58
movq mm4,mm1// mm4进位
pmuludq mm1,mm7
psubq mm2,mm1

movd dword ptr,mm2
lea edi, dword ptr
   lea ebp, dword ptr
   lea esi, dword ptr

dec ecx
   jnz loop02

//--------以上是公共部分之和------------------//

mov       ecx, dword ptr
   mov       eax, dword ptr
cmp       eax,ecx
je       exit1
jl       fffdd
sub       eax,ecx
mov       ecx,eax

loop03:
movd    eax,mm4
   cmp       eax,ebx
je       exit2

movd mm1, dword ptr
paddq mm1,mm4
movq mm2,mm1 // mm2,总和

pmuludq mm1,mm5
psrlq mm1,58
movq mm4,mm1//    mm4 进位

pmuludq mm1,mm7
psubq mm2,mm1

movd dword ptr,mm2

lea esi, dword ptr
lea ebp, dword ptr
dec ecx
jnz loop03//////
exit2:
cmp       ecx,ebx
je       exit1

loop04:
mov    eax,dword ptr
mov    dword ptr,eax

lea esi, dword ptr
lea ebp, dword ptr
dec ecx
jnz loop04/////

   jmp       exit1
///-----------------以上是第一个加数大于第二个加数的求和-------------//
fffdd:

   sub       ecx, eax
loop05:
movd    eax,mm4
   cmp       eax,ebx
je       exit3

movd mm1,dword ptr
paddq mm1,mm4

movq mm2,mm1 //mm2,总和
pmuludq mm1,mm5
psrlq mm1,58
movq mm4,mm1 // mm4 进位
pmuludq mm1,mm7
psubq mm2,mm1


movd dword ptr,mm2

lea edi, dword ptr
lea ebp, dword ptr
dec ecx
jnz loop05//////

exit3:
cmp       ecx,ebx
je       exit1

loop06:
mov    eax,dword ptr
mov    dword ptr,eax

lea edi, dword ptr
lea ebp, dword ptr
dec ecx
jnz loop06/////

//--------------------------------------------//
exit1:

movd dword ptr ,mm4

pop ebp
   pop ebx
   pop edi
   pop esi

emms
   ret
   }
}/////*///

作为对照，我写了一个c语言的大整数加法函数
/*
*这是c语言写的大整数加法，作为对照。
*
*为了简单一点：
*
*1）要求第一个加数src1的实际长度大于等于第二个加数src2的实际长度。
*2）第二个加数的前面补零，与第一个加数的实际长度对齐。（实际处理中只需第二个加数src2的最大空间大于src1的实际长度即可）
*3）参数size1传递第一个加数的实际长度。
*/

void big_add_c(unsigned long *dst,unsigned long *src1,unsigned long *src2,int size1)
{
int i,temp=0;
for(i=0;i<size1;i++)
{
dst=src1+src2+temp;
      temp=0;
   if(dst>=BASE)
{
dst-=BASE;
temp=1;
}
}
dst=temp;
}

只是呼吸 发表于 2016-12-11 01:05:34

本帖最后由只是呼吸于 2016-12-11 03:10 编辑

上面的两个程序，c语言的作为参考。输入两个大整数进行对比。用vc++2008在debug模式下进行实验。

      输入数字规模                                              调用函数次数                                     c 语言运行时间                         汇编语言运行时间
         两个14700位的随机数相加                                     1000 000 次                                        19.157秒                                  5.944秒
         两个14700位的9（全是9）相加                               1000 000 次                                        15.616秒                                  5.881秒
         30000位和 14700位的随机数相加                         1000 000 次                                        30.888秒                                  7.800秒
         30000位和 14700位的9（全是9）相加                   1000 000 次                                        31.800秒                                  12.324秒

上面的数据只记录函数的运行时间，不计输出时间，这是我在电脑上的实测结果，如果机器的性能不同，数字会有变化。

输入两个全是9的大整数相加，主要是考验程序的连续进位。函数big_add_ALU()采用的是两次乘法一次减法来代替判断语句，综合性能较好，在汇编语言中如果采用判断语句来进位，只有在连续进位或连续不进位的情况下占优。例如全部是9相加（连续进位），或者全部是3相加，（连续不进位）这两种情况下的实际运行时间比big_add_ALU()要好得多。而对随机数来说，采用判断语句来进位的函数就比不过big_add_ALU()了。

只是呼吸 发表于 2016-12-11 01:16:05

本帖最后由只是呼吸于 2016-12-11 02:32 编辑

这里给出这个函数的调用方法,用两个很小的数字来说明。
这里我用的是vc++2008，（vc++6.0不能编译）。其它的编译器不知能不能用。

#include <stdio.h>

#defineBASE 100000000 //10^8进制。(共8个零)，32位系统中不能超过10^9。

//char a={0},b={0};
unsigned longc={0},e={0},result={0};
/*
//---------------------------------------
在这里粘一楼的两个函数

//----------------------------------------
*/
int main()
{
inti,m1,m2;

// 大数：c[]138850000000090000099934222009345001230023908770390239
// 大数：e[]                      38462690000660060508019900031415

   c=70390239;// 数组较小的下标对应大整数的低位，采用10^8进制，（可以修改为10^9进制，这样效率更高一些）。
c=239087;
c=34500123;
c=34222009;
c=999;
c=9;
c=138850;

e=31415;
e=5080199;
e=66006;
e=38462690;

m1=7;// c[]的长度
m2=4;// e[]的长度

//for(i=0;i<1000000;i++) // 调用函数100万次，用这种方法估计时间。虽不很准，但可以说明问题。
//big_add_c(result,c,e,m1);// c语言
big_add_ALU(result,c,e,m1,m2);

m1=m1+m2+3;
while((result==0)&&(m1>0))//寻找不为零的第一个数
{
m1--;
}

printf("%d",result);//输出语句。
for(i=m1-1;i>=0;i--)//输出语句。把相加的结果输出。
printf("%08d",result);//输出语句。

printf("\n");

// 计算的结果是： 138850000000090000099972684699345661290531928670421654

system( "pause" );

return 0;

}

liangbch 发表于 2016-12-11 10:40:39

关于加法运算的正常化

1. RAD表示进制，若 a < RAD， b< RAD，c<=1, (a是被加数，b是加数，c是进位）
则, s0=a+b+c, 这时s0<RAD*2, 将s0拆为本位s和进位c的方法(即求s和c,使得s0=c*RAD+s)称为规约，英文直译为“正常化”。
在这里，规约的过程是不需要使用乘法的，一个显然的实现是
if (s0>=RAD)
{s=s0-RAD; c=1; }
else
{s=s0; c=0;}

上面的实现用到了分支指令，由于两个分支的执行机会均等，CPU分支预测机构无论预测那个分支，预测失败的概率总是50%，所以效率不高。一个好的实现方法应该完全不使用分支指令
s= s0- RAD;
tmp = -(s>>31);
c= tmp+1;
s += (tmp & RAD);

下面给出使用MMX指令的规约的代码
若 MM_RES= a+b+c, MM_RAD=R, MM_RAD_S1= R-1, 当执行完下面的指令序列后，MM_RES 为本位，MM_CARRY为进位。
说明，
1. 这里的进位为负数，所以在求 RES=RES+c时，使用的是PSUBD指令而不是PADDD指令，即PSUBD MM_RES,MM_CARRY
2. 为了便于代码阅读，将寄存器起了一个新的名字。
在C语言中，使用下面的定义
#define MM_RAD MM7
#define MM_RAD_S1 MM6
#define MM_CARRYMM5
#define MM_RES MM1
#define MM_TMP MM2

在汇编语言中，应该这样定义
MM_RAD EQU MM7
MM_RAD_S1 EQU MM6
MM_CARRYEQU MM5
MM_RES EQU MM1
MM_TMP EQU MM2

MOVQ MM_TMP,MM_RES
PCMPGTD MM_TMP,MM_RAD_S1 ; if res>=RAD, tmp=-1, else tmp=0
MOVQ MM_CARRY,MM_TMP ; if res>=RAD, carry= -1, else carry=0
PAND MM_TMP,MM_RAD ; if res>=RAD, tmp= RAD, else tmp=0
PSUBD MM_RES,MM_TMP ; res -= RAD or res-=0

liangbch 发表于 2016-12-11 10:42:05

我写了8个版本的基于10^9的大数加法，包括1个c语言版本，和7个汇编语言版本
mpn_D9_add_n_C       : 消除分支的c语言版本
mpn_D9_add_n_ALU1：：消除分支的ALU版本
mpn_D9_add_n_ALU2       ：消除分支的ALU另一个版本，规约方法和ALU1不同
mpn_D9_add_n_CMOV1 ：使用CMOVCC指令消除分支
mpn_D9_add_n_CMOV2 ；另一种使用CMOVCC指令的版本
mpn_D9_add_n_MMX       ：使用MMX指令的版本，一次计算2个DWORD
mpn_D9_add_n_SSE2       ：使用SSE2指令的版本，一次计算4个DWORD
mpn_D9_add_n_SSE4       ；在SSE2指令的版本的基础上，使用了一条SSE4指令 PTEST 来加速。

c语言语言版的源代码
DWORD mpn_D9_add_n_C(DWORD *rp,const DWORD *s1p, const DWORD *s2p, DWORD len)
{
const DWORD *p1End=s1p+len-1;

DWORD carry=0;
DWORD res;

int tmp;
while (s1p <= p1End)
{
res = *s1p + *s2p +carry;
res -= DEC9_RADIX;
tmp = -(res>>31);
carry= tmp+1;
res += (tmp & DEC9_RADIX);
*rp=res;
s1p++; s2p++; rp++;
}

return carry;
}

liangbch 发表于 2016-12-11 10:43:19

性能测试报告：

blk_len mpn_D9_add_n_C mpn_D9_add_n_ALU1 mpn_D9_add_n_ALU2 mpn_D9_add_n_CMOV1 mpn_D9_add_n_CMOV2 mpn_D9_add_n_MMX mpn_D9_add_n_SSE2 mpn_D9_add_n_SSE4
8192 4.602 4.670 3.582 4.642 3.850 5.337 1.445 1.285
16384 4.609 4.640 3.575 4.582 3.831 5.331 1.438 1.271
32768 4.590 4.643 3.549 5.172 3.836 5.333 1.446 1.298
65536 4.600 4.660 3.550 4.586 3.835 5.333 1.449 1.293

说明：
大整数采用10^9进制，每个DWORD（32bit整数）可表示9位10进制数，上表列出当被加数和加数的长度为8K,16K，32K和64K个DWORD时，平均计算每个DWORD所耗费的时间，单位为时钟周期，数字越小表示性能越好。

无心人 发表于 2023-11-23 16:28:50

如果单纯是用mmx寄存器，而不是因为机器太老，用不了sse指令的话，建议mm寄存器换成xmm寄存器，
理论上，现在的CPU，并不存在物理的mm寄存器的，都是SIMD寄存器堆，所以，mm和xmm寄存器很可能物理上是一个东西
两者在早期CPU上确实存在速度差异，现在CPU已经没区别了

lihpb00 发表于 2023-11-24 21:22:26

看不懂

页: [1]

数学研发论坛's Archiver

大整数加法----（10进制，速度快）