快速数论变换
简介
数论变换(Number-theoretic transform, NTT)是 快速傅里叶变换(FFT)在数论基础上的实现。
NTT 解决的是多项式乘法带模数的情况,可以说有些受模数的限制,数也比较大,
定义
数论变换
数论变换 是一种计算卷积(convolution)的快速算法。最常用算法就包括了前文提到的快速傅里叶变换。然而快速傅立叶变换具有一些实现上的缺点,举例来说,资料向量必须乘上复数系数的矩阵加以处理,而且每个复数系数的实部和虚部是一个正弦及余弦函数,因此大部分的系数都是浮点数,也就是说,必须做复数而且是浮点数的运算,因此计算量会比较大,而且浮点数运算产生的误差会比较大。
在数学中,NTT 是关于任意 环 上的离散傅立叶变换(DFT)。在有限域的情况下,通常称为数论变换 (NTT)。
离散傅里叶变换
离散傅里叶变换(Discrete Fourier transform,DFT) 是傅里叶变换在时域和频域上都呈离散的形式,将信号的时域采样变换为其 DTFT 的频域采样。
对于 点序列 ,它的离散傅里叶变换(DFT)为
其中 是自然对数的底数, 是虚数单位。通常以符号 表示这一变换,即
它的 逆离散傅里叶变换(IDFT)为:
可以记为:
实际上,DFT 和 IDFT 变换式中和式前面的归一化系数并不重要。在上面的定义中,DFT 和 IDFT 前的系数分别为 和 。有时我们会将这两个系数都改 。
矩阵公式
由于离散傅立叶变换是一个 线性 算子,所以它可以用矩阵乘法来描述。在矩阵表示法中,离散傅立叶变换表示如下:
生成子群
子群:群 ,满足 ,则 是 的子群
拉格朗日定理: 证明需要用到陪集,得到陪集大小等于子群大小,每个陪集要么不相交要么相等,所有陪集的并是集合 ,那么显然成立。
生成子群: 的生成子群 , 是 的生成元
阶:群 中 的阶是满足 的最小的 ,符号 ,有 ,显然成立。
考虑群
阶就是满足 的最小的 ,
满足 ,对于质数 ,也就是说 结果互不相同。
模 有原根的充要条件 :
离散对数:,
因为 是原根,所以 每 是一个周期,可以取到 的所有元素 对于 是质数时,就是得到 的所有数,就是 到 的映射 离散对数满足对数的相关性质,如
求原根可以证明满足 的最小的 一定是 的约数 对于质数 ,质因子分解 ,若 恒成立, 为 的原根。
NTT
数论变换(NTT)是通过将离散傅立叶变换化为 ,整数模质数 。这是一个 有限域,只要 可除 ,就存在本元 次方根,所以我们有 对于 正整数 。具体来说,对于质数 , 原根 满足 , 将 看做 的等价,则其满足相似的性质,比如
因为这里涉及到数论变化,所以 (为了区分 FFT 中的 ,我们把这里的 称为 )可以比 FFT 中的 大,但是只要把 看做这里的 就行了,能够避免大小问题。
常见的有
就是 的等价
迭代到长度 时 ,或者
下面是一个大数相乘的模板,参考来源。
参考代码
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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102 | #include <algorithm>
#include <bitset>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <iomanip>
#include <iostream>
#include <map>
#include <queue>
#include <set>
#include <string>
#include <vector>
using namespace std;
inline int read() {
int x = 0, f = 1;
char ch = getchar();
while (ch < '0' || ch > '9') {
if (ch == '-') f = -1;
ch = getchar();
}
while (ch <= '9' && ch >= '0') {
x = 10 * x + ch - '0';
ch = getchar();
}
return x * f;
}
void print(int x) {
if (x < 0) putchar('-'), x = -x;
if (x >= 10) print(x / 10);
putchar(x % 10 + '0');
}
const int N = 300100, P = 998244353;
inline int qpow(int x, int y) {
int res(1);
while (y) {
if (y & 1) res = 1ll * res * x % P;
x = 1ll * x * x % P;
y >>= 1;
}
return res;
}
int r[N];
void ntt(int *x, int lim, int opt) {
register int i, j, k, m, gn, g, tmp;
for (i = 0; i < lim; ++i)
if (r[i] < i) swap(x[i], x[r[i]]);
for (m = 2; m <= lim; m <<= 1) {
k = m >> 1;
gn = qpow(3, (P - 1) / m);
for (i = 0; i < lim; i += m) {
g = 1;
for (j = 0; j < k; ++j, g = 1ll * g * gn % P) {
tmp = 1ll * x[i + j + k] * g % P;
x[i + j + k] = (x[i + j] - tmp + P) % P;
x[i + j] = (x[i + j] + tmp) % P;
}
}
}
if (opt == -1) {
reverse(x + 1, x + lim);
register int inv = qpow(lim, P - 2);
for (i = 0; i < lim; ++i) x[i] = 1ll * x[i] * inv % P;
}
}
int A[N], B[N], C[N];
char a[N], b[N];
int main() {
register int i, lim(1), n;
scanf("%s", &a);
n = strlen(a);
for (i = 0; i < n; ++i) A[i] = a[n - i - 1] - '0';
while (lim < (n << 1)) lim <<= 1;
scanf("%s", &b);
n = strlen(b);
for (i = 0; i < n; ++i) B[i] = b[n - i - 1] - '0';
while (lim < (n << 1)) lim <<= 1;
for (i = 0; i < lim; ++i) r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
ntt(A, lim, 1);
ntt(B, lim, 1);
for (i = 0; i < lim; ++i) C[i] = 1ll * A[i] * B[i] % P;
ntt(C, lim, -1);
int len(0);
for (i = 0; i < lim; ++i) {
if (C[i] >= 10) len = i + 1, C[i + 1] += C[i] / 10, C[i] %= 10;
if (C[i]) len = max(len, i);
}
while (C[len] >= 10) C[len + 1] += C[len] / 10, C[len] %= 10, len++;
for (i = len; ~i; --i) putchar(C[i] + '0');
puts("");
return 0;
}
|
参考资料与拓展阅读
build本页面最近更新:2022/5/21 15:16:46,更新历史
edit发现错误?想一起完善? 在 GitHub 上编辑此页!
people本页面贡献者:ChungZH, isdanni, Yukimaikoriya, Enter-tainer, GeZiyue, henryrabbit, Ir1d, ksyx, O-Omega, ouuan, ranwen, shuzhouliu, sshwy, tigerruanyifan, TrisolarisHD, Xeonacid, yifan-r
copyright本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供,附加条款亦可能应用