引入
前置知识:阶与原根
由于 FFT 使用了单位根,导致了精度误差,于是我们考虑用原根代替单位根。
过程和 FFT 是一样的,学过 FFT 的可以跳过,放个 FFT传送门。
取模数 p=998244353,g=3(p−1=223×119)。
设 gn=g(p−1)n。
原根的性质
- gkn=gk−1n⋅gn
- g0n=gnn=1
- g2k2n=gkn
- gn2n=gp−12=−1
- gn2+kn=−gkn
于是发现原根和单位根拥有同样的性质。
快速数论变换
使用 gkn(0≤k<n) 这 n 个不同的数来当作点值表示法中的 n 个 x 。
已知
A(x)=n−1∑i=0ai∗xi
设
A1(x)=a0+a2⋅x+a4⋅x2+⋯+an−2xn2−1
A2(x)=a1+a3⋅x+a5⋅x2+⋯+an−1xn2−1
则
A(x)=A1(x2)+xA2(x2)
将 gkn 和 gn2+kn(0≤k<n2) 代入
A(gkn)=A1(gkn2)+gknA2(gkn2)
A(gn2+kn)=A1(gkn2)−gknA2(gkn2)
两个式子 A1() ,A2() 相同。
当计算 [0,n2) 时,[n2,n) 可以直接计算。
每次计算量减小一半,时间复杂度 O(nlogn)
快速数论逆变换
设 C(x)=c0+c1x+c2x2+⋯+cn−1xn−1
设 bi=gin,yi=C(bi)
C(x) 经过 n 个点,(b0,y0)…(bn−1,yn−1)
引理1
设
S(x)=1+x+x2+⋯+xn−1
因为
S(gkn)=g0n+gkn+g2kn+⋯+g(n−1)kn
又
gknS(gkn)=gkn+g2kn+⋯+g(n−1)kn+gnkn
所以
若 k≠0
则
S(gkn)=01−gkn=0
若 k=0
则
S(gkn)=n
引理2
nck=n−1∑i=0yi(g−kn)i
证明
n−1∑i=0yi(g−kn)i
=n−1∑i=0(n−1∑j=0cj(gin)j)(g−kn)i
=n−1∑i=0(n−1∑j=0cj(gjn)i)(g−kn)i
=n−1∑i=0(n−1∑j=0cj(gj−kn)i)
=n−1∑j=0cjn−1∑i=0(gj−kn)i
=n−1∑j=0cjS(gj−kn)
=nck
证毕
设
D(x)=n−1∑i=0yixi
则
ck=D(g−kn)n
所以再做一次快速数论变换即可
代码
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N = 4000010, M = 25, g = 3, gi = 332748118, mod = 998244353;
ll a[N], b[N];
int qpow(int a, int b) {
int res = 1;
while (b) {
if (b & 1) res = 1ll * res * a % mod;
a = 1ll * a * a % mod;
b >>= 1;
}
return res;
}
struct Poly {
int rev[1 << M];
void init(int bit) {
for (int i = 0; i < (1 << bit); i ++ )
rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (bit - 1));
} void ntt(ll *a, int n, int Inv = 0) {
for (int i = 0; i < n; i ++ ) if (i < rev[i]) swap(a[i], a[rev[i]]);
for (int i = 1; i < n; i <<= 1) {
ll gn = qpow(Inv ? gi : g, (mod - 1) / (i << 1));
for (int j = 0; j < n; j += (i << 1)) {
ll gk = 1;
for (int k = 0; k < i; k ++ , gk = gk * gn % mod) {
ll x = a[j + k], y = gk * a[i + j + k] % mod;
a[j + k] = (x + y) % mod, a[i + j + k] = (x - y + mod) % mod;
}
}
} if (Inv) {
ll inv = qpow(n, mod - 2);
for (int i = 0; i < n; i ++ ) a[i] = a[i] * inv % mod;
}
}
} poly;
int main() {
int n, m;
scanf("%d%d", &n, &m);
for (int i = 0; i <= n; i ++ ) scanf("%lld", &a[i]);
for (int i = 0; i <= m; i ++ ) scanf("%lld", &b[i]);
int bit = 0;
while ((1 << bit) <= n + m) bit ++ ;
poly.init(bit);
poly.ntt(a, 1 << bit);
poly.ntt(b, 1 << bit);
for (int i = 0; i < (1 << bit); i ++ ) a[i] = a[i] * b[i] % mod;
poly.ntt(a, 1 << bit, 1);
for (int i = 0; i <= n + m; i ++ ) printf("%lld ", a[i]);
return 0;
}
一打开acw就看见柚子厨🥰