基于FFT的大整数乘法

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

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

an1an2a2a1a0=a0×100+a1×101+a2×102++an1×10n1\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)=a0x0+a1x1+a2x2++an1xn1A(x)=a_0x^0+a_1x^1+a_2x^2+\cdots+a_{n-1}x^{n-1}

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
// 用于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;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
// 大整数乘法
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秒。大致符合nlog2nnlog_2n的线性对数增长比例

数据/代码下载

下载


基于FFT的大整数乘法
https://blog.hikawa.top/post/fft/
作者
Kayo Hikawa
发布于
2020-01-13
许可协议