无符号大整数相乘优化算法及

Win32下无符号大整数相乘优化算法及其C++实现

Lightning[0GiNr]

1、问题的引出：

两个无符号的大整数相乘是一道实践意味很浓的算法题目，这里的“无符号”(unsigned)指的是相乘的两个数都是正数，不需要考虑符号。由于32位计算机没有指令支持128及以上二进制位数的大整数的运算，所以必须自己设计算法来计算。传统的优化算法基本上都是理论层面上的优化，即尽可能地从理论上减少乘法次数，但是往往不能达到预想的优化效果。比方说二分法优化：

将待相乘的整数M分成相等的左右两个部分M1和M2，另一个相乘整数N也同样地分成N1和N2，然后按这样的方法递归分割，直到最后的元素大小小到可以利用CPU指令直接计算为止。这时利用公式M*N = (M1 + M2) * (N1 + N2) = M1*N1 + M1*N2 + M2*N1 + M2*N2结合移位运算再逐层返回得出最终结果。

显然这种算法理论性过强，一来只有当M和N为2的P次方（P为正整数）时的优化才会节省时间，而实际情况下应对随机数据时则会出现大量位移操作，速度不会得到提升；二来使用的递归算法由于调用栈和跳转指令的开销，浪费大量CPU时间；三来这种方法实际上并没有真正地减少乘法次数，因为除了最后一层递归中的乘法可以直接用CPU指令实现，其余各层的乘法由于数值较大仍得另想办法。

由此，我们须从实际出发，探索一些实用的优化方法。

本程序的测试环境为：Windows XP SP2 32bit + 512MB SDRAM + P4 1.80Ghz + VC

2、朴素的算法思路：

为了简易起见，我们先来设计一个朴素的算法。

使用一个DWORD 类型的数组m_buffer作为缓冲区，大小为64，同时声明一个int类型的变量m_nUsed，记录当前缓冲中DWORD使用的个数（即后面所提到的“位数”）。

类的声明如下：

代码清单：BigNumber.cpp

#include

#include

class CBigNumber

{

public:

CBigNumber()

{

memset(this, 0, sizeof(*this));

m_nUsed = 1;

}

CBigNumber& operator = (DWORD dwData);

CBigNumber& operator *= (const CBigNumber& right);

int GetCount() const { return m_nUsed; }

const DWORD* GetBuffer() const { return m_buffer; }

protected:

VOID OffsetAdd(DWORD dwData, int nOffset);

int m_nUsed;

DWORD m_buffer[64];

};

首先是赋值函数，这个函数将一个DWORD类型的整数转化到CBigNumber中。

CBigNumber& CBigNumber::operator = (DWORD dwData)

{

memset(this, 0, sizeof(*this));

m_nUsed = 1;

m_buffer[0] = dwData;

return *this;

}

下面，重载*=运算符，进行大整数乘法运算。本程序使用双重循环，将被乘数的每一个DWORD位与乘数的每一位相乘，将结果错位累加（类似于手算乘法）。

由于两个DWORD相乘后结果可能溢出DWORD的表示范围，所以可以将一个DWORD分成两个部分，如0x12345678分成0x1234和 0x5678，另一个被乘数假设是 0x44445555，则分成0x4444和0x5555。这样，由于(A+B)*(C+D) = AC + (BC + AD) + BD;结果显然是0x1234*0x4444*0x100000000 + (0x5678*0x4444 + 0x1234*0x5555)*0x10000 + 0x5678*0x5555。

考虑上面的式子，0x1234*0x4444*0x100000000必然是溢出DWORD的，正好可以放到下一个 DWORD。(0x5678*0x4444 + 0x1234*0x5555)*0x10000的前16位是溢出的，需要累加到下一位，后16位则不溢出，放到本位。而0x5678*0x5555绝对不会溢出，也要加到本位。

这种方法比较“正规”，稍后我们会介绍更好的办法。

CBigNumber& CBigNumber::operator *= (const CBigNumber& right)

{

// 朴素的乘法

CBigNumber csTemp;

DWORD A = 0, B = 0, C = 0, D = 0;

for(int i = 0; i < right.m_nUsed; i++)

{

for(int n = 0; n < m_nUsed; n++)

{

A = right.m_buffer[i] >> 16;

B = right.m_buffer[i] & 0x0000FFFF;

C = m_buffer[n] >> 16;

D = m_buffer[n] & 0x0000FFFF;

DWORD AC = A * C;

DWORD BD = B * D;

DWORD BC_AD = B * C + A * D;

csTemp.OffsetAdd(BD + ((BC_AD & 0x0000FFFF) << 16), i + n); // 强行截断

csTemp.OffsetAdd(AC + (BC_AD >> 16), i + n + 1);

}

}

*this = csTemp;

return *this;

}