基于FFT的大整数乘法

关于快速傅里叶变换的基本介绍，可以参考知乎的这篇文章

由十进制数的特点，可以将大整数乘法转换为多项式乘法

$$\overline{a_{n-1}a_{n-2} \cdots a_2a_1a_0} = a_0\times10^0 + a_1\times10^1 + a_2\times10^2 + \cdots + a_{n-1}\times10^{n-1}$$

可将十进制各位数字作为多项式的系数，得到

$$A(x)=a_0x^0+a_1x^1+a_2x^2+\cdots+a_{n-1}x^{n-1}$$

而快速傅里叶变换可以在$O(nlog_2n)$的时间内求出$n$个点组成的时域向量，时域对应相乘，再用逆变换转换为多项式系数表达则可得到结果

完整代码可在文后下载区下载

// 用于fft和ifft的快速傅里叶变换/逆变换基本算法
std::vector<std::complex<double>> BigInteger::fft_algorithm(const std::vector<std::complex<double>> &factors, bool inverse) const
{
	if (factors.size() == 1)
		return factors;

	std::vector<std::complex<double>> odd; // 奇数项系数/结果
	std::vector<std::complex<double>> even; // 偶数项系数/结果

	for (size_t i = 0; i < factors.size(); i++)
	{
		if (i % 2 == 1)
			odd.push_back(factors[i]);
		else
			even.push_back(factors[i]);
	}

	// 递归分别求奇偶项的离散傅里叶变换
	odd = fft_algorithm(odd, inverse);
	even = fft_algorithm(even, inverse);
	size_t n = factors.size();
	std::vector<std::complex<double>> result;

	for (size_t i = 0; i < n; i++)
	{
		std::complex<double> w(cos(2 * PI * i / n), sin(2 * PI * i / n));
		if (inverse) w = std::conj(w);
		result.push_back(even[i % (n / 2)] + w * odd[i % (n / 2)]);
	}

	return result;
}

// 大整数乘法
BigInteger BigInteger::operator*(const BigInteger &number) const
{
	size_t n = 1;
	while (n < number.len + len)
		n <<= 1;

	std::vector<std::complex<double>> a;
	std::vector<std::complex<double>> b;

	// 用两数填满长度为n的vector中，即得到最高次数n-1的多项式系数向量，高次项不足处补0
	for (size_t i = 0; i < len; i++)
		a.push_back(digits[i]);
	for (size_t i = 0; i < number.len; i++)
		b.push_back(number.digits[i]);

	while (a.size() < n)
		a.push_back(0);
	while (b.size() < n)
		b.push_back(0);

	// 离散傅里叶变换
	a = fft(a);
	b = fft(b);

	// 时域相乘，结果存到a中
	for (size_t i = 0; i < n; i++)
	{
		a[i] *= b[i];
	}

	// 离散傅里叶逆变换
	a = ifft(a);

	// 将变换结果转换为BigInteger
	BigInteger result("0", number.len + len);
	int carry = 0;
	for (size_t i = 0; i < number.len + len; i++)
	{
		int temp = (int)(a[i].real() + 0.5) + carry;
		result.digits[i] = temp % 10;
		carry = temp / 10;
	}

	result.len = number.len + len;
	while (result.len > 1 && result.digits[result.len - 1] == 0)
		result.len--;
	result.sign = number.sign * sign;

	return result;
}

参考测试结果

系统环境
CPU：Intel® Core™ i5-8250U CPU @ 1.60GHz
内存：8GB DDR4
OS：Ubuntu 18.04 LTS x64（Windows Subsystem for Linux）
编译器：gcc version 7.4.0 (Ubuntu 7.4.0-1ubuntu1~18.04.1)

使用给定的两个50000位整数相乘的测试数据进行测试，测试通过，结果正确
50000位的两个整数相乘耗时平均为1.93491秒
100000位数字相乘平均耗时3.97488秒；200000位平均耗时8.17346秒。大致符合$nlog_2n$的线性对数增长比例

数据/代码下载

下载