FFT与NTT
学习NTT首先需要一些前置知识,原根及其求解和FFT快速傅里叶变换
NTT与FFT的理论基础完全一样
FFT中求值与插值使用了n次单位根,而NTT则选择了原根
由此与FFT相比,NTT可以取模、避免浮点精度误差、且运算常数小
但也同时因此不能计算浮点系数,且模数也有限制
NTT的膜数一般是
P
=
r
2
k
+
1
P=r2^k+1
P=r2k+1,常取P=998244353,其一个原根为g=3
可以证明,
g
P
−
1
n
g^{\frac{P-1}{n}}
gnP−1与
w
n
=
e
2
π
i
n
w_n=e^{\frac{2\pi i}{n}}
wn=en2πi具有相同的性质
因此直接修改FFT代码中的wn为相应原根即可
const lt P=998244353, G=3, Ginv=332748118;
const int maxn=4000010;
int n,m;
lt A[maxn],B[maxn];
int lim=1,L,R[maxn];
lt qpow(lt a,lt k,lt mod)
{
lt res=1;
while(k){
if(k&1) res=(res*a)%mod;
a=(a*a)%mod; k>>=1;
}
return res;
}
void NTT(lt *a,int opt)
{
for(int i=0;i<lim;++i)
if(i<R[i]) swap(a[i],a[R[i]]);
for(int i=1;i<lim;i<<=1)
{
lt gn=qpow(opt==1?G:Ginv, (P-1)/(i<<1), P);
for(int j=0;j<lim;j+=(i<<1))
{
lt g=1;
for(int k=0;k<i;++k)
{
lt nx=a[j+k],ny=g*a[i+j+k]%P;
a[j+k]=(nx+ny)%P;
a[i+j+k]=(nx-ny+P)%P;
g=g*gn%P;
}
}
}
}
int main()
{
n=read(); m=read();
for(int i=0;i<=n;++i) A[i]=read();
for(int i=0;i<=m;++i) B[i]=read();
while(lim<=n+m) lim<<=1,L++;
for(int i=0;i<lim;++i)
R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
NTT(A,1); NTT(B,1);
for(int i=0;i<lim;++i)
A[i]=A[i]*B[i]%P;
NTT(A,-1);
lt inv=qpow(lim,P-2,P);
for(int i=0;i<=n+m;++i)
printf("%d ",A[i]*inv%P);
return 0;
}
NTT的应用
多项式求逆
洛谷P4238 【模板】多项式乘法逆
题意:给定一个多项式F(x),请求出一个多项式G(x), 满足
F
(
x
)
∗
G
(
x
)
≡
1
(
m
o
d
x
n
)
F(x) * G(x) \equiv 1 \pmod{x^n}
F(x)∗G(x)≡1(modxn),系数对 998244353取模。
假设已知 B ′ B' B′满足 A ( x ) B ′ ( x ) ≡ 1 ( m o d x k ) A(x)B'(x)\equiv 1\pmod{x^k} A(x)B′(x)≡1(modxk),欲求 A ( x ) B ( x ) ≡ 1 ( m o d x 2 k ) A(x)B(x)\equiv 1\pmod{x^{2k}} A(x)B(x)≡1(modx2k)
则有 A ( x ) ( B ( x ) − B ′ ( x ) ) ≡ 0 ( m o d x k ) A(x)(B(x)-B'(x))\equiv 0\pmod{x^{k}} A(x)(B(x)−B′(x))≡0(modxk)
即 B ( x ) − B ′ ( x ) ≡ 0 ( m o d x k ) B(x)-B'(x)\equiv 0\pmod{x^{k}} B(x)−B′(x)≡0(modxk)
两边同时平方得 B ( x ) 2 − 2 B ( x ) B ′ ( x ) + B ′ ( x ) 2 ≡ 0 ( m o d x 2 k ) B(x)^2-2B(x)B'(x)+B'(x)^2\equiv 0\pmod{x^{2k}} B(x)2−2B(x)B′(x)+B′(x)2≡0(modx2k)
两边同时乘A并化简得 B ( x ) − 2 B ′ ( x ) + A ( X ) B ′ ( x ) 2 ≡ 0 ( m o d x 2 k ) B(x)-2B'(x)+A(X)B'(x)^2\equiv 0\pmod{x^{2k}} B(x)−2B′(x)+A(X)B′(x)2≡0(modx2k)
得 B ( x ) ≡ 2 B ′ ( x ) − A ( X ) B ′ ( x ) 2 ( m o d x 2 k ) B(x)\equiv 2B'(x)-A(X)B'(x)^2\pmod{x^{2k}} B(x)≡2B′(x)−A(X)B′(x)2(modx2k)
时间复杂度满足 T ( n ) = T ( n 2 ) + O ( n log n ) = O ( n log n ) T(n)=T(\frac{n}{2})+O(n\log n)=O(n\log n) T(n)=T(2n)+O(nlogn)=O(nlogn)
const lt P=998244353, Grt=3, Ginv=332748118;
const int maxn=500010;
lt F[maxn],G[maxn],t[maxn];
int lim,R[maxn];
lt qpow(lt a,lt k,lt mod)
{
lt res=1;
while(k){
if(k&1) res=(res*a)%mod;
a=(a*a)%mod; k>>=1;
}
return res;
}
void NTT(lt *a,int opt)
{
for(int i=0;i<lim;++i)
if(i<R[i]) swap(a[i],a[R[i]]);
for(int i=1;i<lim;i<<=1)
{
lt gn=qpow(opt==1?Grt:Ginv, (P-1)/(i<<1), P);
for(int j=0;j<lim;j+=(i<<1))
{
lt g=1;
for(int k=0;k<i;++k)
{
lt nx=a[j+k],ny=g*a[i+j+k]%P;
a[j+k]=(nx+ny)%P;
a[i+j+k]=(nx-ny+P)%P;
g=g*gn%P;
}
}
}
lt inv=qpow(lim,P-2,P);
if(opt==-1)
{
for(int i=0;i<lim;++i)
a[i]=a[i]*inv%P;
}
}
void polyInv(lt *A,lt *B,int n)
{
if(n==1)
{
B[0]=qpow(A[0],P-2,P);
return;
}
polyInv(A,B,(n+1)>>1);
lim=1; int L=0;
while(lim<(n<<1)) lim<<=1,L++;
for(int i=0;i<lim;++i)
R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
for(int i=0;i<n;++i) t[i]=A[i];
for(int i=n;i<lim;++i) t[i]=0;
NTT(t,1); NTT(B,1);
for(int i=0;i<lim;++i)
{
B[i]=(((2ll-t[i]*B[i])%P+P)%P) * B[i];
B[i]=(B[i]%P+P)%P;
}
NTT(B,-1);
for(int i=n;i<lim;++i) B[i]=0;
}
int main()
{
int n=read();
for(int i=0;i<n;++i)
F[i]=read();
solve(F,G,n);
for(int i=0;i<n;++i)
printf("%d ",G[i]);
return 0;
}
多项式除法
洛谷P4512 【模板】多项式除法
题意:给定一个 n 次多项式 F(x) 和一个 m 次多项式 G(x) ,请求出多项式 Q(x), R(x),满足以下条件:
Q(x) 次数为 n-m,R(x) 次数小于 m
F(x) = Q(x) * G(x) + R(x)
所有的运算在模998244353意义下进行
定义
A
R
(
x
)
=
x
n
A
(
1
x
)
A_R(x)=x^nA(\frac{1}{x})
AR(x)=xnA(x1),也即将数组翻转
F
(
x
)
=
Q
(
x
)
G
(
x
)
+
R
(
x
)
F(x) = Q(x)G(x) + R(x)
F(x)=Q(x)G(x)+R(x)
F
(
1
x
)
=
Q
(
1
x
)
G
(
1
x
)
+
R
(
1
x
)
F(\frac{1}{x}) = Q(\frac{1}{x})G(\frac{1}{x}) + R(\frac{1}{x})
F(x1)=Q(x1)G(x1)+R(x1)
x
n
F
(
1
x
)
=
x
n
−
m
Q
(
1
x
)
x
m
G
(
1
x
)
+
x
n
−
m
+
1
x
m
−
1
R
(
1
x
)
x^nF(\frac{1}{x}) = x^{n-m}Q(\frac{1}{x})x^mG(\frac{1}{x}) + x^{n-m+1}x^{m-1}R(\frac{1}{x})
xnF(x1)=xn−mQ(x1)xmG(x1)+xn−m+1xm−1R(x1)
F
R
(
x
)
=
Q
R
(
x
)
G
R
(
x
)
+
x
n
−
m
+
1
R
R
(
x
)
F_R(x) = Q_R(x)G_R(x) + x^{n-m+1}R_R(x)
FR(x)=QR(x)GR(x)+xn−m+1RR(x)
F
R
(
x
)
≡
Q
R
(
x
)
G
R
(
x
)
+
x
n
−
m
+
1
R
R
(
x
)
(
m
o
d
x
n
−
m
+
1
)
F_R(x) \equiv Q_R(x)G_R(x) + x^{n-m+1}R_R(x) \pmod { x^{n-m+1}}
FR(x)≡QR(x)GR(x)+xn−m+1RR(x)(modxn−m+1)
Q
R
(
x
)
≡
F
R
(
x
)
G
R
−
1
(
x
)
(
m
o
d
x
n
−
m
+
1
)
Q_R(x) \equiv F_R(x)G_R^{-1}(x) \pmod { x^{n-m+1}}
QR(x)≡FR(x)GR−1(x)(modxn−m+1)
到此可以用多项式求逆求出
G
R
−
1
(
x
)
G_R^{-1}(x)
GR−1(x),然后用多项式乘法计算
Q
R
(
x
)
Q_R(x)
QR(x)
然后
R
(
x
)
=
Q
(
x
)
G
(
x
)
−
F
(
x
)
R(x) = Q(x)G(x) - F(x)
R(x)=Q(x)G(x)−F(x)
const lt P=998244353, Grt=3, Ginv=332748118;
const int maxn=500010;
lt F[maxn],G[maxn],FR[maxn],GR[maxn];
lt GR_inv[maxn],t[maxn],Q[maxn],R[maxn];
int lim,L,rev[maxn];
lt qpow(lt a,lt k,lt mod)
{
lt res=1;
while(k){
if(k&1) res=(res*a)%mod;
a=a*a%mod; k>>=1;
}
return res;
}
void NTT(lt *a,int opt)
{
for(int i=0;i<lim;++i)
if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=1;i<lim;i<<=1)
{
lt gn=qpow(opt==1?Grt:Ginv, (P-1)/(i<<1), P);
for(int j=0;j<lim;j+=(i<<1))
{
lt g=1;
for(int k=0;k<i;++k)
{
lt nx=a[j+k],ny=g*a[i+j+k]%P;
a[j+k]=(nx+ny)%P;
a[i+j+k]=(nx-ny+P)%P;
g=g*gn%P;
}
}
}
lt inv=qpow(lim,P-2,P);
if(opt==-1)
{
for(int i=0;i<lim;++i)
a[i]=a[i]*inv%P;
}
}
void polyInv(lt *A,lt *B,int n)
{
if(n==1)
{
B[0]=qpow(A[0],P-2,P);
return;
}
polyInv(A,B,(n+1)>>1);
lim=1; int L=0;
while(lim<(n<<1)) lim<<=1,L++;
for(int i=0;i<lim;++i)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
for(int i=0;i<n;++i) t[i]=A[i];
for(int i=n;i<lim;++i) t[i]=0;
NTT(t,1); NTT(B,1);
for(int i=0;i<lim;++i)
{
B[i]=(((2ll-t[i]*B[i])%P+P)%P) * B[i];
B[i]=(B[i]%P+P)%P;
}
NTT(B,-1);
for(int i=n;i<lim;++i) B[i]=0;
}
void polyMul(lt* A,lt* B,int n,int m,lt* ans)
{
lim=1; L=0;
while(lim<=n+m) lim<<=1,L++;
for(int i=0;i<lim;++i)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
NTT(A,1); NTT(B,1);
for(int i=0;i<lim;++i)
{
ans[i]=A[i]*B[i]%P;;
A[i]=A[i]*B[i]%P;
}
NTT(A,-1);
NTT(ans,-1);
}
int main()
{
int n=read(),m=read();
for(int i=0;i<=n;++i) F[i]=FR[n-i]=read();
for(int i=0;i<=m;++i) G[i]=GR[m-i]=read();
for(int i=n-m+1;i<=m;++i) GR[i]=0; // mod x^{n-m+1}意义下
polyInv(GR,GR_inv,n-m+1);
polyMul(FR,GR_inv,n,n-m,t);
for(int i=0;i<=n-m;++i)
Q[i]=t[n-m-i];
for(int i=0;i<=n-m;++i)
printf("%lld ",Q[i]);
polyMul(G,Q,m,n-m,t);
for(int i=0;i<m;++i)
R[i]=((F[i]-t[i])%P+P)%P;
printf("\n");
for(int i=0;i<m;++i)
printf("%lld ",R[i]);
return 0;
}