吸毒开始了
多项式求逆
对于一个给定的我们要求一个使得
考虑递归求解。
若则就是的逆元。
否则,假设求出了使得
可以大力推一波式子。
设
则的到项都是0
两边平方,则
对于,发现和一定会有一个小于
故有
移项,再两边乘F,可得
于是发现了递归规律,做下去就好了。不过其实不用递归,循环迭代一下就好了。
复杂度类似于等比数列求和,还是nlogn。
有一个比较简单的优化,就是每次NTT的时候不用卷到 ,只需卷到 即可。这样看起来会错,因为它的确会达到 ,但是我们利用了FFT的性质:循环卷积。考虑我们现在需要的数据其实就是 到 的这部分,而错误的数据只会堆积在 到 的部分,所以一点毛病都没有。代码可以看下面多点求值的求逆代码。
有一个更优美且普适性更强的推法:多项式的牛顿迭代。
这个式子可以让我们想到求函数的零点,从而想到牛顿迭代。
一般地,若想求一个多项式F使得G这个函数取到零点,则
推导过程其实可以统统忘记,因为这就是一个牛顿迭代的形式。注意求导时是对F这个多项式整体求导,其他的已知的多项式等可以看成常数。这个与上面推导的结论是一样的,但是更加简单自然。
有一点要注意,对于求 的零点时用的导数也应是 意义下的。如,逆元等。
#include<bits/stdc++.h>
using namespace std;
const int N=4e5+50,mod=998244353,G=3,Gi=332748118;//多项式求逆时N要开到4倍,t可能是n的2倍,q又是t的2倍
int n,a[N],b[N],p[N],c[N];
int power(int x,int y){
int z=1;
for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
return z;
}
void NTT(int *a,int n,int opt){
for(int i=0;i<n;i++)if(i<p[i])swap(a[i],a[p[i]]);
for(int i=1;i<n;i<<=1)
for(int t=i<<1,Wn=power(opt>0?G:Gi,(mod-1)/t),j=0;j<n;j+=t)
for(int k=0,w=1;k<i;k++,w=1ll*w*Wn%mod){
int x=a[j+k],y=1ll*a[i+j+k]*w%mod;
a[j+k]=(x+y)%mod,a[i+j+k]=(x-y)%mod;
}
if(opt<0)for(int i=0,inv=power(n,mod-2);i<n;i++)a[i]=1ll*a[i]*inv%mod;
}
int main(){//多项式有逆元当且仅当常数项有逆元
scanf("%d",&n);
for(int i=0;i<n;i++)scanf("%d",&a[i]);
b[0]=power(a[0],mod-2);//%x时多项式的逆元就是常数项的逆元
for(int k=1,t=2,q=4;(t>>1)<n;t<<=1,k++,q<<=1){//这个是求出在%x^t时的逆元
for(int i=0;i<q;i++)p[i]=p[i>>1]>>1^(i&1)<<k;
for(int i=0;i<t;i++)c[i]=a[i];for(int i=t;i<q;i++)c[i]=0;//t后面的数值不必要,先扔掉,可以使值域缩小
NTT(b,q,1);NTT(c,q,1);//注意因为这个卷积是3个东西卷起来,肯定只算到t是不行的,算到q=2*t就好了
for(int i=0;i<q;i++)b[i]=1ll*b[i]*(2-1ll*b[i]*c[i]%mod)%mod;//因为b和c都扔掉过一部分,所以最高次项是(t>>1)+(t>>1)+t=2*t=q这样(因为b是在上次循环扔掉的)
NTT(b,q,-1);for(int i=t;i<q;i++)b[i]=0;//t后面的数值不必要,同理扔掉
}
for(int i=0;i<n;i++)printf("%d ",(b[i]+mod)%mod);
}//太优美啦,吸毒真爽!(划掉)
分治fft
有一类dp问题,它的转移方程是类似
那么这个dp可以用分治fft或者多项式求逆优化
先说分治fft,我们发现对于一个i,他的f值都是他前面的f值贡献的。
于是可以考虑cdq分治,每次处理一个区间左半段对右半段的影响。
这个直接卷一卷就好了,复杂度是
#include<bits/stdc++.h>
using namespace std;
const int N=4e5+50,mod=998244353,G=3,Gi=332748118;
int n,a[N],p[N],b[N],c[N],d[N];
int power(int x,int y){
int z=1;
for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
return z;
}
void NTT(int *a,int n,int opt){
for(int i=0;i<n;i++)if(i<p[i])swap(a[i],a[p[i]]);
for(int i=1;i<n;i<<=1)
for(int t=i<<1,j=0,Wn=power(opt>0?G:Gi,(mod-1)/t);j<n;j+=t)
for(int k=0,w=1;k<i;k++,w=1ll*w*Wn%mod){
int x=a[j+k],y=1ll*a[i+j+k]*w%mod;
a[j+k]=(x+y)%mod;a[i+j+k]=(x-y)%mod;
}
if(opt<0)for(int i=0,inv=power(n,mod-2);i<n;i++)a[i]=1ll*a[i]*inv%mod;
}
void solve(int l,int r){
if(l==r)return;
int mid=(l+r)>>1;
solve(l,mid);
int t=1,k=-1;
while((t)<r-l+1)t<<=1,k++;
for(int i=1;i<t;i++)p[i]=p[i>>1]>>1^(i&1)<<k;
for(int i=0;l+i<=mid;i++)d[i]=b[l+i];
for(int i=mid-l+1;i<t;i++)d[i]=0;
for(int i=0;i<=r;i++)c[i]=a[i];
for(int i=r+1;i<t;i++)c[i]=0;
NTT(c,t,1);NTT(d,t,1);
for(int i=0;i<t;i++)c[i]=1ll*d[i]*c[i]%mod;
NTT(c,t,-1);
for(int i=1;mid+i<=r;i++)(b[mid+i]+=c[mid-l+i])%=mod;
solve(mid+1,r);
}
int main(){
scanf("%d",&n);
for(int i=1;i<n;i++)scanf("%d",&a[i]);
b[0]=1;solve(0,n-1);
for(int i=0;i<n;i++)printf("%d ",(b[i]+mod)%mod);
return 0;
}
多项式求逆更优秀。设,则对于
一般有初值。
设,即分别为和的生成函数。
那么,移项可得,这样只需把求下逆就好了,复杂度。
据说还有我卷我自己类型的递推式,比如卡特兰数,也可以用分治fft做,但是好像不能求逆。先咕着233
减法卷积
给定和,求,即c的下标要由a和b下标相减得到。
这个形式很容易想到卷积,其实只要把b这个数组reverse之后再处理就好了。
话说现在觉得FFT真好敲
#include<bits/stdc++.h>
using namespace std;
const int N=4e5+50;
const long double Pi=acos(-1);
int n,m,p[N],lim=1,len=-1;
struct node{
long double x,y;
node operator +(const node &b){return (node){x+b.x,y+b.y};}
node operator -(const node &b){return (node){x-b.x,y-b.y};}
node operator *(const node &b){return (node){x*b.x-y*b.y,x*b.y+y*b.x};}
}a[N],b[N];
void FFT(node *a,int opt){
for(int i=0;i<lim;i++)if(i<p[i])swap(a[i],a[p[i]]);
for(int i=1;i<lim;i<<=1){
node Wn=(node){cos(Pi/i),opt*sin(Pi/i)};
for(int t=i<<1,j=0;j<lim;j+=t){
node w=(node){1,0};
for(int k=0;k<i;k++,w=w*Wn){
node x=a[j+k],y=a[i+j+k]*w;
a[j+k]=x+y;a[i+j+k]=x-y;
}
}
}
if(opt<0)for(int i=0;i<lim;i++)a[i].x/=lim;
}
int main(){
scanf("%d%d",&n,&m);
for(int i=0;i<n;i++)scanf("%Lf",&a[i].x);
for(int i=0;i<m;i++)scanf("%Lf",&b[i].x);
reverse(b,b+m);
while(lim<n+m-1)lim<<=1,len++;
for(int i=1;i<lim;i++)p[i]=p[i>>1]>>1^(i&1)<<len;
FFT(a,1);FFT(b,1);
for(int i=0;i<lim;i++)a[i]=a[i]*b[i];
FFT(a,-1);
for(int i=0;i<n+m-1;i++)printf("%d %d\n",i-m+1,(int)(a[i].x+0.5));
}
多项式ln函数
给一个多项式,求出使得
考虑到求的是函数,那么自然可以想到求导。
两边分别求导,运用复合函数的导数公式,则
然后再对做不定积分(相当于求导的逆运算),就可以得到G了。
多项式的导和不定积分可以求出。
void ln(int *a,int *b,int n){
int lim=1,p=-1;
while(lim<n*2)lim<<=1,p++;
ni(a,b,n);dao(a,n);
for(int i=1;i<lim;i++)r[i]=r[i>>1]>>1^(i&1)<<p;
NTT(a,lim,1);NTT(b,lim,1);
for(int i=0;i<lim;i++)b[i]=1ll*a[i]*b[i]%mod;
NTT(b,lim,-1);jifen(b,n);
}
多项式exp函数
给一个多项式,求出使得
不会求,但是我们会求。那么变换一下式子:
这时又是一个求零点的形式,直接套用牛顿迭代。(导数是在分母,翻上去了)
于是直接倍增做即可。复杂度
多项式除法
对于给定的 次多项式 , 次多项式 , 求出使得 成立的 次多项式 和小于 次多项式 。(求出商和余数)
首先显然对于 次多项式 ,为原多项式系数翻转,设为。
然后是一些巧妙的推导。
我们发现要求的 最高次项是 ,故后面那部分可以直接模掉。
然后就是求个逆的事了。有了 我们就可以轻易计算出 。
#include<bits/stdc++.h>
using namespace std;
const int N=4e5+50,mod=998244353,G=3,Gi=332748118;
int n,m,a[N],b[N],ra[N],rb[N],c[N],irb[N],r[N],lim=1,p=-1,ans1[N],ans2[N];
int power(int x,int y){
int z=1;
for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
return z;
}
void NTT(int *a,int n,int opt){
for(int i=0;i<n;i++)if(r[i]>i)swap(a[i],a[r[i]]);
for(int i=1;i<n;i<<=1)
for(int t=i<<1,Wn=power(opt>0?G:Gi,(mod-1)/t),j=0;j<n;j+=t)
for(int k=0,w=1;k<i;k++,w=1ll*w*Wn%mod){
int x=a[j+k],y=1ll*a[i+j+k]*w%mod;
a[j+k]=(x+y)%mod;a[i+j+k]=(x-y)%mod;
}
if(opt<0)for(int i=0,inv=power(n,mod-2);i<n;i++)a[i]=1ll*a[i]*inv%mod;
}
void getinv(int *a,int *b,int n){
b[0]=power(a[0],mod-2);
for(int k=2,q=4,p=1;k<2*n;k<<=1,q<<=1,p++){
for(int i=0;i<k;i++)c[i]=a[i];
for(int i=1;i<q;i++)r[i]=r[i>>1]>>1^(i&1)<<p;
NTT(b,q,1);NTT(c,q,1);
for(int i=0;i<q;i++)b[i]=b[i]*(2-1ll*b[i]*c[i]%mod)%mod;
NTT(b,q,-1);for(int i=k;i<q;i++)b[i]=0;
}
}
int main(){
scanf("%d%d",&n,&m);
for(int i=0;i<=n;i++)scanf("%d",&a[i]),ra[i]=a[i];
for(int i=0;i<=m;i++)scanf("%d",&b[i]),rb[i]=b[i];
reverse(ra,ra+n+1);reverse(rb,rb+m+1);
getinv(rb,irb,n-m+1);
while(lim<=2*(n-m))lim<<=1,p++;for(int i=n-m+1;i<lim;i++)irb[i]=ra[i]=0;
for(int i=1;i<lim;i++)r[i]=r[i>>1]>>1^(i&1)<<p;
NTT(ra,lim,1);NTT(irb,lim,1);
for(int i=0;i<lim;i++)ans1[i]=1ll*ra[i]*irb[i]%mod;
NTT(ans1,lim,-1);reverse(ans1,ans1+n-m+1);
for(int i=n-m+1;i<lim;i++)ans1[i]=0;
for(int i=0;i<=n-m;i++)printf("%d ",(ans1[i]+mod)%mod);puts("");
lim=1,p=-1;while(lim<=2*m)lim<<=1,p++;
for(int i=m+1;i<lim;i++)ans1[i]=0;
for(int i=1;i<lim;i++)r[i]=r[i>>1]>>1^(i&1)<<p;
NTT(b,lim,1);NTT(ans1,lim,1);
for(int i=0;i<lim;i++)ans2[i]=1ll*b[i]*ans1[i]%mod;
NTT(ans2,lim,-1);
for(int i=0;i<m;i++)printf("%d ",(0ll+a[i]-ans2[i]+mod)%mod);
return 0;
}
我们每时每刻一定要清楚自己想要计算的多项式是有没有翻转,在模什么意义下,值域是什么,有没有把不需要的部分清空避免影响卷积。在这个板子里应该是很好的锻炼。
多项式多点求值
给出一个 次多项式 和 个数 ,求 。
多项式科技其实全都是分治。。考虑分治。
对于要求出 到 的点值,构造多项式 ,求出 ,这样对于 。所以成功缩小规模,递归下去就好了。右半部分同理。复杂度 ,为 。关于 的构造,也是分治一下就好了(考虑线段树的结构)。空间是 的。
值得一提的是,当递归到多项式长度小于一个阈值时可以直接 求解,会快不少。这份代码的NTT经过精细卡常。(参考了神鱼的NTT,Orz NaCly_Fish)
#include<bits/stdc++.h>
using namespace std;
const int N=4e5+50,mod=998244353;
int n,m,a[N],x[N],b[N],c[N],rt[N],irt[N];unsigned long long pp[N];
vector<int>v[N*4],t[N*4];
int power(int x,int y){
int z=1;
for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
return z;
}
int r[N],inv[N],lo[N],W[N],Wi[N];
int Glim(int n){int lim=1;while(lim<n)lim<<=1;return lim;}
void getr(int n){for(int i=1;i<n;i++)r[i]=r[i>>1]>>1^(i&1)<<lo[n];}
void init(int n){
int lim=Glim(n);
for(int i=2,j=0;i<=lim;i<<=1,j++)
inv[i]=power(i,mod-2),lo[i]=j;
int w=power(3,(mod-1)/lim);rt[lim>>1]=1;
for(int i=(lim>>1)+1;i<lim;i++)rt[i]=1ll*rt[i-1]*w%mod;
for(int i=(lim>>1)-1;i;i--)rt[i]=rt[i<<1];
for(int i=1;i<lim;i++)irt[i]=power(rt[i],mod-2);
}
void NTT(int *a,int n,int op){
for(int i=0;i<n;i++)pp[i]=a[r[i]];
for(int i=1,*w=op>0?rt:irt;i<n;i<<=1)
for(int t=i<<1,j=0;j<n;j+=t)
for(int k=0;k<i;k++){
int y=pp[i|j|k]*w[i|k]%mod;
pp[i|j|k]=pp[j|k]+mod-y;pp[j|k]+=y;
}
for(int i=0;i<n;i++)a[i]=pp[i]%mod;
if(op<0)for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv[n]%mod;
}
void build(int x,int l,int r){
if(l==r){v[x].resize(2);v[x][0]=mod-::x[l];v[x][1]=1;return;}
int mid=(l+r)>>1,lc=x<<1,rc=x<<1|1;
build(lc,l,mid);build(rc,mid+1,r);
int lim=Glim(r-l+2);
for(int i=0;i<lim;i++)
b[i]=i<v[lc].size()?v[lc][i]:0,
c[i]=i<v[rc].size()?v[rc][i]:0;
getr(lim);NTT(b,lim,1);NTT(c,lim,1);
for(int i=0;i<lim;i++)b[i]=1ll*b[i]*c[i]%mod;
NTT(b,lim,-1);v[x].resize(r-l+2);
for(int i=0;i<v[x].size();i++)v[x][i]=b[i];
}
int d[N];
void Inv(int *a,int *b,int n){
if(n==1){b[0]=power(a[0],mod-2);return;}
int m=(n+1)>>1;Inv(a,b,m);int lim=Glim(n+m);
for(int i=0;i<lim;i++)c[i]=i<n?a[i]:0,d[i]=i<m?b[i]:0;
getr(lim);NTT(c,lim,1);NTT(d,lim,1);
for(int i=0;i<lim;i++)d[i]=1ll*d[i]*(2+mod-1ll*c[i]*d[i]%mod)%mod;
NTT(d,lim,-1);for(int i=m;i<n;i++)b[i]=d[i];
}
int ra[N],rb[N],irb[N],e[N];
void Div(int *a,int *b,int n,int m,int *c,int *d){
if(n<m){for(int i=0;i<m;i++)d[i]=i<=n?a[i]:0;return;}
for(int i=0;i<=n;i++)ra[i]=a[i];reverse(ra,ra+n+1);
for(int i=0;i<=m;i++)rb[i]=b[i];reverse(rb,rb+m+1);
for(int i=m+1;i<=n-m;i++)rb[i]=0;
Inv(rb,irb,n-m+1);int lim=Glim(2*(n-m+1));
for(int i=n-m+1;i<lim;i++)ra[i]=0,irb[i]=0;
getr(lim);NTT(ra,lim,1);NTT(irb,lim,1);
for(int i=0;i<lim;i++)c[i]=1ll*ra[i]*irb[i]%mod;
NTT(c,lim,-1);reverse(c,c+n-m+1);
for(int i=n-m+1;i<max(lim,m);i++)c[i]=0;
lim=Glim(2*m);getr(lim);
for(int i=0;i<lim;i++)d[i]=i<m?b[i]:0,e[i]=i<m?c[i]:0;
NTT(d,lim,1);NTT(e,lim,1);
for(int i=0;i<lim;i++)d[i]=1ll*d[i]*e[i]%mod;
NTT(d,lim,-1);for(int i=0;i<lim;i++)d[i]=i<m?(a[i]-d[i]+mod)%mod:0;
}
int f[N],g[N];
void solve(int x,int l,int r){
if(r-l<=100&&t[x].size()<100){
for(int i=l;i<=r;i++){
int ret=0;
for(int j=0,d=1;j<t[x].size();j++,d=1ll*d*::x[i]%mod)ret=(ret+1ll*d*t[x][j])%mod;
printf("%d\n",(ret+mod)%mod);
}return;
}
if(l==r){printf("%d\n",(t[x][0]+mod)%mod);return;}
int mid=(l+r)>>1,lc=x<<1,rc=x<<1|1;
for(int i=0;i<t[x].size();i++)a[i]=t[x][i];
for(int i=0;i<v[lc].size();i++)b[i]=v[lc][i];
Div(a,b,t[x].size()-1,v[lc].size()-1,f,g);
t[lc].resize(v[lc].size()-1);
for(int i=0;i<t[lc].size();i++)t[lc][i]=g[i];
for(int i=0;i<v[rc].size();i++)b[i]=v[rc][i];
Div(a,b,t[x].size()-1,v[rc].size()-1,f,g);
t[rc].resize(v[rc].size()-1);
for(int i=0;i<t[rc].size();i++)t[rc][i]=g[i];
solve(lc,l,mid);solve(rc,mid+1,r);
}
int main(){
scanf("%d%d",&n,&m);init(2*max(n,m));
for(int i=0;i<=n;i++)scanf("%d",&a[i]);
for(int i=1;i<=m;i++)scanf("%d",&x[i]);
build(1,1,m);t[1].resize(n+1);
for(int i=0;i<=n;i++)t[1][i]=a[i];
solve(1,1,m);
return 0;
}
多项式快速插值
给出 组点值,求出一个 次多项式使得这些点在多项式上。输出系数。
首先最暴力就是高斯消元 。
考虑拉格朗日插值。
我们平时是拿这个求点值的。但是稍微动动脑子就能知道这公式拆开不就是一个多项式么。。于是我们考虑暴力拆开。
枚举外面那层求和,然后发现y和分母构成的系数真的可以暴力算。然而分子呢?好像可以背包?然而每次都背一遍的话就 了。因为每次只是相当于少了一个物品,那么可以考虑求出全部物品的背包,每次再去把这个物品撤回。具体看代码吧
#include<bits/stdc++.h>
using namespace std;
const int N=5050,mod=998244353;
int n,m,x[N],y[N],f[N],dat,g[N],inv[N],ans[N],anss;
int power(int x,int y){
int z=1;
for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
return z;
}
int main(){
scanf("%d%d",&n,&m);
for(int i=1;i<=n;i++)scanf("%d%d",&x[i],&y[i]);
for(int i=1;i<=n;i++)inv[i]=power(x[i],mod-2);
f[0]=1;
for(int i=1;i<=n;i++)for(int j=i;~j;j--)
f[j]=((j?f[j-1]:0)-1ll*x[i]*f[j])%mod;
for(int i=1;i<=n;i++){
if(!x[i])for(int j=0;j<n;j++)g[i]=f[i+1];
else for(int j=0;j<n;j++)g[j]=1ll*((j?g[j-1]:0)-f[j])*inv[i]%mod;
int dat=y[i];
for(int j=1;j<=n;j++)if(j!=i)dat=1ll*dat*power(x[i]-x[j],mod-2)%mod;
for(int j=0;j<n;j++)ans[j]=(ans[j]+1ll*dat*g[j])%mod;
}
for(int i=0,dat=1;i<n;i++,dat=1ll*dat*m%mod)anss=(anss+1ll*dat*ans[i])%mod;
printf("%d",(anss+mod)%mod);
// for(int i=0;i<n;i++)printf("%d ",(ans[i]+mod)%mod);
return 0;
}
注意这份代码其实是 的,因为每次都求了逆元。。。其实可以线性求逆元的。我懒
发现这个做法暴力就暴力在它每次都暴力背包。这种多个单项式相乘的形式显然是可以用多项式科技优化的。。
先考虑分母。我们构造多项式 ,于是 $$\prod_{j!=i}(x_i-x_j)=\lim_{x\to x_i}\frac{G(x)}{x-x_i}$$
因为 ,故也可写成
发现就是 在 处的导。于是求出 的导之后把所有 扔进去多点求值就好了。。好了,现在我们可以轻松地算出前面的系数,拉格朗日的式子变成
其中,前面已经求出。
还是考虑分治,然后化式子。。
设 ,则
连乘的那部分前面已经处理好了。于是皆大欢喜。复杂度还是 。
不过这个其实已经不是人写的了。
不过需要插值的题目一般是元素是多项式,直接计算复杂度太高,就会考虑代入足够的点值计算,再插值还原多项式。
如果瓶颈不在插值就可以写好写的拉格朗日背包。