并不对劲的拉格朗日插值
题目大意
给\(n\)(\(n\leq 2000\))个点\((x_1,y_1),...,(x_n,y_n)\)(\(x,y\leq 998244353\)),求多项式\(f(x)\)使\(\forall i\in [1,n],f(x_i) mod 998244353=y_i\)。
题解
结论:\(f(x)=\sum\limits_{i=1}^{n}y_i\frac{\prod\limits_{j=1,j\neq i}^{n}(x-x_j)}{\prod\limits_{j=1,j\neq i}^{n}(x_i-x_j)}\)
该式可以保证当\(x=x_k\)时,\(i=k\) 时 \(\frac{\prod\limits_{j=1,j\neq i}^{n}(x-x_j)}{\prod\limits_{j=1,j\neq i}^{n}(x_i-x_j)}=1\),\(i\neq k\) 时 \(\frac{\prod\limits_{j=1,j\neq i}^{n}(x-x_j)}{\prod\limits_{j=1,j\neq i}^{n}(x_i-x_j)}=0\),符合题目要求。
考虑怎么求\(f(x)\)每一项系数。
\(f(x)=\sum\limits_{i=1}^{n}\frac{y_i}{\prod\limits_{j=1,j\neq i}^{n}(x_i-x_j)}\times(\prod\limits_{j=1,j\neq i}^{n}(x-x_j))\)
前一部分可以直接求,对于后一部分\(\prod\limits_{j=1,j\neq i}^{n}(x-x_j)\):
先求出\(g(x)=\prod\limits_{j=1}^{n}(x-x_j)\),即求出\((x-x_1)(x-x_2)(x-x_3)...(x-x_n)\),\(k\)次项系数相当于在\((-x_1),(-x_2),...,(-x_n)\)中选\(n-k\)个数的所有方案的积的和,这部分可以设\(dp(i,j)\)表示考虑前\(i\)个数,选了\(j\)个数;
接下来考虑如何从\(g(x)\)中去掉\((x-x_i)\)得到\(\prod\limits_{j=1j\neq i}^{n}(x-x_j)\),相当于已知\(g(x),x_i\),求\(h(x)\)使\(h(x)\times(x-x_i)=g(x)\),设\(g(x)=\sum\limits_{j=0}^n a_j\times x^j\),\(h(x)=\sum\limits_{j=0}^{n-1}b_j\times x^j\),就有\(\sum\limits_{j=0}^{n-1}b_j\times x^{j+1}-\sum\limits_{j=0}^{n-1}b_j\times x_j\times x^j=\sum\limits_{j=0}^n a_j\times x^j\),该式无论\(x\)取何值恒成立,就可以得出\(x\)对应次项的系数相等,即\(b_0=-\frac{a_0}{x_i}\),\(\forall j\in[1,n],b_j=-\frac{a_j-b_{j-1}}{x_i}\)。
时间复杂度\(\Theta(n^2)\)。
代码
#include<algorithm>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<ctime>
#include<iomanip>
#include<iostream>
#include<map>
#include<queue>
#include<set>
#include<stack>
#include<vector>
#define rep(i,x,y) for(register int i=(x);i<=(y);++i)
#define dwn(i,x,y) for(register int i=(x);i>=(y);--i)
#define view(u,k) for(int k=fir[u];~k;k=nxt[k])
#define LL long long
#define maxn 2007
using namespace std;
int read()
{
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)&&ch!='-')ch=getchar();
if(ch=='-')f=-1,ch=getchar();
while(isdigit(ch))x=(x<<1)+(x<<3)+ch-'0',ch=getchar();
return x*f;
}
void write(int x)
{
if(x==0){putchar('0'),putchar('\n');return;}
int f=0;char ch[20];
if(x<0)putchar('-'),x=-x;
while(x)ch[++f]=x%10+'0',x/=10;
while(f)putchar(ch[f--]);
putchar('\n');
return;
}
const LL mod=998244353;
int qx[maxn],qy[maxn],a[maxn],n,K,f[maxn];
int mo(int x){return x>=mod?x-mod:x;}
int mul(int x,int y){int res=1;while(y){if(y&1)res=(LL)res*x%mod;x=(LL)x*x%mod,y>>=1;}return res;}
int getf(int x)
{
int ans=0,now=1;
rep(i,0,n-1)ans=mo(ans+(LL)now*a[i]%mod),now=(LL)now*x%mod;
return ans;
}
int main()
{
n=read(),K=read();
rep(i,1,n)qx[i]=read(),qy[i]=read();
f[0]=1;
rep(i,1,n)dwn(j,i,1)
f[j]=mo(f[j]+(LL)f[j-1]*(mod-qx[i])%mod);
reverse(f,f+n+1);
rep(i,1,n)
{
int num=1,nyx=mul(mod-qx[i],mod-2);
rep(j,1,n)if(j!=i)num=(LL)num*mo(qx[i]-qx[j]+mod)%mod;
num=(LL)mul(num,mod-2)*qy[i]%mod;
int lst=0;
rep(j,0,n-1)
{
lst=(LL)mo(f[j]-lst+mod)*nyx%mod,a[j]=mo(a[j]+(LL)lst*num%mod);
}
}
write(getf(K));
return 0;
}
一些感想
终于填上这个坑了。。。