从拉插到快速插值求值

· · 个人记录

标题原来很长的,简写却之后变得很拗口

老是忘记这个东西,干脆写个文章好了……

为了显示点值和未知数的区别,自由元 X 一律大写。

广告 : 多项式计数杂谈

1.拉格朗日插值

P4781 【模板】拉格朗日插值

首先由代数基本定理, n 个不同的点值,在 n-1 次以内确定唯一的一个多项式。(也可以理解为待定系数法的推论)

对于第 k 个点,我们构造一个多项式 F_k(X) ,使得其满足:

①F_k(x_k)=y_k 且 ②(k≠i)\ F_k(x_i)=0。

也就是说,这个多项式专门用来对付第 k 个点值,而在其他的“观测点”上值都为 0 ,无影响。

那么 G(X)=\sum\limits_{i=1}^nF_k(X) 就是满足要求的多项式。(不难验证其满足所有点值的要求)

考虑先用因式定理构造 0 点,如 T_k(X)\prod\limits_{i=1}^n[k≠i](X-x_i) 容易验证其符合要求②。

然后还要令满足 ①F_k(x_k)=y_k ,考虑先把 x_k 代入 T_k(X) ,那么得到:

我们要修正这个多项式,只需要令其除以这个常数,在 $x_k$ 处的函数值就变为了 $1$ ,再乘上 $y_k$ 即可。 令 $F_k(X)=\dfrac{y_kT_k(X)}{T_k(x_k)}$ 就能满足要求了。 总的来说就是 $G(X)=\sum\limits_{k=1}^n\left(y_k\prod\limits_{i=1,k≠i}^n\dfrac{X-x_i}{x_k-x_i}\right)

这些多项式均不超过 n-1 次,那么这个多项式还是唯一的。

下面我们来看看一些具体问题:

CF622F The Sum of the k-th Powers

这是一道很好的例题,题意不再赘述。

容易感性认知,自然数 m 次方和必然是 m+1 次多项式。

为啥嘞?数列微积分

你可以这样理解, \sum\limits_{i=1}^ni^m 类似把 x^m 积分一下,得到的都是 m+1 次多项式。

我们随便找 m+2 个点值就好了,那么选 \Big(p,\sum\limits_{i=1}^pi^m\Big)(p∈[0,m+1]) 就好了。

( 为了方便,下面m是原题中的m+2 )

那么直接套用上述插值,即可达到O(m^2)的复杂度。这样子是跑不过去的。

我们观察这个问题的特点 : 点值都是连续的,而且不需要真的求出这个多项式。

F(X)=\sum\limits_{k=1}^m\left(y_k\prod\limits_{i=1,k≠i}^m\dfrac{X-x_i}{x_k-x_i}\right)

变一下:

{\rm Ans}=F(n)=\sum\limits_{k=1}^m\left(y_k\prod\limits_{i=1,k≠i}^m\dfrac{n-i}{k-i}\right) \dfrac{\prod\limits_{i=1,k≠i}^mn-i}{\prod\limits_{i=1,k≠i}^mk-i}=\dfrac{\prod\limits_{i=1}^{k-1}(n-i)\prod\limits_{i=k+1}^m(n-i)}{\prod\limits_{i=1}^{k-1}(k-i)\prod\limits_{i=k+1}^m(k-i)}

分母上面的东西维护一个前缀积和后缀积即可。

\prod\limits_{i=1}^{k-1}(k-i)=(k-1)(k-2)...2*1=(k-1)! \prod\limits_{i=k+1}^m(k-i)=(-1)(-2)...(k-m+1)(k-m) =(-1)^{n-k}1*2...(m-k-1)(n-k)=(-1)^{m-k}(m-k)!

这样子就可以 O(m) 做了,不过需要线性筛 y ,嫌麻烦的话直接 O(m\log m) 也可以。

比隔壁O(mlogm)MTT第二类斯特林数不知高到哪里去了!

当然,我这个蒟蒻比较懒,所以要线性筛代码还是去膜拜别的大佬吧……

Code:

#include<cstdio>
#define ll long long
#define MaxN 1005000
#define mod 1000000007
using namespace std;
ll powM(ll a,ll t){
  ll ans=1;
  while(t){
    if (t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
ll n,m,ans,y[MaxN],lf[MaxN],rf[MaxN],ifac[MaxN];
int main()
{
  scanf("%I64d%I64d",&n,&m);
  for (int i=1;i<=m+2;i++)
    y[i]=(y[i-1]+powM(i,m))%mod;
  m+=2;
  ifac[0]=ifac[1]=1;
  for (int i=2;i<=m;i++)
    ifac[i]=ifac[mod%i]*(mod-mod/i)%mod;
  for (int i=2;i<=m;i++)
    ifac[i]=ifac[i]*ifac[i-1]%mod;
  lf[0]=1;
  for (int i=1;i<=m;i++)
    lf[i]=lf[i-1]*(n-i)%mod;
  rf[m+1]=1;
  for (int i=m;i;i--)
    rf[i]=rf[i+1]*(n-i)%mod;
  for (int i=1;i<=m;i++)
    ans=(ans+y[i]*lf[i-1]%mod*
         rf[i+1]%mod*ifac[i-1]%mod*
         ((m-i)&1 ? mod-ifac[m-i] : ifac[m-i]))%mod;
  printf("%I64d",(ans+mod)%mod);
  return 0;
}

上面的例题是不需要具体的多项式系数的,一旦需要怎么办呢?

还是看着式子 : F(X)=\sum\limits_{k=1}^m\left(y_k\prod\limits_{i=1,k≠i}^m\dfrac{X-x_i}{x_k-x_i}\right)

设 c_k=y_k*\prod\limits_{i=1,k≠i}^m\dfrac{1}{x_k-x_i} 这是容易求的。

F(X)=\sum\limits_{k=1}^m\left(c_k\prod\limits_{i=1,k≠i}^m(X-x_i)\right)

我们先计算 \prod\limits_{i=1}^m(X-x_i) ,对于每个 \prod\limits_{i=1,k≠i}^m(X-x_i) 只要除去 (X-x_k) 就好。

具体细节请见代码,验板子可以去下方模板的 40' 部分分。

复杂度 O(n^2)。

#include<algorithm>
#include<cstdio>
#define mod 998244353
#define Maxn 5050
#define ll long long 
using namespace std;
ll powM(ll a,ll t=mod-2){
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
int n;
ll x[Maxn],y[Maxn],c[Maxn],fs[Maxn],g[Maxn],f[Maxn];
int main()
{
  scanf("%d",&n);
  for (int i=0;i<n;i++)
    scanf("%lld%lld",&x[i],&y[i]);
  for (int i=0;i<n;i++){
    c[i]=1;
    for (int j=0;j<n;j++)
      if (i!=j)
        c[i]=c[i]*(x[i]-x[j])%mod;
    c[i]=powM(c[i])*y[i]%mod;
  }
  fs[0]=1;
  for (int i=0;i<n;i++){
    for (int j=n;j;j--)
      fs[j]=(fs[j-1]+fs[j]*(mod-x[i]))%mod;
    fs[0]=fs[0]*(mod-x[i])%mod;
  }
  for (int i=0;i<n;i++){
    ll buf=powM(mod-x[i]);
    g[0]=fs[0]*buf%mod;
    for (int j=1;j<n;j++)
      g[j]=(fs[j]-g[j-1])*buf%mod;
    for (int j=0;j<n;j++)
      f[j]=(f[j]+c[i]*g[j])%mod;
  }
  for (int i=0;i<n;i++)
    printf("%lld ",(f[i]+mod)%mod);
  return 0;
}

以下内容可能引起不适……

请先掌握 NTT与多项式全家桶 ,Ln+Exp可以先忽略。

P5050 【模板】多项式多点求值

题意 : 给一个多项式和 x_{1...m} ,欲求 y_{1...m}=F(x_{1...m}).

首先,设 G_i(X)=(X-x_i),把原多项式改写成 F(X)=F_1(X)G_i(X)+F_2(X) ,即多项式取模。

那么,当代入 x_i 的时候, G_i(x_1)=0 ,那么 F(x_i)=F_1(x_i)*0+F_2(x_i)=F_2(x_i)

但是对于每个点值都暴力取模岂不是要 $O(n^2\log n)$ ,比暴力代入还慢呐…… 我们考虑设多项式 $G_{[l,r]}(X)=\prod\limits_{i=l}^rG_i(X)$ 即区间积。 那么, $i∈[l,r]\rightarrow G_i(X)\ |\ G_{[l,r]}(X)$ ,因式关系。 也就是说, $F\bmod G_{[l,r]}\bmod G_i=F\bmod G_i$,先膜上 $G_{[l,r]}$ 对答案没有影响。 但是取模之后次数减少,复杂度就降低了。 先把多项式膜上 $G_{[1,m]}$ ,这样子 $F$ 就是小于 $m$ 次的多项式了,相应地设为 $F_{[1,m]}(X)$。 然后把点值均分成两部分,令 $mid=\frac{l+r}{2}$ ,分别构造 $G_{[1,mid]}$ 和 $G_{[mid+1,m]}$. 然后把当前的 $F$ 分别膜上 $G_{[1,mid]}$ 和 $G_{[mid+1,m]}$ ,产生两个多项式 $F_{[1,mid]};\ F_{[mid+1,n]}$. 那么就变成了两个规模为一半的子问题咯。可以称之为分治膜法。 时间复杂度为 $T(n)=2T(n/2)+O(n\log n)=O(n\log^2n)

代码呢,不那么好写。先要分治NTT先预处理出每一层的 G ,存下来,然后再一层层取膜,所以空间复杂度是 O(n\log n)。

这是我第一次写分治FFT啊……我这种老是爆边界的人自然会把代码写的很丑啦TAT

此外这个模板毒瘤卡常数,记得范围小的时候采用暴力代入减小常数(重要),NTT一定要弄个常数小的来。

#include<algorithm>
#include<cstring>
#include<cstdio>
#define ll long long
#define ull unsigned long long
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
using namespace std;
const int _G=3,mod=998244353,Maxn=66666;
int read(){
  int X=0;char ch=0;
  while(ch<48||ch>57)ch=getchar();
  while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
  return X;
}
ll powM(ll a,ll t=mod-2){
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1],tf;
void tpre(int n){
  if (tf==n)return ;tf=n;
  for(int i=0;i<n;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
  static ull f[Maxn<<1],w[Maxn]={1};
  tpre(n);
  for (int i=0;i<n;i++)f[i]=(((ll)mod<<5)+g[tr[i]])%mod;
  for(int l=1;l<n;l<<=1){
    ull tG=powM(op?_G:invG,(mod-1)/(l+l));
    for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        int tt=w[p]*f[k|l|p]%mod;
        f[k|l|p]=f[k|p]+mod-tt;
        f[k|p]+=tt;
      }      
    if (l==(1<<10))
      for (int i=0;i<n;i++)f[i]%=mod;
  }if (!op){
    ull invn=powM(n);
    for(int i=0;i<n;++i)
      g[i]=f[i]%mod*invn%mod;
  }else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
int _g1[Maxn<<1];
#define sav _g1
void times(int *f,int *g,int len,int lim)
{
  int n=1;while(n<len+len)n<<=1;
  cpy(sav,g,n);clr(sav+len,n-len);
  NTT(f,1,n);NTT(sav,1,n);
  px(f,sav,n);NTT(f,0,n);
  clr(f+lim,n-lim);clr(sav,n);
}
void invp(int *f,int m)
{
  static int w[Maxn<<1],r[Maxn<<1];
  int n;for (n=1;n<m;n<<=1);
  w[0]=powM(f[0]);
  for (int len=2;len<=n;len<<=1){
    for (int i=0;i<(len>>1);i++)
      r[i]=(w[i]<<1)%mod;
    cpy(sav,f,len);
    NTT(w,1,len<<1);px(w,w,len<<1);
    NTT(sav,1,len<<1);px(w,sav,len<<1);
    NTT(w,0,len<<1);clr(w+len,len);
    for (int i=0;i<len;i++)
      w[i]=(r[i]-w[i]+mod)%mod;
  }cpy(f,w,m);clr(sav,n+n);clr(w,n+n);clr(r,n+n);
}
void fan(int *f,int m){
  cpy(sav,f,m);
  for (int i=0;i<m;i++)f[i]=sav[m-i-1];
  clr(sav,m);
}
#undef sav
void savtimes(int *f,int *g,int len1,int len2,int lim)
{
  static int s1[Maxn<<1],s2[Maxn<<1];
  cpy(s1,f,len1);cpy(s2,g,len2);
  times(s1,s2,max(len1,len2),lim);
  cpy(f,s1,lim);clr(s1,lim);clr(s2,len2);
}
void mof(int *f,int *g,int n,int m)
{
  static int q[Maxn<<1],t[Maxn<<1];
  int L=n-m+1;
  fan(g,m);cpy(q,g,L);fan(g,m);
  fan(f,n);cpy(t,f,L);fan(f,n);
  invp(q,L);times(q,t,L,L);fan(q,L);
  savtimes(q,g,L,m,m-1);
  for (int i=0;i<m-1;i++)
    f[i]=(f[i]-q[i]+mod)%mod;
  clr(f+m-1,L);clr(q,n);clr(t,n);
}
int gl[17][Maxn<<1];
void qfpre(int lev,int l,int r,int *x)
{
  if (l==r){
    gl[lev][l<<1]=mod-x[l];
    gl[lev][l<<1|1]=1;
    return ;
  }int mid=(l+r)>>1;
  qfpre(lev+1,l,mid,x);
  qfpre(lev+1,mid+1,r,x);
  cpy(&gl[lev][l<<1],&gl[lev+1][l<<1],mid-l+2);
  savtimes(&gl[lev][l<<1],&gl[lev+1][(mid+1)<<1],mid-l+2,r-mid+1,r-l+2);
}
int _s1[Maxn<<1],_s2[Maxn<<1];
#define s1 _s1
#define s2 _s2
void _queryf(int lev,int l,int r,int *f,int *x,int *y)
{
  if (r-l<=350){
    for (int i=l;i<=r;i++){
      ll buf=1;
      for (int j=l+l;j<l+r+1;j++){
        y[i]=(y[i]+buf*f[j])%mod;
        buf=buf*x[i]%mod;
      }
    }return ;
  }int mid=(l+r)>>1;
  cpy(s1,&f[l<<1],r-l+1);
  mof(s1,&gl[lev+1][l<<1],r-l+1,mid-l+2);
  cpy(s2,&f[l<<1],r-l+1);
  mof(s2,&gl[lev+1][(mid+1)<<1],r-l+1,r-mid+1);
  for (int i=(l<<1);i<(r<<1|1);i++)f[i]=0;
  cpy(&f[l<<1],_s1,mid-l+1);
  cpy(&f[(mid+1)<<1],_s2,r-mid);
  clr(s1,r-l+1);clr(s2,r-l+1);
  _queryf(lev+1,l,mid,f,x,y);
  _queryf(lev+1,mid+1,r,f,x,y);
}
#undef s1
#undef s2
void queryf(int *f,int *x,int *y,int n,int m)
{
  qfpre(0,0,m-1,x);
  if (n>m)mof(f,gl[0],n,m+1);
  _queryf(0,0,m-1,f,x,y);
}
int n,m,f[Maxn<<1],x[Maxn],y[Maxn];
int main()
{
  n=read()+1;m=read();
  for (int i=0;i<n;i++)f[i]=read();
  for (int i=0;i<m;i++)x[i]=read();
  queryf(f,x,y,n,m);
  for (int i=0;i<m;i++)printf("%d\n",y[i]);
  return 0;
}

先来缓口气(?),做做这道卡常题:P5282 【模板】快速阶乘算法

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

P5158 【模板】多项式快速插值

显然直接暴力 O(n^2) 插值是过不去的。

拉格朗日插值公式还是要拿过来用的:

F=\sum\limits_{k=1}^m\left(y_k\prod\limits_{i=1,k≠i}^m\dfrac{X-x_i}{x_k-x_i}\right)

给变得好看一点:

F(X)=\sum\limits_{k=1}^m\left(\dfrac{y_k}{\prod\limits_{i=1,k≠i}^m(x_k-x_i)}\prod\limits_{i=1,k≠i}^m(X-x_i)\right)

现在就能得到F(X)=\sum\limits_{k=1}^m\left(\dfrac{y_k}{G'(x_k)}\prod\limits_{i=1,k≠i}^m(X-x_i)\right)

设V_k=\dfrac{y_k}{G'(x_i)}(反正都是常数嘛)

F(X)=\sum\limits_{k=1}^m\left(V_k\prod\limits_{i=1,k≠i}^m(X-x_i)\right)

这部分考虑分治,设F_{[l,r]}=\sum\limits_{k=1}^r\left(V_k\prod\limits_{i=l,k≠i}^r(X-x_i)\right)

令mid=(l+r)/2,开始分治推式子:

F_{[l,r]}=\sum\limits_{k=l}^r\left(V_k\prod\limits_{i=l,k≠i}^r(X-x_i)\right) =\sum\limits_{k=l}^{mid}\left(V_k\prod\limits_{i=l,k≠i}^r(X-x_i)\right) +\sum\limits_{k=mid+1}^{r}\left(V_k\prod\limits_{i=l,k≠i}^r(X-x_i)\right) =\left[\prod\limits_{i=mid+1}^r(X-x_i)\right]\sum\limits_{k=l}^{mid}\left(V_k\prod\limits_{i=l,k≠i}^{mid}(X-x_i)\right) +\left[\prod\limits_{i=l}^{mid}(X-x_i)\right]\sum\limits_{k=mid+1}^{r}\left(V_k\prod\limits_{i=mid+1,k≠i}^{r}(X-x_i)\right) 边界就是$F_{[k,k]}=V_k=\dfrac{y_k}{G'(x_k)}

那么就变成两个规模为一半的子问题了,分治NTT上!

复杂度O(nlog^2n),常数极大。

如何减小常数呢?G',多点求值,分治NTT中的G都是同一个G,那么可以只算一次。

此外,分治范围小到一定程度的时候可以暴力拉插公式(这个优化我没写,因为其实求值才是瓶颈)。

卡到最后当然是预处理单位根啦。

#include<algorithm>
#include<cstdio>
#define mod 998244353
#define G 3
#define Maxn 140000
#define ll long long 
using namespace std;
inline int read()
{
  register int X=0;
  register char ch=0;
  while(ch<48||ch>57)ch=getchar();
  while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
  return X;
}
int r[Maxn<<1],rn;
ll invn,invG;
ll powM(ll a,ll t=mod-2)
{
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;
    t>>=1;
  }return ans;
}
//f=g
inline void cop(ll *f,ll *g,int len)
{for (int i=0;i<len;i++)f[i]=g[i];}
ll pG[2][19][Maxn<<2];
void preG()
{
  for(int d=1;d<=17;d++){
    ll buf=powM(G,(mod-1)/(1<<d));
    pG[1][d][0]=1;
    for(int i=1;i<(1<<d);i++)
      pG[1][d][i]=pG[1][d][i-1]*buf%mod;
  }invG=powM(G);
  for(int d=1;d<=17;d++){
    ll buf=powM(invG,(mod-1)/(1<<d));
    pG[0][d][0]=1;
    for(int i=1;i<(1<<d);i++)
      pG[0][d][i]=pG[0][d][i-1]*buf%mod;
  }
}
void NTT(ll *f,int n,int op)
{
  for (int i=0;i<n;i++)
    if (r[i]<i)swap(f[r[i]],f[i]);
  ll *buf;
  for (int p=2,j=1;p<=n;p<<=1,j++){
    int len=p>>1;
    for (int k=0;k<n;k+=p){
      buf=pG[op][j];
      for (int i=k;i<k+len;i++,buf++){
        int sav=f[len+i]*(*buf)%mod;
        f[len+i]=f[i]-sav;
        if (f[len+i]<0)f[len+i]+=mod;
        f[i]=f[i]+sav;
        if (f[i]>=mod)f[i]-=mod;
      }
    }
  }
}
ll _g[Maxn<<1];//f*=g
void times(ll *f,ll *gg,int len1,int len2,int limit)
{
  int n=1;for(;n<len1+len2;n<<=1);
  ll *g=_g;
  for (int i=0;i<len2;i++)g[i]=gg[i];
  if (rn!=n){
    rn=n;
    for(int i=0;i<n;i++)
      r[i]=(r[i>>1]>>1)|((i&1)?n>>1:0);
  }NTT(f,n,1);NTT(g,n,1);
  for(int i=0;i<n;++i)f[i]=f[i]*g[i]%mod;
  NTT(f,n,0);invn=powM(n);
  for(int i=0;i<limit;++i)f[i]=f[i]*invn%mod;
  for(int i=limit;i<n;++i)f[i]=0;
  for (int i=0;i<n;i++)g[i]=0;
}
ll _sav[Maxn<<1];//f*=g
inline void savtimes(ll *f,ll *g,int len1,int len2,int limit)
{
  for (int i=0;i<len1;i++)_sav[i]=f[i];
  times(_sav,g,len1,len2,limit);
  for(int i=0;i<limit;++i)f[i]=_sav[i];
  for(int i=0;i<limit;++i)_sav[i]=0;
}
void dao(ll *f,int n)
{for (int i=1;i<n;i++)f[i-1]=f[i]*i%mod;f[n]=0;}
ll _r[Maxn<<1],_rr[Maxn<<1];
void invp(ll *f,int len)
{
  ll *r=_r,*rr=_rr;
  int n=1;for(;n<len;n<<=1);
  rr[0]=powM(f[0]);
  for (int len=2;len<=n;len<<=1){
    for (int i=0;i<len;i++)
      r[i]=(rr[i]<<1)%mod;
    times(rr,rr,len>>1,len>>1,len);
    times(rr,f,len,len,len);
    for (int i=0;i<len;i++)
      rr[i]=(r[i]-rr[i]+mod)%mod;
  }for (int i=0;i<len;i++)f[i]=rr[i];
  for (int i=0;i<n;i++)r[i]=rr[i]=0;
}
inline void fan(ll *f,int m)
{
  for (int i=0;i<m;i++)_sav[i]=f[i];
  for (int i=0;i<m;i++)f[i]=_sav[m-i-1];
  for (int i=0;i<m;i++)_sav[i]=0;
}
ll _q[Maxn<<1],_t[Maxn<<1];//f%=g
void mof(ll *f,ll *g,int n,int m)
{
  ll *q=_q,*t=_t;
  fan(g,m);cop(q,g,n-m+1);fan(g,m);
  invp(q,n-m+1);
  fan(f,n);cop(t,f,n-m+1);fan(f,n);
  times(q,t,n-m+1,n-m+1,n-m+1);
  fan(q,n-m+1);
  times(q,g,n-m+1,m,n);
  for (int i=0;i<m-1;i++)f[i]=(f[i]-q[i]+mod)%mod;
  for (int i=m-1;i<n;i++)f[i]=0;
  for (int i=0;i<n;i++)q[i]=t[i]=0;
}
ll gl[18][Maxn<<1];
void qfpre(int lev,int l,int r,ll *x)
{
  if (l==r){
    gl[lev][l<<1]=mod-x[l];
    gl[lev][l<<1|1]=1;
    return ;
  }int mid=(l+r)>>1;
  qfpre(lev+1,l,mid,x);
  qfpre(lev+1,mid+1,r,x);
  cop(&gl[lev][l<<1],&gl[lev+1][l<<1],mid-l+2);
  savtimes(&gl[lev][l<<1],&gl[lev+1][(mid+1)<<1],mid-l+2,r-mid+1,r-l+2);
}
ll _s1[Maxn<<1],_s2[Maxn<<1];
void queryf(int lev,int l,int r,ll *f,ll *x,ll *y)
{
  if (r-l<=800){
    for (int i=l;i<=r;i++){
      register ll buf=1;
      for (int j=l+l;j<l+r+1;j++){
        y[i]=(y[i]+buf*f[j])%mod;
        buf=buf*x[i]%mod;
      }
    }return ;
  }int mid=(l+r)>>1;
  cop(_s1,&f[l<<1],r-l+1);
  mof(_s1,&gl[lev+1][l<<1],r-l+1,mid-l+2);
  cop(_s2,&f[l<<1],r-l+1);
  mof(_s2,&gl[lev+1][(mid+1)<<1],r-l+1,r-mid+1);
  for (int i=(l<<1);i<(r<<1|1);i++)f[i]=0;
  cop(&f[l<<1],_s1,mid-l+1);
  cop(&f[(mid+1)<<1],_s2,r-mid);
  for (int i=0;i<r-l+1;i++)_s1[i]=_s2[i]=0;
  queryf(lev+1,l,mid,f,x,y);
  queryf(lev+1,mid+1,r,f,x,y);
}
void polf(int lev,int l,int r,ll *f,ll *x,ll *v)
{
  if (l==r){f[l<<1]=v[l];return ;}
  int mid=(l+r)>>1;
  polf(lev+1,l,mid,f,x,v);
  polf(lev+1,mid+1,r,f,x,v);
  cop(_s1,&f[l<<1],mid-l+1);
  times(_s1,&gl[lev+1][(mid+1)<<1],mid-l+1,r-mid+1,r-l+1);
  cop(_s2,&f[(mid+1)<<1],r-mid);
  times(_s2,&gl[lev+1][l<<1],r-mid,mid-l+2,r-l+1);
  for (int i=(l<<1);i<(r<<1|1);i++)f[i]=0;
  for (int i=0;i<r-l+1;i++){
    f[i+(l<<1)]=(_s1[i]+_s2[i])%mod;
    _s1[i]=_s2[i]=0;
  }
}
int n,m;
ll x[Maxn],y[Maxn],y2[Maxn],g[Maxn<<1],f[Maxn<<1];
int main()
{
  preG();
  m=read();
  for (int i=0;i<m;i++)
    {x[i]=read();y[i]=read();}
  qfpre(0,0,m-1,x);
  cop(g,gl[0],m+1);
  dao(g,m+1);
  queryf(0,0,m-1,g,x,y2);
  for (int i=0;i<m;i++)y[i]=y[i]*powM(y2[i])%mod;
  polf(0,0,m-1,f,x,y);
  for (int i=0;i<m;i++)printf("%lld ",f[i]);
  return 0;
}

快速插值非常难写,而且跑得很慢。

很多时候点值可以自己找,为了计算方便,通常找x取值连续的点值(比如说下降幂多项式)

这时的插值虽然不用写求值了,直接连乘搞搞就好,但是两个 \log 还是很慢。