并不对劲的拉格朗日插值

题目大意

给\(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;
}

一些感想

终于填上这个坑了。。。

转载于:https://www.cnblogs.com/xzyf/p/11592260.html