吸毒开始了

多项式求逆

对于一个给定的FF我们要求一个GG使得 FG1(modxn)F*G\equiv1\pmod{x^n}
考虑递归求解。
n=1n=1GG就是F0F_0的逆元。
否则,假设求出了HH使得 FH1(modxn2)F*H\equiv1\pmod{x^{\lceil\frac{n}{2}\rceil}}
可以大力推一波式子。
FH1(modxn2)F*H\equiv1\pmod{x^{\lceil\frac{n}{2}\rceil}}
FG1(modxn)F*G\equiv1\pmod{x^n}
F(GH)0(modxn2)F*(G-H)\equiv0\pmod{x^{\lceil\frac{n}{2}\rceil}}
GH0(modxn2)G-H\equiv0\pmod{x^{\lceil\frac{n}{2}\rceil}}
P=GHP=G-H
PP00n21\lceil\frac{n}{2}\rceil-1项都是0
两边平方,则(P2)x=i=0xPiPxi(P^2)_x=\sum_{i=0}^xP_i*P_{x-i}
对于x<nx<n,发现iixix-i一定会有一个小于n2\lceil\frac{n}{2}\rceil
故有G22GH+H20(modxn)G^2-2*G*H+H^2\equiv0\pmod{x^n}
移项,再两边乘F,可得GH(2FH)(modxn)G\equiv H*(2-F*H)\pmod{x^n}
于是发现了递归规律,做下去就好了。不过其实不用递归,循环迭代一下就好了。
复杂度类似于等比数列求和,还是nlogn。

有一个比较简单的优化,就是每次NTT的时候不用卷到 2n2n,只需卷到 1.5n1.5n 即可。这样看起来会错,因为它的确会达到 2n2n,但是我们利用了FFT的性质:循环卷积。考虑我们现在需要的数据其实就是 0.5n0.5nnn 的这部分,而错误的数据只会堆积在 000.5n0.5n 的部分,所以一点毛病都没有。代码可以看下面多点求值的求逆代码。

有一个更优美且普适性更强的推法:多项式的牛顿迭代。
F1G0(modxn)F^{-1}-G\equiv0\pmod{x^n}这个式子可以让我们想到求函数的零点,从而想到牛顿迭代。
一般地,若想求一个多项式F使得G这个函数取到零点,则

推导过程其实可以统统忘记,因为这就是一个牛顿迭代的形式。注意求导时是对F这个多项式整体求导,其他的已知的多项式等可以看成常数。这个与上面推导的结论是一样的,但是更加简单自然。
有一点要注意,对于求 mod xnmod\ x^n 的零点时用的导数也应是 mod xnmod\ x^n 意义下的。如lnln,逆元等。

#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问题,它的转移方程是类似Fi=j=0i1FjGijF_i=\sum_{j=0}^{i-1}F_jG_{i-j}
那么这个dp可以用分治fft或者多项式求逆优化
先说分治fft,我们发现对于一个i,他的f值都是他前面的f值贡献的。
于是可以考虑cdq分治,每次处理一个区间左半段对右半段的影响。
这个直接卷一卷就好了,复杂度是O(nlog2n)O(nlog^2n)

#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;
}

多项式求逆更优秀。设G0=0G_0=0,则对于i>0i>0

Fi=j=0iFjGijF_i=\sum_{j=0}^iF_jG_{i-j}

F0F_0一般有初值。
A=i0Fixi,B=i0GixiA=\sum_{i\geq0}F_ix^i,B=\sum_{i\geq0}G_ix^i,即分别为FFGG的生成函数。
那么A=AB+F0A=A*B+F_0,移项可得A=F01BA=\frac{F_0}{1-B},这样只需把1B1-B求下逆就好了,复杂度O(nlogn)O(nlogn)
据说还有我卷我自己类型的递推式,比如卡特兰数,也可以用分治fft做,但是好像不能求逆。先咕着233

减法卷积

给定aia_ibib_i,求ci=ai+jbjc_i=\sum a_{i+j}b_j,即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函数

给一个多项式FF,求出GG使得GlnF(modxn)G\equiv lnF\pmod{x^n}
考虑到求的是lnln函数,那么自然可以想到求导。
两边分别求导,运用复合函数的导数公式,则G=FFG'=\frac{F'}{F}
然后再对GG'做不定积分(相当于求导的逆运算),就可以得到G了。
多项式的导和不定积分可以O(n)O(n)求出。

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函数

给一个多项式FF,求出GG使得GeF(modxn)G\equiv e^F\pmod{x^n}
expexp不会求,但是我们会求lnln。那么变换一下式子:
lnGF0(modxn)lnG-F\equiv0\pmod{x^n}
这时又是一个求零点的形式,直接套用牛顿迭代。(导数是1G0\frac{1}{G_0}在分母,翻上去了)
G=G0(1lnG0+F)G=G_0(1-lnG_0+F)
于是直接倍增做即可。复杂度O(nlogn)O(nlogn)

多项式除法

对于给定的 nn 次多项式 A(x)A(x), mm 次多项式 B(x)B(x), 求出使得 A(x)=B(x)C(x)+D(x)A(x)=B(x)*C(x)+D(x) 成立的 nmn-m 次多项式 C(x)C(x) 和小于 mm 次多项式 D(x)D(x)。(求出商和余数)
首先显然对于 nn 次多项式 f(x)f(x),xnf(x1)x^nf(x^{-1})为原多项式系数翻转,设为fR(x)f_R(x)
然后是一些巧妙的推导。

A(x)=B(x)C(x)+D(x)A(x)=B(x)*C(x)+D(x)

xnA(x1)=xmB(x1)xnmC(x1)+xnm+1xm1D(x1)x^nA(x^{-1})=x^mB(x^{-1})*x^{n-m}C(x^{-1})+x^{n-m+1}\cdot x^{m-1}D(x^{-1})

AR(x)=BR(x)CR(x)+xnm+1DR(x)A_R(x)=B_R(x)C_R(x)+x^{n-m+1}D_R(x)

我们发现要求的 C(x)C(x) 最高次项是 xnmx^{n-m},故后面那部分可以直接模掉。
然后就是求个逆的事了。有了 C(x)C(x) 我们就可以轻易计算出 D(x)D(x)

#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;
}

我们每时每刻一定要清楚自己想要计算的多项式是有没有翻转,在模什么意义下,值域是什么,有没有把不需要的部分清空避免影响卷积。在这个板子里应该是很好的锻炼。

多项式多点求值
给出一个 nn 次多项式 A(x)A(x)mm 个数 x0,x1xmx_0,x_1\dots x_m,求 A(x0),A(x1)A(xm)A(x_0),A(x_1)\dots A(x_m)

多项式科技其实全都是分治。。考虑分治。
对于要求出 x1x_1xnx_n 的点值,构造多项式 B(x)=i=1mid(xxi)B(x)=\prod_{i=1}^{mid}(x-x_i),求出 C(x)=A(x)%B(x)C(x)=A(x)\%B(x),这样对于 x1xmid,C(x)=A(x)x_1\dots x_{mid},C(x)=A(x)。所以成功缩小规模,递归下去就好了。右半部分同理。复杂度 T(n)=2T(n/2)+nlognT(n)=2T(n/2)+nlogn,为 O(nlog2n)O(nlog^2n)。关于 B(x)B(x) 的构造,也是分治一下就好了(考虑线段树的结构)。空间是 O(nlogn)O(nlogn) 的。
值得一提的是,当递归到多项式长度小于一个阈值时可以直接 n2n^2 求解,会快不少。这份代码的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;
}

多项式快速插值

给出 nn 组点值,求出一个 n1n-1 次多项式使得这些点在多项式上。输出系数。

首先最暴力就是高斯消元 n3n^3

考虑拉格朗日插值。

f(x)=i=1nyij!=ixxjxixjf(x)=\sum_{i=1}^ny_i\prod_{j!=i}\frac{x-x_j}{x_i-x_j}

我们平时是拿这个求点值的。但是稍微动动脑子就能知道这公式拆开不就是一个多项式么。。于是我们考虑暴力拆开。
枚举外面那层求和,然后发现y和分母构成的系数真的可以暴力算。然而分子呢?好像可以背包?然而每次都背一遍的话就 n3n^3 了。因为每次只是相当于少了一个物品,那么可以考虑求出全部物品的背包,每次再去把这个物品撤回。具体看代码吧

#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;
}

注意这份代码其实是 O(n2logn)O(n^2logn) 的,因为每次都求了逆元。。。其实可以线性求逆元的。我懒

发现这个做法暴力就暴力在它每次都暴力背包。这种多个单项式相乘的形式显然是可以用多项式科技优化的。。

先考虑分母。我们构造多项式 G(x)=i=1n(xxi)G(x)=\prod_{i=1}^n(x-x_i),于是 $$\prod_{j!=i}(x_i-x_j)=\lim_{x\to x_i}\frac{G(x)}{x-x_i}$$
因为 G(xi)=0G(x_i)=0,故也可写成

limxxiG(x)G(xi)xxi\lim_{x\to x_i}\frac{G(x)-G(x_i)}{x-x_i}

发现就是 G(x)G(x)xix_i 处的导。于是求出 G(x)G(x) 的导之后把所有 xix_i扔进去多点求值就好了。。好了,现在我们可以轻松地算出前面的系数,拉格朗日的式子变成

i=1naij!=i(xxj)\sum_{i=1}^na_i\prod_{j!=i}(x-x_j)

其中ai=yij!=i(xixj)a_i=\frac{y_i}{\prod_{j!=i}(x_i-x_j)},前面已经求出。
还是考虑分治,然后化式子。。

fl,r=i=lrailjr&j!=i(xxj)f_{l,r}=\sum_{i=l}^ra_i\prod_{l\leq j\leq r\&j!=i}(x-x_j)

m=l+r2m=\frac{l+r}{2},则

fl,r=i=1maij=m+1n(xxj)+i=m+1naij=1m(xxj)+fl,m+fm+1,rf_{l,r}=\sum_{i=1}^ma_i\prod_{j=m+1}^n(x-x_j)+\sum_{i=m+1}^na_i\prod_{j=1}^m(x-x_j)+f_{l,m}+f_{m+1,r}

连乘的那部分前面已经处理好了。于是皆大欢喜。复杂度还是 O(nlog2n)O(nlog^2n)

不过这个其实已经不是人写的了。
不过需要插值的题目一般是元素是多项式,直接计算复杂度太高,就会考虑代入足够的点值计算,再插值还原多项式。
如果瓶颈不在插值就可以写好写的拉格朗日n2n^2背包。