多项式全家桶
\(\color{#E74C3C}\text{warning: }\) 目前为止本博客仍是半半半半半成品。
以下是本文参考以及推荐阅读的相关博客\(\%\%\%\)
\(p.s.\) 博客都写得很好！没懂是蒟蒻个人问题
点击查看相关博客$Link$
OI-wiki 快速傅里叶变换：难得感觉有点抽象。
OI - wiki 拉格朗日插值：配合其他博客食用更佳。
傅里叶变换(FFT)学习笔记：从前置芝士讲起的详细学习笔记。
小学生都能看懂的FFT！！！：比较通俗并且不枯燥的讲解。
快速傅里叶变换(FFT)详解：证明很详细。
多项式初步与傅里叶变换（Fourier Transform） : 清晰但是初看有一丢没看懂，可以康康迭代部分。
多项式 I：拉格朗日插值与快速傅里叶变换：感觉已经算是全家桶了（）通过这篇学会了拉插！
从拉插到快速插值求值：超级详细的插值笔记。
概念
-
\(\operatorname{DFT:}\) 系数表达转换为点值表达。
-
\(\operatorname{IDFT:}\) 点值表达转换为系数表达。（ \(\operatorname{DFT}\)的逆运算 ）
拉格朗日插值
考虑构造 \(f_i(x)\) ， 使得 \(f_i(x_i)=1\) 且 \(f_i(x_j)=0\space(j\ne i)\)。
首先考虑满足后者，则有: \(f_i(x)=\prod_{i\ne j}(x-x_j)\)。
再考虑前者，则在两边同乘 \(\frac{y_i}{f_i(x_i)}\) ，可以得到：\(f_i(x)=y_i\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}\)。
多项式 \(f(x)\) 即为：
注意最后一起算分母的逆元，而不是分别求解。
模版题：luogu P4781
点击查看代码
#include<bits/stdc++.h>
#define cs const
#define il inline
#define LL long long
#define pc(i) putchar(i)
#define fo(i,j,k) for(int i=(j);(i)<=(k);++(i))
#define of(i,j,k) for(int i=(j);(i)>=(k);--(i))
using namespace std;
mt19937 srd(time(0));
cs int inf=0x3f3f3f3f,Mod=998244353,N=2003;
int FL,CH;
template<typename T> bool in(T&a)
{
for(FL=1;!isdigit(CH)&&CH!=EOF;CH=getchar())
if(CH=='-')FL=-1;
for(a=0;isdigit(CH);CH=getchar())
a=a*10+CH-'0';
return a*=FL,CH==EOF?0:1;
}
template<typename T,typename...Args>
int in(T&a,Args&...args){return in(a)+in(args...);}
void wt(int x){if(x<0)x=-x,pc('-');if(x>9)wt(x/10);pc(x%10|48);}
int x[N],y[N],n,k; LL ans;
il int qwq(cs int x,cs int y){return (x-y+Mod)%Mod;}
il int qpow(int b,int p)
{
int r=1;
while(p>0)
{
if(p&1) r=1ll*r*b%Mod;
p>>=1,b=1ll*b*b%Mod;
}
return r;
}
il int inv(cs int x){return qpow(x,Mod-2);}
signed main()
{
in(n,k); fo(i,1,n) in(x[i],y[i]);
fo(i,1,n)
{
int s1=y[i]%Mod,s2=1;
fo(j,1,n) if(i!=j)
s1=1ll*s1*qwq(k,x[j])%Mod,s2=1ll*s2*qwq(x[i],x[j])%Mod;
ans=(ans+1ll*s1*inv(s2)%Mod)%Mod;
}
wt(ans);
return 0;
}
横坐标是连续整数的拉格朗日插值
省流：可以 \(O(n)\)。
横坐标为 \(1,\dots ,n\)的插值公式：
例题：CF622F The Sum of the k-th Powers
题目大意：求 \(\sum_{i=1}^ni^k,n\le10^9,k\le10^6\)
答案是一个 \(k+1\) 次的多项式，找 \(k+2\) 个点带进去拉插，证明在这里找。
预处理阶乘和前后缀积。
点击查看代码
#include<bits/stdc++.h>
#define cs const
#define il inline
#define LL long long
#define pc(i) putchar(i)
#define fo(i,j,k) for(int i=(j);(i)<=(k);++(i))
#define of(i,j,k) for(int i=(j);(i)>=(k);--(i))
using namespace std;
mt19937 srd(time(0));
cs int inf=0x3f3f3f3f,Mod=1e9+7,N=1e6+7;
int FL,CH;
template<typename T> bool in(T&a)
{
for(FL=1;!isdigit(CH)&&CH!=EOF;CH=getchar())
if(CH=='-')FL=-1;
for(a=0;isdigit(CH);CH=getchar())
a=a*10+CH-'0';
return a*=FL,CH==EOF?0:1;
}
template<typename T,typename...Args>
int in(T&a,Args&...args){return in(a)+in(args...);}
void wt(int x){if(x<0)x=-x,pc('-');if(x>9)wt(x/10);pc(x%10|48);}
int n,k,ans,fac[N],L[N],R[N],y,num,den;
il int qpow(int b,int p)
{
int re=1;
while(p>0)
{
if(p&1) re=1ll*re*b%Mod;
p>>=1,b=1ll*b*b%Mod;
}
return re;
}
il int inv(cs int x){ return qpow(x,Mod-2); }
signed main()
{
in(n,k),L[0]=R[k+3]=fac[0]=1;
fo(i,1,k+2) L[i]=1ll*L[i-1]*(n-i)%Mod,fac[i]=1ll*fac[i-1]*i%Mod;
of(i,k+2,1) R[i]=1ll*R[i+1]*(n-i)%Mod;
fo(i,1,k+2)
{
y=(y+qpow(i,k))%Mod,num=1ll*L[i-1]*R[i+1]%Mod,
den=1ll*fac[i-1]*((k-i)&1?-1:1)*fac[k+2-i]%Mod;
ans=(ans+1ll*y*num%Mod*inv(den)%Mod)%Mod;
}
wt((ans+Mod)%Mod);
return 0;
}
\(\operatorname{FFT}\)
复数乘法
单位根
- \(\omega_n^k = \cos k \times \frac{2\pi}{n} + i \sin k \times\frac{2\pi}{n}\)
- \(\omega_{2n}^{2k} =\omega_{n}^{k}\)
- \(\omega_n^0=\omega_n^n=1\)
快速傅里叶变换
\(A(x)=a_0+a_1\times x+a_2\times x^2+a_3\times x^3+a_4\times x^4+\dots+a_{n-2}\times x^{n-2}+a_{n-1}\times x^{n-1}\)
按照下标的奇偶性分类：
\(A(x)=(a_0+a_2\times x^2 +a_4\times x^4+\dots+a_{n-2}\times x^{n-2})+(a_1\times x+a_3\times x^3+a_5\times x^5\dots a_{n-1}\times x^{n-1})\)
设：
\(A_1(x)=a_0+a_2\times x+a_4\times x^2+\dots+a_{n-2}\times x^{\frac{n}{2}-1}\)
\(A_2(x)=a_1+a_3\times x+a_5\times x^2+\dots+a_{n-1}\times x^{\frac{n}{2}-1}\)
则有： \(A(x)=A_1(x^2)+xA_2(x^2)\)
将 \(\omega_n^k(k<\frac{n}{2})\) 代入：\(A(\omega_n^k)=A_1({\omega_n^{2k}})+\omega_n^kA_2(\omega_n^{2k})=A_1(\omega_\frac{n}{2}^k)+\omega_n^kA_2(\omega_\frac{n}{2}^{k})\)
将 \(\omega_n^{k+\frac{n}{2}}\) 代入： \(A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})=A_1(\omega_n^{2k})-\omega_n^k A_2(\omega_n^{2k})\)
快速傅里叶逆变换
点值表示法转化为系数表示法。
\((y_0,y_1,\dots,y_{n-1})\) 为 \((a_0,a_1,\dots,a_{n-1})\) 的点值表示法。
设有一向量 \((c_0,c_1,\dots,c_{n-1})\) 满足 \(c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i\) ，即多项式 \(B(x)=y_0+y_1\times x+y_2\times x^2+\dots+y_{n-1}\times x^{n-1}\) 在 \(\omega_n^0,\omega_n^{-1},\omega_n^{-2},\dots,\omega_n^{-(n-1)}\) 处的点值表示。
设 \(S(x)=\sum_{i=0}^{n-1}x^i\)，将 \(\omega_n^k\) 代入：\(S(\omega_n^k)=1+(\omega_n^k)+(\omega_n^k)^2+\dots+(\omega_n^k)^{n-1}\)
当 \(k \neq 0\) 时，同乘 \(\omega_n^k\) ：\(\omega_n^kS(\omega_n^k)=\omega_n^k+(\omega_n^k)+(\omega_n^k)^2+(\omega_n^k)^3+\dots+(\omega_n^k)^n\)
上下两式相减可得：
因此，\(k\neq 0\)时，\(S(\omega_n^k)=0\)，而 \(k=0\) 时 \(S(\omega_n^0)=n\) 。
对于
当 \(j\neq k\) 时，值为 \(0\) ，而 \(j=k\) 时值为 \(n\)。
因此，\(c_k=na_k\) ，\(a_k=\frac{c_k}{n}\)。
迭代实现
在原序列中 \(i\) 与 \(\frac{i}{2}\) 的关系是 ： \(i\) 可以看做是 \(\frac{i}{2}\) 的二进制上的每一位左移一位得来，那么在反转后的数组中就需要右移一位，同时特殊处理一下奇数。
记为 \(rev[i]=(rev[i>>1]>>1)|((rev[i] \& 1)<<(bit-1)),bit=\lceil \log_2n\rceil\)。
\(\operatorname{FFT}\)常数优化论文
P3803 【模板】多项式乘法（FFT）
点击查看代码
#include<bits/stdc++.h>
#define cs const
#define il inline
#define fo(i,j,k) for(int i=(j);(i)<=(k);++(i))
#define of(i,j,k),for(int i=(j);(i)>=(k);--(i))
#define LL long long
using namespace std;
cs int N=3e6+9;
cs double pi=acos(-1);
int FL,CH;
template<typename T>bool in(T&a)
{
for(FL=1;!isdigit(CH)&&CH!=EOF;CH=getchar())
if(CH=='-')FL=-1;
for(a=0;isdigit(CH);CH=getchar())
a=a*10+CH-'0';
return a*=FL,CH==EOF?0:1;
}
template<typename T,typename...Args>
int in(T&a,Args&...args){return in(a)+in(args...);}
struct Complex
{
double x,y;
Complex operator+(const Complex &a) const
{ return {x+a.x,y+a.y}; }
Complex operator-(const Complex &a) const
{ return {x-a.x,y-a.y}; }
Complex operator*(const Complex &a) const
{ return {x*a.x-y*a.y,x*a.y+y*a.x}; }
}A[N],B[N];
int rev[N],bit,tot,n,m;
il void FFT(Complex a[],int inv)
{
fo(i,0,tot-1) if(i<rev[i]) swap(a[i],a[rev[i]]); //求出要迭代的序列
for(int mid=1;mid<tot;mid<<=1)//待合并区间的长度的一半
{
Complex w1=(Complex){cos(pi/mid),inv*sin(pi/mid)};//单位根
for(int i=0;i<tot;i+=mid*2)
{
Complex wk=(Complex){1,0};
for(int j=0;j<mid;++j,wk=wk*w1)//枚举左半部分
{
Complex x=a[i+j],y=wk*a[i+j+mid];
a[i+j]=x+y,a[i+j+mid]=x-y;
}
}
}
}
signed main()
{
in(n,m);
fo(i,0,n) scanf("%lf",&A[i].x);
fo(i,0,m) scanf("%lf",&B[i].x);
while((1<<bit)<n+m+1) ++bit; //bit=ceil(log_2^n)
tot=1<<bit;
fo(i,0,tot-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
FFT(A,1),FFT(B,1);
fo(i,0,tot-1) A[i]=A[i]*B[i];
FFT(A,-1);
fo(i,0,n+m) printf("%d ",(int)(A[i].x/tot+0.5));
return 0;
}
P1919 【模板】A*B Problem 升级版（FFT 快速傅里叶变换）
点击查看代码
#include<bits/stdc++.h>
#define cs const
#define il inline
#define fo(i,j,k) for(int i=(j);(i)<=(k);++(i))
#define of(i,j,k) for(int i=(j);(i)>=(k);--(i))
#define LL long long
using namespace std;
cs int N=3e6+9;
cs double pi=acos(-1);
int FL,CH;
template<typename T>bool in(T&a)
{
for(FL=1;!isdigit(CH)&&CH!=EOF;CH=getchar())
if(CH=='-')FL=-1;
for(a=0;isdigit(CH);CH=getchar())
a=a*10+CH-'0';
return a*=FL,CH==EOF?0:1;
}
template<typename T,typename...Args>
int in(T&a,Args&...args){return in(a)+in(args...);}
struct Complex
{
double x,y;
Complex operator+(cs Complex &a) cs
{ return {x+a.x,y+a.y}; }
Complex operator-(cs Complex &a) cs
{ return {x-a.x,y-a.y}; }
Complex operator*(cs Complex &a) cs
{ return {a.x*x-y*a.y,x*a.y+y*a.x}; }
}A[N],B[N];
char ca[N],cb[N];
int rev[N],bit,tot,n,m,ans[N],id;
il void FFT(Complex a[],cs int inv)
{
fo(i,0,tot-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int mid=1;mid<tot;mid<<=1)
{
Complex w1=(Complex){cos(pi/mid),inv*sin(pi/mid)};
for(int i=0;i<tot;i+=mid*2)
{
Complex wk=(Complex){1,0};
for(int j=0;j<mid;++j,wk=wk*w1)
{
Complex x=a[i+j],y=wk*a[i+j+mid];
a[i+j]=x+y,a[i+j+mid]=x-y;
}
}
}
}
signed main()
{
scanf("%s%s",ca,cb);
n=strlen(ca)-1,m=strlen(cb)-1;
fo(i,0,n) A[i].x=(ca[n-i]^48);
fo(i,0,m) B[i].x=(cb[m-i]^48);
while((1<<bit)<n+m-1) ++bit;
tot=1<<bit;
fo(i,0,tot-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
FFT(A,1),FFT(B,1);
fo(i,0,tot-1) A[i]=A[i]*B[i];
FFT(A,-1);
for(int i=0,x=0;i<tot||x;++i)
x+=A[i].x/tot+0.5,ans[++id]=x%10,x/=10;
while(id>2&&(!ans[id-1])) --id;
of(i,id-1,1) printf("%d",ans[i]);
return 0;
}
浙公网安备 33010602011771号