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

123456789

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

                http://bbs.emath.ac.cn/thread-216-11-1.html     在这篇文章的第102楼，学习了汇编语言的接口技术。
               
                http://bbs.emath.ac.cn/thread-521-3-1.html       在这篇文章，学习了把单精度除法转化为乘法的技术。

                下面是我写的大整数加法程序：

1. /* *  声明：本函数的接口技术、单精度乘法代替单精度除法技术，学习自管理员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){#undef  BASE10000_MASK#define BASE10000_MASK BASE        _asm        {                push                esi        push                edi        push                ebx        push                ebp                mov         edi, dword ptr [esp+0Ch+16]                mov         ebp, dword ptr [esp+04h+16]                xor         eax,eax                mov         ebx, BASE10000_MASK                movd        mm4,eax        mov         esi, dword ptr [esp+08h+16]                mov         eax, 0xabcc7712                movd        mm7,ebx                movd        mm5,eax                        xor         ebx,ebx                mov         ecx, dword ptr [esp+14h+16]        mov         eax, dword ptr [esp+10h+16]                cmp         eax, ecx                jl          fffd                                jmp         jjj        fffd:                  mov         ecx,eax   //  < 则跳转jjj:         loop02:                 movd    mm1, dword ptr [esi]                movd    mm2, dword ptr [edi]                                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[ebp],mm2                lea                        edi, dword ptr [edi+4]        lea                        ebp, dword ptr [ebp+4]        lea                        esi, dword ptr [esi+4]                        dec                        ecx        jnz                        loop02//--------以上是公共部分之和------------------//                mov         ecx, dword ptr [esp+14h+16]        mov         eax, dword ptr [esp+10h+16]                cmp         eax,ecx                je          exit1                jl          fffdd                sub         eax,ecx                mov         ecx,eaxloop03:                movd        eax,mm4        cmp         eax,ebx                je          exit2                movd                mm1, dword ptr [esi]                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[ebp],mm2                lea                        esi, dword ptr [esi+4]                lea                        ebp, dword ptr [ebp+4]                dec                        ecx                jnz                        loop03//////exit2:                cmp         ecx,ebx                je          exit1                loop04:                 mov        eax,dword ptr [esi]                mov        dword ptr[ebp],eax                                lea                        esi, dword ptr [esi+4]                lea                        ebp, dword ptr [ebp+4]                dec                        ecx                jnz                        loop04/////                                jmp         exit1///-----------------以上是第一个加数大于第二个加数的求和-------------//fffdd:                sub         ecx, eaxloop05:                movd        eax,mm4        cmp         eax,ebx                je          exit3                movd                mm1,dword ptr [edi]                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[ebp],mm2                        lea                        edi, dword ptr [edi+4]                lea                        ebp, dword ptr [ebp+4]                dec                        ecx                jnz                        loop05//////exit3:                cmp         ecx,ebx                je          exit1                loop06:                 mov        eax,dword ptr [edi]                mov        dword ptr[ebp],eax                                lea                        edi, dword ptr [edi+4]                lea                        ebp, dword ptr [ebp+4]                dec                        ecx                jnz                        loop06/////        //--------------------------------------------//exit1:                        movd                dword ptr [ebp],mm4              pop                        ebp        pop                        ebx        pop                        edi        pop                        esi                                emms        ret        }}/////*///

复制代码

作为对照，我写了一个c语言的大整数加法函数

1. /* *  这是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[i]=src1[i]+src2[i]+temp;              temp=0;            if(dst[i]>=BASE)                {                        dst[i]-=BASE;                        temp=1;                }        }        dst[i]=temp;}

复制代码

uzmedafacof

上面的两个程序，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()了。

obibohazek

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

1. #include <stdio.h>#define  BASE 100000000 //10^8进制。(共8个零)，32位系统中不能超过10^9。//char a[210000]={0},b[200000]={0};unsigned long  c[30000]={0},e[30000]={0},result[2000000]={0};/*//---------------------------------------在这里粘一楼的两个函数//----------------------------------------*/int main(){        int  i,m1,m2;            // 大数：c[]  138850000000090000099934222009345001230023908770390239// 大数：e[]                        38462690000660060508019900031415                         c[0]=70390239;  // 数组较小的下标对应大整数的低位，采用10^8进制，（可以修改为10^9进制，这样效率更高一些）。        c[1]=239087;        c[2]=34500123;        c[3]=34222009;        c[4]=999;        c[5]=9;        c[6]=138850;        e[0]=31415;        e[1]=5080199;        e[2]=66006;        e[3]=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[m1]==0)&&(m1>0))//寻找不为零的第一个数        {                m1--;        }                printf("%d",result[m1]);//输出语句。        for(i=m1-1;i>=0;i--)//输出语句。把相加的结果输出。                printf("%08d",result[i]);//输出语句。        printf("\n");  // 计算的结果是： 138850000000090000099972684699345661290531928670421654          system( "pause" );        return 0;}

复制代码

AzertVab

关于加法运算的正常化

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)称为规约，英文直译为“正常化”。
   在这里，规约的过程是不需要使用乘法的，一个显然的实现是

1.    if (s0>=RAD)   {  s=s0-RAD; c=1; }   else   {  s=s0; c=0;  }

复制代码

上面的实现用到了分支指令，由于两个分支的执行机会均等，CPU分支预测机构无论预测那个分支，预测失败的概率总是50%，所以效率不高。一个好的实现方法应该完全不使用分支指令

1.   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语言中，使用下面的定义

1.   #define   MM_RAD    MM7  #define   MM_RAD_S1 MM6  #define   MM_CARRY  MM5  #define   MM_RES    MM1  #define   MM_TMP    MM2

复制代码

在汇编语言中，应该这样定义

1.    MM_RAD    EQU   MM7   MM_RAD_S1 EQU   MM6   MM_CARRY  EQU   MM5   MM_RES    EQU   MM1   MM_TMP    EQU   MM2

复制代码

1. MOVQ    MM_TMP,MM_RESPCMPGTD MM_TMP,MM_RAD_S1        ; if res>=RAD, tmp=-1, else tmp[i]=0MOVQ    MM_CARRY,MM_TMP                ; if res>=RAD, carry[i]= -1, else carry[i]=0PAND    MM_TMP,MM_RAD                ; if res[i]>=RAD, tmp[i]= RAD, else tmp[i]=0PSUBD   MM_RES,MM_TMP                ; res -= RAD or res-=0

复制代码

uucepekoqijew

我写了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语言语言版的源代码

1. 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;}

复制代码

Dennislem

性能测试报告：

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所耗费的时间，单位为时钟周期，数字越小表示性能越好。

anaboewako

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

eqajoefiwae

看不懂