就最基础的3个多项式算法:FFT,NTT 和 拉格朗日插值
普通多项式乘法
FFT
把多项式用DFT转化成点值表示法,然后再用IDFT转化回来。
我们把多项式的项数写成2的幂次(不够就再补几位),DFT利用单位根的优良性质,有(直接写结论)
可以递归求解。然后 IDFT 的时候把 (omega_{n}^{k}) 换成 (omega{n}^{-k}) 即可。最终的答案要除以 (n)。
考虑优化用递推优化常数。
如上是FFT的递归树。观察到底下的所有数的排列时候特征的(二进制位等于其原下标的反转),考虑从下向上递推计算。
其中对于每个数的二进制反转,可以递推 rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1))
,其中 l
为 (n) 的位数。
具体算法步骤
for 每一层(记录长度len)
for 每一块(i)
for 每一个单位根(ω_j)
x=a[i+j], y=ω*a[i+j+len]
a[i+j]=x+y, a[i+j+len]=x-y
#include<bits/stdc++.h>
#define rep(i,a,b) for(register int i=(a);i<=(b);i++)
#define per(i,a,b) for(register int i=(a);i>=(b);i--)
using namespace std;
const int N=3e6+9;
const double pi=acos(-1);
inline int read() {
register int x=0, f=1; register char c=getchar();
while(c<'0'||c>'9') {if(c=='-') f=-1; c=getchar();}
while(c>='0'&&c<='9') {x=(x<<3)+(x<<1)+c-48,c=getchar();}
return x*f;
}
struct cp {
double x,y;
cp(double xx=0,double yy=0) {x=xx,y=yy;}
};
cp operator + (cp a,cp b) {return cp(a.x+b.x,a.y+b.y);}
cp operator - (cp a,cp b) {return cp(a.x-b.x,a.y-b.y);}
cp operator * (cp a,cp b) {return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int n,m,l,lim,r[N];
cp f[N],g[N];
void fft(cp *a,int t) {
rep(i,0,lim-1) if(i<r[i]) swap(a[i],a[r[i]]);
for(int len=1;len<lim;len<<=1) {
cp dw(cos(pi/len),t*sin(pi/len));
for(int i=0;i<lim;i+=(len<<1)) {
cp w(1,0);
for(int j=0;j<len;j++,w=w*dw) {
cp x=a[i+j],y=w*a[i+j+len];
a[i+j]=x+y,a[i+j+len]=x-y;
}
}
}
}
int main() {
n=read(), m=read();
rep(i,0,n) f[i].x=read();
rep(i,0,m) g[i].x=read();
for(lim=1;lim<=n+m;l++,lim<<=1);
rep(i,0,lim-1) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
fft(f,1), fft(g,1);
rep(i,0,lim) f[i]=f[i]*g[i];
fft(f,-1);
rep(i,0,n+m) printf("%d ",(int)(0.5+f[i].x/lim));
return 0;
}
NTT
用原根 (g) 代替 (omega)。(omega_n o g^{frac{p-1}{n}}),(omega_{n}^{k} o g^{frac{p-1}{n} imes k})。逆变换的时候 (g o g^{-1}) 就行了。
常用的 NTT 模数:(998244353, 1004535809, 469762049),原根都是 (3)。
#include<bits/stdc++.h>
#define int long long
#define rep(i,a,b) for(register int i=(a);i<=(b);i++)
#define per(i,a,b) for(register int i=(a);i>=(b);i--)
using namespace std;
const int N=3e6+9,gg=3,mod=998244353,ig=332748118;
inline int read() {
register int x=0, f=1; register char c=getchar();
while(c<'0'||c>'9') {if(c=='-') f=-1; c=getchar();}
while(c>='0'&&c<='9') {x=(x<<3)+(x<<1)+c-48,c=getchar();}
return x*f;
}
int n,m,f[N],g[N],l,lim,r[N];
int ksm(int x,int y,int res=1) {
while(y) {
if(y&1) res=res*x%mod;
y>>=1, x=x*x%mod;
}
return res;
}
int add(int x,int y) {x+=y; return x>=mod?x-mod:x;}
int mns(int x,int y) {x-=y; return x<0?x+mod:x;}
void ntt(int *a,int t) {
rep(i,0,lim-1) if(i<r[i]) swap(a[i],a[r[i]]);
for(int len=1;len<lim;len<<=1) {
int dg=ksm((t>0?gg:ig),(mod-1)/(len*2));
for(int i=0;i<lim;i+=(len<<1)) {
int w=1;
for(int j=0;j<len;j++,w=w*dg%mod) {
int x=a[i+j],y=w*a[i+j+len]%mod;
a[i+j]=add(x,y), a[i+j+len]=mns(x,y);
}
}
}
}
signed main() {
n=read(), m=read();
rep(i,0,n) f[i]=read();
rep(i,0,m) g[i]=read();
for(lim=1;lim<=n+m;l++,lim<<=1);
rep(i,0,lim) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
ntt(f,1), ntt(g,1);
rep(i,0,lim) f[i]=f[i]*g[i]%mod;
ntt(f,-1);
rep(i,0,n+m) printf("%lld ",(f[i]+mod)*ksm(lim,mod-2)%mod);
return 0;
}
P3338 [ZJOI2014]力
我们发现左边是一个卷积形式。然后左右两边分开讨论。令 (E_i=L_i-R_i)。令 (g_i=frac{1}{i^2})。
一个trick,将 q 反转。设 (p_i=q_{n-i}),则有
是卷积形式。设 (R_{i}=S_{n-i}),则有
code: https://www.luogu.com.cn/record/46475365
P4173 残缺的字符串
对于串 (a[0dots m-1]) 和串 (b[0dots m-1]),如果要求匹配,则有
考虑带入式子。设 (f_x) 表示最终以第 (x) 位结尾的答案。设 (a) 为反转后的 (a),则有
记 (A_p(x)=sum a^p_ix^i, B_p(x)=sum b_ix^{i}, C(x)=sum f_ix^i),则有
做 3 次多项式乘法即可。
https://www.luogu.com.cn/record/47534859
拉格朗日插值
普通的拉格朗日插值
给出一个多项式的点值表示法,然后希望能在 (O(n^2)) 的时间内进行插值(即计算出这个多项式)
用高斯消元可以 (O(n^3)) 求解(待定系数法)
拉格朗日插值定理说,设点对为 ((x_i,y_i)),那么多项式为
很妙!!1
洛谷模板 code: https://www.luogu.com.cn/record/46488893
重心拉格朗日插值法
可以动态 (O(n)) 时间加入一个新的点。具体做法就是把第 (i) 个点得贡献提取出来,做到能在新进来一个点时 (O(1)) 增量求出。
维护 (w_i=prod_{i eq j} frac{y_i}{(x_i-x_j)}) 这样每个部分就都能增量求解了。
LOJ模板
code: https://loj.ac/s/1063338
CF622F The Sum of the k-th Powers
求 (sum_{i=0}^{k} i^k)。
有一个简单的小性质,对于一个通项为 (k) 次多项式 (f(x)) 的序列,差分的通项 (g(x)) 是一个 (k-1) 次整式。其逆定理同样适用。考虑 (a_i=sum_{j=1}^{n} i^k) 的差分序列 (a_i=i^k),其通项 (g(x)=x^k) 显然是一个 (k) 次整式。那么我们知道原序列 (a) 的通项 (f(x)) 为一个 (k+1) 次多项式。 我们考虑用拉插来求出这个多项式。由于 (k) 很小,我们插 (k+2) 组 ((i,a_i)) 即可得到这个通项。为了方便,我们选择 (i=1,2,...,n+2),这样可以 (O(k)) 处理出 (a_i) 然后由于插值连续,拉插可以做到 (O(k))。
取值连续的拉插优化:
我们设关于 (n-j) 的前后缀积 (p_i=prod_{j=1}^{i}n-j),(s_i=prod_{j=i}^{k+2}n-j)。
预处理出 (p,s) 和阶乘逆元,然后线性筛求 (i^k),可以做到 (O(k))。但是我懒,所以在处理 (a) 时就不筛了。
code: https://codeforces.com/contest/622/submission/107185130
P4593 [TJOI2018]教科书般的亵渎
考虑把所有怪物的血量画成一张柱形图拼起来。我们发现每一张亵渎等于说消掉一个梯形。再随便玩玩可以发现一共要使用 (m+1) 张亵渎。我们枚举每一张亵渎的使用。
易推除式子
前半部分式子用拉插,后半部分暴力算,(O(m^2))。
code: https://www.luogu.com.cn/record/46532239
CF995F Cowmpany Cowmpensation
插值优化 DP。
考虑设计 (dp) 方程 (f(u,i)) 表示 (u) 子树且 (u) 的值为 (i) 的方案数。考虑前缀和优化,(g(u,i)) 为 (sum_{j=1}^{i} f(u,i)),则有
观察合并的过程,我们发现 (g) 为 (f) 的部分和,然后 (f) 为 (g) 的乘积。所以 (g) 是一个多项式!考虑证明。(f(leaf,*)) 是一个 (0) 次式。然后 (g(x,i)) 是在 (i) 处的点值,转移相当于点值表达式相乘,即多项式的乘法,乘出来是一个 (f) 为 (sz_x-1) 次多项式。所以我们计算出 满足 (xle n) 的 (g(1,x)) 然后插值即可。复杂度 (O(n^2))。
code: https://codeforces.com/contest/995/submission/107289330