从拉插到快速插值求值
command_block
2019-07-25 16:05:03
~~标题原来很长的,简写却之后变得很拗口~~
老是忘记这个东西,干脆写个文章好了……
为了显示点值和未知数的区别,自由元 $X$ **一律大写**。
**广告** : [多项式计数杂谈](https://www.luogu.com.cn/blog/command-block/sheng-cheng-han-shuo-za-tan)
## 1.拉格朗日插值
[P4781 【模板】拉格朗日插值](https://www.luogu.org/problem/P4781)
- **题意** :
有$n$个不同的点值$(x_i,y_i)$,让你求出一个低于$n$次的多项式,满足经过上述$n$个点。
首先由代数基本定理, $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)$ ,那么得到:
$T_k(x_k)=\prod\limits_{i=1}^n[k≠i](x_k-x_i)$ ,这是个常数。
我们要修正这个多项式,只需要令其除以这个常数,在 $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](https://www.luogu.org/problem/CF622F)
这是一道很好的例题,题意不再赘述。
容易感性认知,自然数 $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:**
```cpp
#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)$。
```cpp
#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与多项式全家桶](https://www.luogu.org/blog/command-block/ntt-yu-duo-xiang-shi-quan-jia-tong) ,Ln+Exp可以先忽略。
[P5050 【模板】多项式多点求值](https://www.luogu.org/problem/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)$
$F_2$ 是膜了一个**一次**多项式得到的,肯定只有一个常数,**这个常数就是答案**。
但是对于每个点值都暴力取模岂不是要 $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一定要弄个常数小的来。
```cpp
#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 【模板】快速阶乘算法](https://www.luogu.org/problem/P5282)
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
------------
[P5158 【模板】多项式快速插值](https://www.luogu.org/problem/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)$
- 首先算每个$\prod\limits_{i=1,k≠i}^m(x_k-x_i)$
我们学了那么多多项式技巧,不能被这么一个常数骗过去啊!
我们设$G(X)=\prod\limits_{i=1}^m(X-x_i)$,这可以分治NTT求出来。
$\lim\limits_{X→x_k}\dfrac{G}{(X-x_k)}=\prod\limits_{i=1,k≠i}^m(x_k-x_i)$就是我们要求的东西。
这个式子就是$\dfrac{0}{0}$形式的**洛必达**,由$\lim\limits_{X→x_k}G(X)=\lim\limits_{X→x_k}(X-x_k)=0$
那么$\lim\limits_{X→x_k}\dfrac{G(X)}{(X-x_k)}=\lim\limits_{X→x_k}\dfrac{G'(X)}{(X-x_k)'}$(上下同时求导)
$=\lim\limits_{X→x_k}G'(x)=G'(x_k)$
那么我们写个多项式多点求值,就能从$G'(x)$里弄出上述插值常数了。
现在就能得到$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)$
$=G_{[mid+1,r]}F_{[l,mid]}+G_{[l,mid]}F_{[mid+1,r]}$(注意是交叉乘哦)
边界就是$F_{[k,k]}=V_k=\dfrac{y_k}{G'(x_k)}$
那么就变成两个规模为一半的子问题了,分治NTT上!
复杂度$O(nlog^2n)$,常数极大。
如何减小常数呢?$G'$,多点求值,分治NTT中的$G$都是同一个$G$,那么可以只算一次。
此外,分治范围小到一定程度的时候可以暴力拉插公式(这个优化我没写,因为其实**求值才是瓶颈**)。
卡到最后当然是预处理单位根啦。
```cpp
#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$ 还是很慢。