comba乘法（性能优良）

只是呼吸 发表于 2023-11-6 19:04:46

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

硬乘法：    乘数 //[]表示一个limb
            另一个乘数

                  相乘的结果



-----------------------------------------------------------
                  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位十进制数。

_declspec(naked)
void re_mul_abs(UINT32 *result, UINT32 *left, UINT32 *right, int LL, int LR)
{
#undef LEFT_REG
#undef RIGHT_REG
#undef RESULT_REG
#define LEFT_REG    esi
#define RIGHT_REG    edi
#define RESULT_REG    ebp

__asm
{
push    esi
push    edi
push    ebx
push    ebp

mov    RIGHT_REG, dword ptr;right
mov    ebx, dword ptr; LL
loop01:
         mov    eax, dword ptr
sub    eax, ebx

mov    ecx, dword ptr; LL//ecx中的值是内循环的次数，每次进入内循环时取新值。
mov    LEFT_REG, dword ptr; left// LEFT--REGk中的值是内循环的乘数，每次进入内循环时取新值。
mov    RESULT_REG, dword ptr; result//RESULT_REG，结果，每次进入内循环时移位，
lea    RESULT_REG,
push ebx
loop00 :
push ecx
mov    eax, dword ptr// 第一次乘
mul    dword ptr

add    dword ptr, eax
adc    dword ptr, edx
//-------------------------------------------------------------------
mov eax, dword ptr // 第二次乘（一次）
mul dword ptr

mov ebx, edx
mov ecx, eax

mov eax, dword ptr/// 第二次乘（二次）
mul dword ptr
add ecx, eax
adc ebx, edx
//------------------------------------------------------------------
mov    eax, dword ptr/////第三次乘( 1次）
mul    dword ptr

add dword ptr, ecx //
adc dword ptr, ebx// 第二次乘的结果
mov ebx, edx
mov ecx, eax

mov    eax, dword ptr/////第三次乘（二次）
mul    dword ptr
add ecx, eax
adc ebx, edx

mov    eax, dword ptr/////第三次乘（三次）
mul    dword ptr
add ecx, eax
adc ebx, edx
//---------------------------------------------------------
mov    eax, dword ptr/////第四次乘( 1次）
mul    dword ptr

add dword ptr, ecx //
adc dword ptr, ebx// 第三次乘的结果
mov ebx, edx
mov ecx, eax

mov    eax, dword ptr/////第四次乘（ 2次）
mul    dword ptr
add ecx, eax
adc ebx, edx

mov    eax, dword ptr/////第四次乘（ 3次）
mul    dword ptr
add ecx, eax
adc ebx, edx

mov    eax, dword ptr/////第四次乘（ 4次）
mul    dword ptr
add ecx, eax
adc ebx, edx
//------------------------------------------------------
mov    eax, dword ptr/////第五次乘( 1次）
mul    dword ptr

add dword ptr, ecx //
adc dword ptr, ebx// 第四次乘的结果
mov ebx, edx
mov ecx, eax

mov    eax, dword ptr/////第五次乘（ 2次）
mul    dword ptr
add ecx, eax
adc ebx, edx

mov    eax, dword ptr/////第五次乘（ 3次）
mul    dword ptr
add ecx, eax
adc ebx, edx
//-----------------------------------------------------------
mov    eax, dword ptr/////第六次乘( 1次）
mul    dword ptr

add dword ptr, ecx //
adc dword ptr, ebx// 第五次乘的结果
mov ebx, edx
mov ecx, eax

mov    eax, dword ptr/////第六次乘（ 2次）
mul    dword ptr
add ecx, eax
adc ebx, edx
//-----------------------------------------------------------
mov    eax, dword ptr/////第七次乘( 1次）
mul    dword ptr

add dword ptr, ecx //
adc dword ptr, ebx// 第六次乘的结果

add    dword ptr, eax
adc    dword ptr, edx//第七次乘的结果
//-----------------------------------------------------------
lea    LEFT_REG,
lea    RESULT_REG,

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

lea    RIGHT_REG,

pop ebx
sub    ebx, 4
ja    loop01//////外循环，ja如果ebx为正数则跳转

//====================================================================//
   xor    edi, edi
mov    ebx, BASE

mov    ecx, dword ptr
add    ecx, dword ptr
mov    RESULT_REG, dword ptr; result
loop02 :
xor    edx, edx
xor    eax, eax

add    edi, dword ptr
adc    eax, dword ptr///取本位的高32位做除法
div    ebx

add    dword ptr , eax///商加在下一个数的高32位
mov    eax, edi///取本位数的低32位，与edx中的余数做除法
div    ebx

mov    dword ptr, edx///本位的余数就是结果。
mov edi, eax///本位的商，加在下一位的低32位

lea    RESULT_REG,
dec    ecx
jnz    loop02

///---------------------------------------------------------------------------------------------
mov    ecx, dword ptr
add    ecx, dword ptr
mov    esi, ecx
mov    edx, 1
xor ebx, ebx
mov    RESULT_REG, dword ptr; result
loop04 :
mov    eax, dword ptr
mov    dword ptr, eax

add    edx, 1
dec    ecx
jnz    loop04////////

lea    RESULT_REG,

loop05:
mov    dword ptr, ebx
lea    RESULT_REG,

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

mov    dword ptr, ebx
//-----------------------------------------------------///*///

pop    ebp
pop    ebx
pop    edi
pop    esi

ret
}/////*////
}///*////

只是呼吸 发表于 2023-11-6 20:42:32

程序第147行。BASE是10^8,进位制的基数。

nyy 发表于 2023-11-7 07:56:55

现在有现成的大整数算法库，比hugecalc还牛逼很多，所以我建议你不要浪费时间搞这个了

nyy 发表于 2023-11-7 09:50:00

你用什么编辑器来写汇编的？

nyy 发表于 2023-11-7 10:15:22

只是呼吸发表于 2023-11-6 20:42
程序第147行。BASE是10^8,进位制的基数。

估计你也有60了，少搞大整数乘法，浪费时间

gxqcn 发表于 2023-11-7 11:32:33

泼冷水的坏习性

ShuXueZhenMiHu 发表于 2023-11-8 10:44:25

NB!

无心人 发表于 2023-11-23 16:21:34

还是用x86-64位汇编写吧，32位机器面临淘汰了，而且64下寄存器多不少呢，还默认可以用SSE2
当然，vs 64位不支持内嵌汇编可能麻烦些

只是呼吸 发表于 2023-12-1 10:21:49

无心人发表于 2023-11-23 16:21
还是用x86-64位汇编写吧，32位机器面临淘汰了，而且64下寄存器多不少呢，还默认可以用SSE2
当然，vs 64位不 ...

正在转为\(10^{18}\)进制。我个人不太喜欢\(2^{64}\)进制，一直用10进制。
64位中的寄存器的确好用。

我正在用的vs2022的确不能在c文件中嵌入汇编，但是可以单独添加一个.asm文件，在.asm文件中写汇编实现，与.c文件一起编译，能通过编译器，结果正确。我只需要64位乘法汇编，加法不要汇编问题不大。

我个人经验，64x64to128的乘法一定要用汇编。

无心人 发表于 2023-12-1 11:33:13

只是呼吸发表于 2023-12-1 10:21
正在转为\(10^{18}\)进制。我个人不太喜欢\(2^{64}\)进制，一直用10进制。
64位中的寄存器的确好用。

你正说错了，加法用不用汇编区别很大，因为编译生成的代码不用分支指令做不到获取进位，自己写的优化代码不但可以获取进位，还能做进位加法

页: [1]

数学研发论坛's Archiver

comba乘法（性能优良）