找回密码
 欢迎注册
查看: 206|回复: 13

[原创] comba乘法(性能优良)

[复制链接]
发表于 2023-11-6 19:04:46 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?欢迎注册

×
本帖最后由 只是呼吸 于 2023-11-6 20:38 编辑

乘法是进行大数运算的基础,一个好的乘法实现,是我们梦寐以求的。
最简单的乘法就是硬乘法了,它也是最基础的乘法,然而硬乘法的时间复杂度是O(n2),当n较大时,这个时间消耗根本不能忍受,于是出现了karatsuba乘法,能将时间复杂度压缩到O(n1.57)。
    由于karatsuba乘法需要调用硬乘法,所以就要想一些办法,把硬乘法做好。
    通常的方案是:
        1)采用comba乘法。
        2)采用汇编语言。
        3)采用较大的进位制。
我在本站大整数乘法代码的帖子,本质就是满足上述3个条件的一个较好的实现。
现在来看,大整数乘法代码还能有提升的空间,经过实践,把[url=]大整数乘法代码[/url]优化如下:

     硬乘法:        乘数   [3] [2] [1] [0]   //[]表示一个limb
               另一个乘数   [3] [2] [1] [0]

                       相乘的结果       [0][3] [0][2] [0][1] [0][0]
                                     [1][3] [1][2] [1][1] [1][0]
                           [2][3] [2][2] [2][1] [2][0]
                  [3][3] [3][2] [3][1] [3][0]
-----------------------------------------------------------
                     7        6       5       4         3       2        1

每一次乘完的结果,不立即写在内存上,而是保存在寄存器中累加。
下面就是优化后的性能良好的comba乘法实现。



编译环境:Windows 10 ,visual studio 2022 ,debug+x86模式。切记是x86模式!
      
      函数原形: void Big_mul( UINT32 *result,const UINT32 * const left,const UINT32 * const right,int LL,int LR)
            其中:*result 乘积的结果。 * left 乘数。*right 另一个乘数。LL, *left的长度,LR, *right的长度。
            进制:采用10000 0000进制。十进制数高位对应数组较大的下标,十进制数的低位对应数组较小的下标。
      特别声明:*result 要用到的存贮空间,是* left 和*right 两个乘数的空间之和的两倍。
                       * left 和*right 两个乘数的空间,要至少保证有3位冗余。即LL+3、LR+3是数组的有效长度,防止访问越界。
                     * left 和*right 两个乘数的长度不能同时超过14700位十进制数。

  1. _declspec(naked)
  2. void re_mul_abs(UINT32 *result, UINT32 *left, UINT32 *right, int LL, int LR)
  3. {
  4. #undef LEFT_REG
  5. #undef RIGHT_REG
  6. #undef RESULT_REG
  7. #define LEFT_REG        esi
  8. #define RIGHT_REG        edi
  9. #define RESULT_REG        ebp

  10.         __asm
  11.         {
  12.                 push        esi
  13.                 push        edi
  14.                 push        ebx
  15.                 push        ebp


  16.                 mov        RIGHT_REG, dword ptr[esp + 0Ch + 16];  right
  17.                 mov        ebx, dword ptr[esp + 14h + 16]; LL
  18. loop01:
  19.                         mov        eax, dword ptr[esp + 14h + 16]
  20.                         sub        eax, ebx

  21.                         mov        ecx, dword ptr[esp + 10h + 16]; LL//ecx中的值是内循环的次数,每次进入内循环时取新值。
  22.                         mov        LEFT_REG, dword ptr[esp + 08h + 16]; left// LEFT--REGk中的值是内循环的乘数,每次进入内循环时取新值。
  23.                         mov        RESULT_REG, dword ptr[esp + 04h + 16]; result//RESULT_REG,结果,每次进入内循环时移位,
  24.                         lea        RESULT_REG, [RESULT_REG + 8 * eax]
  25.                         push                ebx
  26.                         loop00 :
  27.                 push                ecx
  28.                         mov        eax, dword ptr[LEFT_REG]  // 第一次乘
  29.                         mul        dword ptr[RIGHT_REG]

  30.                         add        dword ptr[RESULT_REG], eax
  31.                         adc        dword ptr[RESULT_REG + 4], edx
  32.                         //-------------------------------------------------------------------
  33.                         mov                        eax, dword ptr[LEFT_REG + 4] // 第二次乘(一次)
  34.                         mul                        dword ptr[RIGHT_REG]

  35.                         mov                        ebx, edx
  36.                         mov                        ecx, eax

  37.                         mov                        eax, dword ptr[LEFT_REG]///   第二次乘(二次)
  38.                         mul                        dword ptr[RIGHT_REG + 4]
  39.                         add                        ecx, eax
  40.                         adc                        ebx, edx
  41.                         //------------------------------------------------------------------
  42.                         mov        eax, dword ptr[LEFT_REG]/////  第三次乘( 1次)
  43.                         mul        dword ptr[RIGHT_REG + 8]

  44.                         add                        dword ptr[RESULT_REG + 8], ecx   //
  45.                         adc                        dword ptr[RESULT_REG + 12], ebx  // 第二次乘的结果
  46.                         mov                        ebx, edx
  47.                         mov                        ecx, eax

  48.                         mov        eax, dword ptr[LEFT_REG + 4]/////  第三次乘(二次)
  49.                         mul        dword ptr[RIGHT_REG + 4]
  50.                         add                        ecx, eax
  51.                         adc                        ebx, edx

  52.                         mov        eax, dword ptr[LEFT_REG + 8]/////  第三次乘(三次)
  53.                         mul        dword ptr[RIGHT_REG]
  54.                         add                        ecx, eax
  55.                         adc                        ebx, edx
  56.                         //---------------------------------------------------------
  57.                         mov        eax, dword ptr[LEFT_REG]/////  第四次乘( 1次)
  58.                         mul        dword ptr[RIGHT_REG + 12]

  59.                         add                        dword ptr[RESULT_REG + 16], ecx   //
  60.                         adc                        dword ptr[RESULT_REG + 20], ebx  // 第三次乘的结果
  61.                         mov                        ebx, edx
  62.                         mov                        ecx, eax

  63.                         mov        eax, dword ptr[LEFT_REG + 4]/////  第四次乘( 2次)
  64.                         mul        dword ptr[RIGHT_REG + 8]
  65.                         add                        ecx, eax
  66.                         adc                        ebx, edx

  67.                         mov        eax, dword ptr[LEFT_REG + 8]/////  第四次乘( 3次)
  68.                         mul        dword ptr[RIGHT_REG + 4]
  69.                         add                        ecx, eax
  70.                         adc                        ebx, edx

  71.                         mov        eax, dword ptr[LEFT_REG + 12]/////  第四次乘( 4次)
  72.                         mul        dword ptr[RIGHT_REG]
  73.                         add                        ecx, eax
  74.                         adc                        ebx, edx
  75.                         //------------------------------------------------------
  76.                         mov        eax, dword ptr[LEFT_REG + 4]/////  第五次乘( 1次)
  77.                         mul        dword ptr[RIGHT_REG + 12]

  78.                         add                        dword ptr[RESULT_REG + 24], ecx   //
  79.                         adc                        dword ptr[RESULT_REG + 28], ebx  // 第四次乘的结果
  80.                         mov                        ebx, edx
  81.                         mov                        ecx, eax

  82.                         mov        eax, dword ptr[LEFT_REG + 8]/////  第五次乘( 2次)
  83.                         mul        dword ptr[RIGHT_REG + 8]
  84.                         add                        ecx, eax
  85.                         adc                        ebx, edx

  86.                         mov        eax, dword ptr[LEFT_REG + 12]/////  第五次乘( 3次)
  87.                         mul        dword ptr[RIGHT_REG + 4]
  88.                         add                        ecx, eax
  89.                         adc                        ebx, edx
  90.                         //-----------------------------------------------------------
  91.                         mov        eax, dword ptr[LEFT_REG + 8]/////  第六次乘( 1次)
  92.                         mul        dword ptr[RIGHT_REG + 12]

  93.                         add                        dword ptr[RESULT_REG + 32], ecx   //
  94.                         adc                        dword ptr[RESULT_REG + 36], ebx  // 第五次乘的结果
  95.                         mov                        ebx, edx
  96.                         mov                        ecx, eax

  97.                         mov        eax, dword ptr[LEFT_REG + 12]/////  第六次乘( 2次)
  98.                         mul        dword ptr[RIGHT_REG + 8]
  99.                         add                        ecx, eax
  100.                         adc                        ebx, edx
  101.                         //-----------------------------------------------------------
  102.                         mov        eax, dword ptr[LEFT_REG + 12]/////  第七次乘( 1次)
  103.                         mul        dword ptr[RIGHT_REG + 12]

  104.                         add                        dword ptr[RESULT_REG + 40], ecx   //
  105.                         adc                        dword ptr[RESULT_REG + 44], ebx  // 第六次乘的结果

  106.                         add        dword ptr[RESULT_REG + 48], eax
  107.                         adc        dword ptr[RESULT_REG + 52], edx  //第七次乘的结果
  108.                         //-----------------------------------------------------------   
  109.                         lea        LEFT_REG, [LEFT_REG + 16]
  110.                         lea        RESULT_REG, [RESULT_REG + 32]

  111.                         pop        ecx
  112.                         sub        ecx, 4
  113.                         ja                        loop00;;///内循环。ja  如果ecx为正数则跳转


  114.                 lea        RIGHT_REG, [RIGHT_REG + 16]

  115.                         pop                        ebx
  116.                         sub        ebx, 4
  117.                         ja        loop01//////外循环,ja  如果ebx为正数则跳转

  118.                         //====================================================================//
  119.                     xor                   edi, edi
  120.                         mov        ebx, BASE
  121.                        
  122.                         mov        ecx, dword ptr[esp + 10h + 16]
  123.                         add        ecx, dword ptr[esp + 14h + 16]
  124.                         mov        RESULT_REG, dword ptr[esp + 04h + 16]; result
  125.                         loop02 :
  126.                         xor                   edx, edx
  127.                         xor                   eax, eax

  128.                         add                   edi, dword ptr[RESULT_REG]
  129.                         adc                   eax, dword ptr[RESULT_REG + 4]  ///取本位的高32位做除法
  130.                         div        ebx

  131.                         add        dword ptr [RESULT_REG+12], eax///商加在下一个数的高32位
  132.                         mov        eax, edi  ///取本位数的低32位,与edx中的余数做除法
  133.                         div        ebx

  134.                         mov        dword ptr[RESULT_REG], edx  ///本位的余数就是结果。
  135.                         mov                        edi, eax///本位的商,加在下一位的低32位

  136.                         lea        RESULT_REG, [RESULT_REG + 8]
  137.                         dec        ecx
  138.                         jnz        loop02

  139.         ///---------------------------------------------------------------------------------------------               
  140.                                 mov        ecx, dword ptr[esp + 10h + 16]
  141.                                 add        ecx, dword ptr[esp + 14h + 16]
  142.                                 mov        esi, ecx
  143.                                 mov        edx, 1
  144.                                 xor ebx, ebx
  145.                                 mov        RESULT_REG, dword ptr[esp + 04h + 16]; result
  146.                                 loop04 :
  147.                                 mov       eax, dword ptr[RESULT_REG + 8 * edx]
  148.                                 mov       dword ptr[RESULT_REG + 4 * edx], eax

  149.                                 add       edx, 1
  150.                                 dec       ecx
  151.                                 jnz       loop04////////

  152.                                 lea       RESULT_REG, [RESULT_REG + 4 * edx]

  153.                                 loop05:
  154.                         mov       dword ptr[RESULT_REG], ebx
  155.                                 lea        RESULT_REG, [RESULT_REG + 4]

  156.                                 dec       esi
  157.                                 jnz       loop05//(///  loop 04\05循环,把结果移位,将多余的数清零

  158.                                 mov       dword ptr[RESULT_REG], ebx
  159.                         //-----------------------------------------------------///*///

  160.                         pop        ebp
  161.                         pop        ebx
  162.                         pop        edi
  163.                         pop        esi

  164.                         ret
  165.         }/////*////
  166. }///*////
复制代码

评分

参与人数 1威望 +2 金币 +2 贡献 +2 鲜花 +2 收起 理由
ShuXueZhenMiHu + 2 + 2 + 2 + 2 很给力!

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2023-11-6 20:42:32 | 显示全部楼层
程序第147行。BASE是10^8,进位制的基数。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-11-7 07:56:55 | 显示全部楼层
现在有现成的大整数算法库,比hugecalc还牛逼很多,所以我建议你不要浪费时间搞这个了

点评

就像我的亲弟弟,五十七岁,还整晚在河边钓鱼。  发表于 2023-11-7 09:33
明白,我只是把写程序当作业余爱好,不会去比那些大数库。  发表于 2023-11-7 09:21
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-11-7 09:50:00 | 显示全部楼层
你用什么编辑器来写汇编的?

点评

vs 2022  发表于 2023-11-7 11:57
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-11-7 10:15:22 | 显示全部楼层
只是呼吸 发表于 2023-11-6 20:42
程序第147行。BASE是10^8,进位制的基数。

估计你也有60了,少搞大整数乘法,浪费时间
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-11-7 11:32:33 | 显示全部楼层

点评

nyy
https://bbs.emath.ac.cn/forum.php?mod=redirect&goto=findpost&ptid=129&pid=627 我就找个你给别人泼冷水的帖子让你自己看看  发表于 2023-11-8 09:38
nyy
多赚钱,或者多享受生活不好吗?难道死磕大整数乘法的软件就好?  发表于 2023-11-8 09:34
nyy
浪费时间就是好习性?  发表于 2023-11-8 07:52
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-11-8 10:44:25 | 显示全部楼层
NB!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
回复

使用道具 举报

发表于 2023-11-23 16:21:34 | 显示全部楼层
还是用x86-64位汇编写吧,32位机器面临淘汰了,而且64下寄存器多不少呢,还默认可以用SSE2
当然,vs 64位不支持内嵌汇编可能麻烦些
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2023-12-1 04:55 , Processed in 0.073409 second(s), 17 queries .

Powered by Discuz! X3.5

© 2001-2023 Discuz! Team.

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