找回密码
 欢迎注册
楼主: gxqcn

[讨论] 不用64位乘除法,求 uint64_t 对 10^9 的余数

[复制链接]
发表于 2017-1-16 14:59:14 | 显示全部楼层
不错。使用模算术的方法,代码简洁。
我的算法仍然是普通的除法运算,很难化简到这么简单。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2017-1-16 17:27:25 | 显示全部楼层
经过变换后,问题的核心转化为计算一个 \(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\) 位下无溢出),
虽然减少了一次乘法运算,但提取数据比之前的变困难了,需要更多的位运算。

点评

方法很赞,共用体用的很巧妙。只需一次移位。最关键的还是转化为5的9次方,这步很难想到。  发表于 2017-1-16 18:40
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-2-1 11:35:35 | 显示全部楼层
受到各位高手的刺激,加上本人对大数运算的爱好,经过深入的思考,不断的摸索,终于得到了我自己比较满意的结果。

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

只用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.                 mov                        eax,2606387915   //-------这是10^(-9)的低32位
  23.                 mul                        dword ptr[esi+4]// 64位的高位 乘 10^(-9)的低32位
  24.                
  25.                 add                        ebx,eax
  26.                 adc                        edi,edx         //edi 将作为进位,edi  ebx  存储结果
  27.                

  28.                 mov                        eax,1152921504//-----这是10^(-9)的高32位
  29.                 mul                        dword ptr[esi+4]// 64位的高位 乘 10^(-9)的高32位
  30.        
  31.                 add                        eax,edi         //把edi加进来,丢弃ebx中的低32位
  32.                 adc                        edx,0                        //eax  edx 中包含了  原始64位/10^9的商,其中 edx 存储高32位,eax 存储低32位。
  33.      
  34.                 shrd                        eax,edx,28//整体右移28位,右移方向 edx->eax->28,右移后 eax 中是商的低32位
  35.                
  36.                 //shr                                edx,28     //得到商的高4位
  37.                 mov                                ecx,eax   //保留商的低32位
  38.                 //mov                                ebx,edx   //保留商的高32位         

  39.                 mov                        eax,1000000000
  40.                 mul                        ecx                                //ecx 商的低32位
  41.                                                                         //这一步执行后,我们将丢弃edx中的数,只取 eax 中的值。
  42.                
  43.                 mov                        esi,dword ptr[esp+08h+16]
  44.                 mov                        edi,dword ptr[esi]  //edi 原始64位数的低32位,

  45.                 sub                        edi,eax  // a64 mov 10^9 的结果,保存在 edi 中

  46.                 //以下语句对 edi 中的余数和 ebx ecx 中的商进行调整,这个方法由管理员 liangbch 在http://bbs.emath.ac.cn/thread-9245-1-1.html 中提供。
  47.                 mov                        esi,1000000000
  48.                 mov                        edx,0xffffffff
  49.                 xor                        eax,eax

  50.                 cmp                        edi,esi
  51.                 SETB                        eax
  52.                 add                        edx,eax
  53.                 mov                        eax,edx//传送进位
  54.                 and                        edx,esi
  55.                 sub                        edi,edx

  56.                 mov                        dword ptr[ebp],edi  //得到 a64 mod 10^9 (余数)的准确值,返回主函数

  57.                 //neg                        eax
  58.                 //add                        ecx,eax
  59.                 //adc                        ebx,0  // 得到 a64/10^9 (商)的准确值,其中 ebx中存储商的高32位,ecx中存储商的低32位

  60.                 //如果需要 a64/10^9 (商) 的值,把第 38、40、65、66、67   行的注释取消即可。
  61.                 pop                        ebp
  62.                 pop                        ebx
  63.                 pop                        edi
  64.                 pop                        esi
  65.                

  66.                 ret
  67.         }
  68. }






  69. int main()
  70. {
  71.         UINT64 a64;
  72.         UINT32 scr[2],dst[1];

  73.         a64=30789000000000;
  74.                
  75.            scr[1]=(UINT64)a64/4294967296;
  76.            scr[0]=(UINT64)a64%4294967296;
  77.       
  78.        div_mod(dst,scr);/// dst[] 返回的是 a64 mod 10^9 (余数)。
  79.                                                    /// 若有需要,可以同时返回 a64/10^9 (商)。
  80.         printf("dst=%d\n",dst[0]);
  81.        
  82.         system( "pause" );

  83.         return 0;

  84. }
复制代码

点评

原来的(UINT64)等于没有写  发表于 2017-2-10 11:10
抱歉:今天检查了一下,第92行应改为:scr[1]=(UINT32)(a64/4294967296); 第93行应改为:scr[0]=(UINT32)(a64/4294967296);  发表于 2017-2-10 11:09
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-3-29 14:46 , Processed in 0.060070 second(s), 16 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表