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

只是呼吸

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

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

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

1. /*
   
2. *  声明：本函数的接口技术、单精度乘法代替单精度除法技术，学习自管理员liangbch的两篇帖子：http://bbs.emath.ac.cn/thread-216-11-1.html第102楼 和 http://bbs.emath.ac.cn/thread-521-3-1.html
   
3. *         liangbch依法享有著作权和解释权
   
4. *\
   
5. *  dst:   目的数组，存放和。
   
6. *  src1:  第一个加数，（操作数）
   
7. *  src2:  第二个加数，（操作数）
   
8. *  size1; 第一个加数的长度。
   
9. *  size2: 第二个加数的长度。
   
10. *         所有数组均采用10^8进制，数组较大的下标对应10进制数的高位，数组较小的下标对应10进制数的低位。
    
11. */
    
12. 
    
13. 
    
14. _declspec(naked)////这个程序能处理任意长度的两个无符号大整数之和，速度快。
    
15. void big_add_ALU(unsigned long  *dst,unsigned long *src1,unsigned long *src2,int size1,int size2)
    
16. {
    
17. #undef  BASE10000_MASK
    
18. #define BASE10000_MASK BASE
    
19.         _asm
    
20.         {
    
21.                 push                esi
    
22.         push                edi
    
23.         push                ebx
    
24.         push                ebp
    
25.         
    
26.         mov         edi, dword ptr [esp+0Ch+16]
    
27.                 mov         ebp, dword ptr [esp+04h+16]
    
28.                 xor         eax,eax
    
29.                 mov         ebx, BASE10000_MASK
    
30.                 movd        mm4,eax
    
31.         mov         esi, dword ptr [esp+08h+16]
    
32.                 mov         eax, 0xabcc7712
    
33.                 movd        mm7,ebx
    
34.                 movd        mm5,eax
    
35.                
    
36.         xor         ebx,ebx
    
37.                 mov         ecx, dword ptr [esp+14h+16]
    
38.         mov         eax, dword ptr [esp+10h+16]
    
39.                 cmp         eax, ecx
    
40.                 jl          fffd
    
41.                
    
42.                 jmp         jjj        
    
43. fffd:  
    
44.                 mov         ecx,eax   //  < 则跳转
    
45. jjj:
    
46.        
    
47. loop02:
    
48. 
    
49.                 movd    mm1, dword ptr [esi]
    
50.                 movd    mm2, dword ptr [edi]
    
51.                
    
52.                 paddq   mm1, mm2
    
53.                 paddq   mm1, mm4
    
54.                
    
55.                 movq                mm2, mm1 // mm2, 总和
    
56.                 pmuludq                mm1,mm5
    
57.                 psrlq                mm1,58
    
58.                 movq                mm4,mm1//   mm4  进位
    
59.                 pmuludq                mm1,mm7
    
60.                 psubq                mm2,mm1
    
61.             
    
62.                 movd                dword ptr[ebp],mm2
    
63.                 lea                        edi, dword ptr [edi+4]
    
64.         lea                        ebp, dword ptr [ebp+4]
    
65.         lea                        esi, dword ptr [esi+4]
    
66.         
    
67.                 dec                        ecx
    
68.         jnz                        loop02
    
69. 
    
70. //--------以上是公共部分之和------------------//
    
71. 
    
72.                 mov         ecx, dword ptr [esp+14h+16]
    
73.         mov         eax, dword ptr [esp+10h+16]
    
74.                 cmp         eax,ecx
    
75.                 je          exit1
    
76.                 jl          fffdd
    
77.                 sub         eax,ecx
    
78.                 mov         ecx,eax
    
79. 
    
80. loop03:
    
81.                 movd        eax,mm4
    
82.         cmp         eax,ebx
    
83.                 je          exit2
    
84. 
    
85.                 movd                mm1, dword ptr [esi]
    
86.                 paddq                mm1,mm4
    
87.                 movq                mm2,mm1 // mm2,总和
    
88. 
    
89.                 pmuludq                mm1,mm5
    
90.                 psrlq                mm1,58
    
91.                 movq                mm4,mm1//          mm4 进位
    
92.                
    
93.                 pmuludq                mm1,mm7
    
94.                 psubq                mm2,mm1
    
95. 
    
96.                 movd                dword ptr[ebp],mm2
    
97. 
    
98.                 lea                        esi, dword ptr [esi+4]
    
99.                 lea                        ebp, dword ptr [ebp+4]
    
100.                 dec                        ecx
     
101.                 jnz                        loop03//////
     
102. exit2:
     
103.                 cmp         ecx,ebx
     
104.                 je          exit1
     
105.                
     
106. loop04:
     
107.                 mov        eax,dword ptr [esi]
     
108.                 mov        dword ptr[ebp],eax
     
109.                
     
110.                 lea                        esi, dword ptr [esi+4]
     
111.                 lea                        ebp, dword ptr [ebp+4]
     
112.                 dec                        ecx
     
113.                 jnz                        loop04/////        
     
114.                
     
115. 
     
116.         jmp         exit1
     
117. ///-----------------以上是第一个加数大于第二个加数的求和-------------//
     
118. fffdd:       
     
119. 
     
120.         sub         ecx, eax
     
121. loop05:
     
122.                 movd        eax,mm4
     
123.         cmp         eax,ebx
     
124.                 je          exit3
     
125. 
     
126.                 movd                mm1,dword ptr [edi]
     
127.                 paddq                mm1,mm4
     
128.                
     
129.                 movq                mm2,mm1 //mm2,总和
     
130.                 pmuludq                mm1,mm5   
     
131.                 psrlq                mm1,58
     
132.                 movq                mm4,mm1 //   mm4 进位
     
133.                 pmuludq                mm1,mm7
     
134.                 psubq                mm2,mm1
     
135.                
     
136. 
     
137.             
     
138.                 movd                dword ptr[ebp],mm2       
     
139. 
     
140.                 lea                        edi, dword ptr [edi+4]
     
141.                 lea                        ebp, dword ptr [ebp+4]
     
142.                 dec                        ecx
     
143.                 jnz                        loop05//////
     
144. 
     
145. exit3:
     
146.                 cmp         ecx,ebx
     
147.                 je          exit1
     
148.                
     
149. loop06:
     
150.                 mov        eax,dword ptr [edi]
     
151.                 mov        dword ptr[ebp],eax
     
152.                
     
153.                 lea                        edi, dword ptr [edi+4]
     
154.                 lea                        ebp, dword ptr [ebp+4]
     
155.                 dec                        ecx
     
156.                 jnz                        loop06/////        
     
157. 
     
158. //--------------------------------------------//
     
159. exit1:
     
160.                
     
161.         movd                dword ptr [ebp],mm4
     
162.       
     
163.         pop                        ebp
     
164.         pop                        ebx
     
165.         pop                        edi
     
166.         pop                        esi
     
167.                
     
168. 
     
169. 
     
170.                 emms
     
171.         ret
     
172.         }
     
173. }/////*///
     

复制代码

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

1. /*
   
2. *  这是c语言写的大整数加法，作为对照。
   
3. *   
   
4. *  为了简单一点：
   
5. *
   
6. *  1）要求第一个加数src1的实际长度大于等于第二个加数src2的实际长度。
   
7. *  2）第二个加数的前面补零，与第一个加数的实际长度对齐。（实际处理中只需第二个加数src2的最大空间大于src1的实际长度即可）
   
8. *  3）参数size1传递第一个加数的实际长度。
   
9. */
   
10. 
    
11. void big_add_c(unsigned long *dst,unsigned long *src1,unsigned long *src2,int size1)
    
12. {
    
13.         int i,temp=0;
    
14.         for(i=0;i<size1;i++)
    
15.         {
    
16.                 dst[i]=src1[i]+src2[i]+temp;
    
17.               temp=0;
    
18.             if(dst[i]>=BASE)
    
19.                 {
    
20.                         dst[i]-=BASE;
    
21.                         temp=1;
    
22.                 }
    
23.         }
    
24.         dst[i]=temp;
    
25. }
    

复制代码

只是呼吸

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

只是呼吸

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

1. #include <stdio.h>
   
2. 
   
3. #define  BASE 100000000 //10^8进制。(共8个零)，32位系统中不能超过10^9。
   
4. 
   
5. 
   
6. //char a[210000]={0},b[200000]={0};
   
7. unsigned long  c[30000]={0},e[30000]={0},result[2000000]={0};
   
8. /*
   
9. //---------------------------------------
   
10. 在这里粘一楼的两个函数
    
11. 
    
12. //----------------------------------------
    
13. */
    
14. int main()
    
15. {
    
16.         int  i,m1,m2;
    
17.        
    
18.     // 大数：c[]  138850000000090000099934222009345001230023908770390239
    
19. // 大数：e[]                        38462690000660060508019900031415
    
20.                  
    
21.         c[0]=70390239;  // 数组较小的下标对应大整数的低位，采用10^8进制，（可以修改为10^9进制，这样效率更高一些）。
    
22.         c[1]=239087;
    
23.         c[2]=34500123;
    
24.         c[3]=34222009;
    
25.         c[4]=999;
    
26.         c[5]=9;
    
27.         c[6]=138850;
    
28. 
    
29.         e[0]=31415;
    
30.         e[1]=5080199;
    
31.         e[2]=66006;
    
32.         e[3]=38462690;
    
33. 
    
34.         m1=7;// c[]的长度
    
35.         m2=4;// e[]的长度
    
36. 
    
37.        
    
38.         //for(i=0;i<1000000;i++) // 调用函数100万次，用这种方法估计时间。虽不很准，但可以说明问题。       
    
39.                 //big_add_c(result,c,e,m1);// c  语言
    
40.                 big_add_ALU(result,c,e,m1,m2);
    
41. 
    
42.         m1=m1+m2+3;
    
43.     while((result[m1]==0)&&(m1>0))//寻找不为零的第一个数
    
44.         {
    
45.                 m1--;
    
46.         }
    
47. 
    
48.        
    
49. 
    
50.         printf("%d",result[m1]);//输出语句。
    
51.         for(i=m1-1;i>=0;i--)//输出语句。把相加的结果输出。
    
52.                 printf("%08d",result[i]);//输出语句。
    
53. 
    
54.         printf("\n");
    
55. 
    
56. 
    
57.   // 计算的结果是： 138850000000090000099972684699345661290531928670421654
    
58. 
    
59.   
    
60. 
    
61.         system( "pause" );
    
62. 
    
63.         return 0;
    
64. 
    
65. }
    

复制代码

liangbch

关于加法运算的正常化

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)
   
2.    {  s=s0-RAD; c=1; }
   
3.    else
   
4.    {  s=s0; c=0;  }
   

复制代码

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

1.   s= s0- RAD;
   
2.   tmp = -(s>>31);
   
3.   c= tmp+1;
   
4.   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
   
2.   #define   MM_RAD_S1 MM6
   
3.   #define   MM_CARRY  MM5
   
4.   #define   MM_RES    MM1
   
5.   #define   MM_TMP    MM2
   

复制代码

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

1.    MM_RAD    EQU   MM7
   
2.    MM_RAD_S1 EQU   MM6
   
3.    MM_CARRY  EQU   MM5
   
4.    MM_RES    EQU   MM1
   
5.    MM_TMP    EQU   MM2
   

复制代码

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

复制代码

liangbch

我写了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)
   
2. {
   
3.         const DWORD *p1End=s1p+len-1;
   
4. 
   
5.         DWORD carry=0;
   
6.         DWORD res;
   
7. 
   
8.         int tmp;
   
9.         while (s1p <= p1End)
   
10.         {
    
11.                 res = *s1p + *s2p +carry;
    
12.                 res -= DEC9_RADIX;
    
13.                 tmp = -(res>>31);
    
14.                 carry= tmp+1;
    
15.                 res += (tmp & DEC9_RADIX);
    
16.                 *rp=res;
    
17.                 s1p++;        s2p++;         rp++;
    
18.         }
    
19. 
    
20.         return carry;
    
21. }
    
22. 
    

复制代码

liangbch

性能测试报告：

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

无心人

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

lihpb00

看不懂