不用64位乘除法，求 uint64_t 对 10^9 的余数

liangbch

不错。使用模算术的方法，代码简洁。
我的算法仍然是普通的除法运算，很难化简到这么简单。

gxqcn

经过变换后，问题的核心转化为计算一个 \(56\) 比特的整数对 \(2*5^9\) 的余数（最后再附加一次与 \(5^9\) 的比较运算）。
在 30# 里的 Algo_3 中，从低位到高位分成了五段（\(24|8|8|8|8\)），分别对 \(2*5^9\) 同余，最后累加后取余数。

依此原理，还可以考虑再精简一下，分成四段：\(29|9|9|9\)，
即低 \(29\) 比特直接累加，而后每 \(9\) 比特依次乘以 \(1714662\)、\(2906944\)、\(74078\) 再累加（易证明， \(\sum \leqslant 2 936 365 435 \lt 2^{32}\)，在 \(32\) 位下无溢出），
虽然减少了一次乘法运算，但提取数据比之前的变困难了，需要更多的位运算。

只是呼吸

受到各位高手的刺激，加上本人对大数运算的爱好，经过深入的思考，不断的摸索，终于得到了我自己比较满意的结果。

由于站长已经把路堵死了，本人只有把站长的题目改动一下：

只用4次32位乘法，求 uint64_t 对 10^9 的余数及商。
（这里，我就不管LED灯的问题了）

程序中用到了管理员 liangbch 在 http://bbs.emath.ac.cn/thread-9245-1-1.html 中提供的“规约”技术。


我用的是vc++2008,        在debug环境下测试通过。
以下是程序的源码

1. #include <stdio.h>
   
2. #include <windows.h>
   
3. _declspec(naked)
   
4. int div_mod(UINT32 *dst,UINT32 *scr)
   
5. {
   
6.         _asm
   
7.         {
   
8.                 push                esi
   
9.                 push                edi
   
10.                 push                ebx
    
11.                 push                ebp
    
12.                
    
13.                
    
14.                 mov                        esi,dword ptr[esp+08h+16]
    
15.                 mov                        ebp,dword ptr[esp+04h+16]               
    
16.                
    
17.                 mov                         eax,1152921504///-----这是10^(-9)的高32位
    
18.                 mul                        dword ptr[esi]//64位的低位 乘 10^(-9)的高32位
    
19.                
    
20.                 mov                        edi,edx       //乘积的结果,暂时保存在 edi ebx 中。
    
21.                 mov                        ebx,eax
    
22. 
    
23.                 mov                        eax,2606387915   //-------这是10^(-9)的低32位
    
24.                 mul                        dword ptr[esi+4]// 64位的高位 乘 10^(-9)的低32位
    
25.                
    
26.                 add                        ebx,eax
    
27.                 adc                        edi,edx         //edi 将作为进位，edi  ebx  存储结果
    
28.                
    
29. 
    
30.                 mov                        eax,1152921504//-----这是10^(-9)的高32位
    
31.                 mul                        dword ptr[esi+4]// 64位的高位 乘 10^(-9)的高32位
    
32.        
    
33.                 add                        eax,edi         //把edi加进来，丢弃ebx中的低32位
    
34.                 adc                        edx,0                        //eax  edx 中包含了  原始64位/10^9的商，其中 edx 存储高32位,eax 存储低32位。
    
35.      
    
36.                 shrd                        eax,edx,28//整体右移28位，右移方向 edx->eax->28,右移后 eax 中是商的低32位
    
37.                
    
38.                 //shr                                edx,28     //得到商的高4位
    
39.                 mov                                ecx,eax   //保留商的低32位
    
40.                 //mov                                ebx,edx   //保留商的高32位         
    
41. 
    
42.                 mov                        eax,1000000000
    
43.                 mul                        ecx                                //ecx 商的低32位
    
44.                                                                         //这一步执行后，我们将丢弃edx中的数，只取 eax 中的值。
    
45.                
    
46.                 mov                        esi,dword ptr[esp+08h+16]
    
47.                 mov                        edi,dword ptr[esi]  //edi 原始64位数的低32位,
    
48. 
    
49.                 sub                        edi,eax  // a64 mov 10^9 的结果，保存在 edi 中
    
50. 
    
51.                 //以下语句对 edi 中的余数和 ebx ecx 中的商进行调整,这个方法由管理员 liangbch 在http://bbs.emath.ac.cn/thread-9245-1-1.html 中提供。
    
52.                 mov                        esi,1000000000
    
53.                 mov                        edx,0xffffffff
    
54.                 xor                        eax,eax
    
55. 
    
56.                 cmp                        edi,esi
    
57.                 SETB                        eax
    
58.                 add                        edx,eax
    
59.                 mov                        eax,edx//传送进位
    
60.                 and                        edx,esi
    
61.                 sub                        edi,edx
    
62. 
    
63.                 mov                        dword ptr[ebp],edi  //得到 a64 mod 10^9 (余数)的准确值，返回主函数
    
64. 
    
65.                 //neg                        eax
    
66.                 //add                        ecx,eax
    
67.                 //adc                        ebx,0  // 得到 a64/10^9 (商)的准确值，其中 ebx中存储商的高32位，ecx中存储商的低32位
    
68. 
    
69.                 //如果需要 a64/10^9 (商) 的值，把第 38、40、65、66、67   行的注释取消即可。
    
70.                 pop                        ebp
    
71.                 pop                        ebx
    
72.                 pop                        edi
    
73.                 pop                        esi
    
74.                
    
75. 
    
76.                 ret
    
77.         }
    
78. }
    
79. 
    
80. 
    
81. 
    
82. 
    
83. 
    
84. 
    
85. int main()
    
86. {
    
87.         UINT64 a64;
    
88.         UINT32 scr[2],dst[1];
    
89. 
    
90.         a64=30789000000000;
    
91.                
    
92.            scr[1]=(UINT64)a64/4294967296;
    
93.            scr[0]=(UINT64)a64%4294967296;
    
94.       
    
95.        div_mod(dst,scr);/// dst[] 返回的是 a64 mod 10^9 (余数）。
    
96.                                                    /// 若有需要，可以同时返回 a64/10^9 (商）。
    
97.         printf("dst=%d\n",dst[0]);
    
98.        
    
99.         system( "pause" );
    
100. 
     
101.         return 0;
     
102. 
     
103. }

复制代码