LG5050 多项式多点求值 和 LG5158 多项式快速插值
LG5050 多项式多点求值
给定一个 \(n\) 次多项式 \(f(x)\) ，现在请你对于 \(i \in [1,m]\) ，求出 \(f(a_i)\)。
答案对 \(998244353\) 取模。
\(n,m∈[1,64000]\)
题解
我们将要求值的点均分成两份，构造多项式\(P_0(x)=\prod\limits_{i=1}^{\lfloor\frac n 2\rfloor}(x-x_i)\)，\(P_1(x)=\prod\limits_{i=\lfloor\frac n 2\rfloor+1}^{n}(x-x_i)\)。
显然，对于\(\forall i\in[1,\lfloor\frac n 2 \rfloor]\)，有\(P_0(x_i)=0\)。\(P_1\)同理。
我们假设多项式\(D(x),R(x)\)满足：\(D(x)\)是一个\(n-\lfloor\frac n 2\rfloor\)次多项式，\(R(x)\)是一个次数不超过\(\lfloor\frac n 2\rfloor-1\)的多项式，且\(A(x)=P_0(x)D(x)+R(x)\)。
那么对于\(\forall i\in[1,\lfloor\frac n 2 \rfloor]\)，有\(A(x_i)=R(x_i)\)。\(P_1\)同理可得。
由于\(R(x)\)的次数是\(A(x)\)的次数的一半，所以我们把一个规模为\(n\)的问题，用\(O(n\log n)\)的复杂度（多项式取模，详见多项式除法），转化为两个规模为\(\frac n 2\)的问题。
分治计算即可，由主定理得时间复杂度\(O(n\log^2 n)\)。
求每次的\(P_0(x),P_1(x)\)，可以开始用一次分治NTT预处理。此处时间复杂度也为\(O(n\log^2 n)\)。
typedef vector<int> polynomial;
void num_trans(polynomial&a,int dir){
int lim=a.size();
static vector<int> rev,w[2];
if(rev.size()!=lim){
rev.resize(lim);
int len=log2(lim);
for(int i=0;i<lim;++i) rev[i]=rev[i>>1]>>1|(i&1)<<(len-1);
for(int dir=0;dir<2;++dir){
w[dir].resize(lim);
w[dir][0]=1,w[dir][1]=fpow(g[dir],(mod-1)/lim);
for(int i=2;i<lim;++i) w[dir][i]=mul(w[dir][i-1],w[dir][1]);
}
}
for(int i=0;i<lim;++i)if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int step=1;step<lim;step<<=1){
int quot=lim/(step<<1);
for(int i=0;i<lim;i+=step<<1){
int j=i+step;
for(int k=0;k<step;++k){
int t=mul(w[dir][quot*k],a[j+k]);
a[j+k]=add(a[i+k],mod-t),a[i+k]=add(a[i+k],t);
}
}
}
if(dir){
int ilim=fpow(lim,mod-2);
for(int i=0;i<lim;++i) a[i]=mul(a[i],ilim);
}
}
polynomial poly_inv(polynomial a,int n){
polynomial b(1,fpow(a[0],mod-2));
if(n==1) return b;
int lim=2;
for(;lim<n;lim<<=1){
polynomial a1(a.begin(),a.begin()+lim);
a1.resize(lim<<1),num_trans(a1,0);
b.resize(lim<<1),num_trans(b,0);
for(int i=0;i<lim<<1;++i) b[i]=mul(add(2,mod-mul(a1[i],b[i])),b[i]);
num_trans(b,1),b.resize(lim);
}
a.resize(lim<<1),num_trans(a,0);
b.resize(lim<<1),num_trans(b,0);
for(int i=0;i<lim<<1;++i) b[i]=mul(add(2,mod-mul(a[i],b[i])),b[i]);
num_trans(b,1),b.resize(n);
return b;
}
polynomial poly_div(polynomial f,polynomial g){
int n=f.size()-1,m=g.size()-1;
reverse(g.begin(),g.end()),g.resize(n-m+1),g=poly_inv(g,n-m+1);
reverse(f.begin(),f.end()),f.resize(n-m+1);
int lim=1<<int(ceil(log2(((n-m)<<1)+1)));
f.resize(lim),num_trans(f,0);
g.resize(lim),num_trans(g,0);
for(int i=0;i<lim;++i) f[i]=mul(f[i],g[i]);
num_trans(f,1),f.resize(n-m+1);
return reverse(f.begin(),f.end()),f;
}
polynomial poly_mod(polynomial f,polynomial g){
int n=f.size()-1,m=g.size()-1;
polynomial q=poly_div(f,g);
int lim=1<<int(ceil(log2(n+1)));
q.resize(lim),num_trans(q,0);
g.resize(lim),num_trans(g,0);
for(int i=0;i<lim;++i) q[i]=mul(q[i],g[i]);
num_trans(q,1);
for(int i=0;i<m;++i) f[i]=add(f[i],mod-q[i]);
return f.resize(m),f;
}
co int N=64000+1;
int a[N];
polynomial s[N<<2];
#define lc (x<<1)
#define rc (x<<1|1)
polynomial poly_mul(polynomial a,polynomial b){
int n=a.size()-1,m=b.size()-1;
int lim=1<<int(ceil(log2(n+m+1)));
a.resize(lim),num_trans(a,0);
b.resize(lim),num_trans(b,0);
for(int i=0;i<lim;++i) a[i]=mul(a[i],b[i]);
num_trans(a,1),a.resize(n+m+1);
return a;
}
void build(int x,int l,int r){
if(l==r){
s[x].resize(2);
s[x][0]=mod-a[l],s[x][1]=1;
return;
}
int mid=(l+r)>>1;
build(lc,l,mid),build(rc,mid+1,r);
s[x]=poly_mul(s[lc],s[rc]);
}
void solve(int x,int l,int r,polynomial f){
if(l==r){
printf("%d\n",f[0]);
return;
}
int mid=(l+r)>>1;
solve(lc,l,mid,poly_mod(f,s[lc])),solve(rc,mid+1,r,poly_mod(f,s[rc]));
}
int main(){
int n=read<int>(),m=read<int>();
if(!m) return 0;
polynomial f(n+1);
for(int i=0;i<=n;++i) read(f[i]);
for(int i=1;i<=m;++i) read(a[i]);
build(1,1,m);
if(n>=m) f=poly_mod(f,s[1]);
solve(1,1,m,f);
return 0;
}
卡常数做法
定义减法卷积\(\text{mul}(F(x),G(x))=\sum_i\sum_{j\leq i}f_ig_jx^{i-j}\)
注意到\(F(a_i)=[x^0]\text{mul}(F(x),\frac{1}{1-a_ix})\)。
减法卷积有性质\(\text{mul}(F(x),G(x)H(x))=\text{mul}(\text{mul}(F(x),G(x)),H(x))\)，证明考虑\(i-(j+k)=i-j-k\)。
于是类似分治，根节点传入\(F_{1,m}(x)=\text{mul}(F(x),\frac{1}{\prod_{i=1}^m1-a_ix})\)
从\([l,r]\)递归的时候，\([l,\text{mid}]\)传入\(F_{l,\text{mid}}(x)=\text{mul}(F_{l,r}(x),\prod_{i=\text{mid}+1}^r(1-a_ix))\)，\([\text{mid}+1,r]\)传入\(F_{\text{mid+1,r}(x)}=\text{mul}(F_{l,r}(x),\prod_{i=l}^{\text{mid}}(1-a_ix))\)。
区间\([l,r]\)的多项式只用保留到\(x^{r-l+1}\)，因为之后的项不可能通过减法卷积贡献到\([x^0]\)。
时间复杂度\(O(n\log^2 n)\)，因为没有多项式取模所以常数小。
CO int N=1<<17;
int omg[2][N],rev[N];
void NTT(poly&a,int dir){
int lim=a.size(),len=log2(lim);
for(int i=0;i<lim;++i){
int r=rev[i]>>(17-len);
if(i<r) swap(a[i],a[r]);
}
for(int i=1;i<lim;i<<=1)
for(int j=0;j<lim;j+=i<<1)for(int k=0;k<i;++k){
int t=mul(omg[dir][N/(i<<1)*k],a[j+i+k]);
a[j+i+k]=add(a[j+k],mod-t),a[j+k]=add(a[j+k],t);
}
if(dir){
int ilim=fpow(lim,mod-2);
for(int i=0;i<lim;++i) a[i]=mul(a[i],ilim);
}
}
poly operator~(poly a){
int n=a.size();
poly b={fpow(a[0],mod-2)};
if(n==1) return b;
a.resize(1<<(int)ceil(log2(n)));
for(int lim=2;lim<2*n;lim<<=1){
poly c(a.begin(),a.begin()+lim);
c.resize(lim<<1),NTT(c,0);
b.resize(lim<<1),NTT(b,0);
for(int i=0;i<lim<<1;++i) b[i]=mul(2+mod-mul(c[i],b[i]),b[i]);
NTT(b,1),b.resize(lim);
}
return b.resize(n),b;
}
poly operator*(poly a,poly b){
if(a.size()<=20 or b.size()<=20){
int n=a.size()-1,m=b.size()-1;
a.resize(n+m+1);
for(int i=n;i>=0;--i){
for(int j=m;j>=1;--j) a[i+j]=add(a[i+j],mul(a[i],b[j]));
a[i]=mul(a[i],b[0]);
}
return a;
}
int n=a.size()+b.size()-1,lim=1<<(int)ceil(log2(n));
a.resize(lim),NTT(a,0);
b.resize(lim),NTT(b,0);
for(int i=0;i<lim;++i) a[i]=mul(a[i],b[i]);
NTT(a,1),a.resize(n);
return a;
}
poly operator/(poly a,poly b){
if(a.size()<=20 or b.size()<=20){
int n=a.size()-1,m=b.size()-1;
for(int i=0;i<=n;++i){
for(int j=min(i,m);j>=1;--j) a[i-j]=add(a[i-j],mul(a[i],b[j]));
a[i]=mul(a[i],b[0]);
}
return a;
}
int n=a.size()-1,m=b.size()-1;
if(n<m) b.resize(n+1),m=n;
reverse(b.begin(),b.end());
int lim=1<<(int)ceil(log2(n+m+1));
a.resize(lim),NTT(a,0);
b.resize(lim),NTT(b,0);
for(int i=0;i<lim;++i) a[i]=mul(a[i],b[i]);
NTT(a,1);
for(int i=0;i<=n;++i) a[i]=a[i+m];
return a.resize(n+1),a;
}
int a[N];
poly g[N];
#define lc (x<<1)
#define rc (x<<1|1)
#define mid ((l+r)>>1)
void build(int x,int l,int r){
if(l==r) {g[x]={1,mod-read<int>()}; return;}
build(lc,l,mid),build(rc,mid+1,r);
g[x]=g[lc]*g[rc];
}
void solve(int x,int l,int r,CO poly&f){
if(l==r) {printf("%d\n",f[0]); return;}
poly lf=f/g[rc];lf.resize(mid-l+2);
solve(lc,l,mid,lf);
poly rf=f/g[lc];rf.resize(r-mid+1);
solve(rc,mid+1,r,rf);
}
#undef lc
#undef rc
#undef mid
int main(){
omg[0][0]=1,omg[0][1]=fpow(3,(mod-1)/N);
omg[1][0]=1,omg[1][1]=fpow(omg[0][1],mod-2);
rev[0]=0,rev[1]=1<<16;
for(int i=2;i<N;++i){
omg[0][i]=mul(omg[0][i-1],omg[0][1]);
omg[1][i]=mul(omg[1][i-1],omg[1][1]);
rev[i]=rev[i>>1]>>1|(i&1)<<16;
}
int n=read<int>(),m=read<int>();
poly f(n+1);
for(int i=0;i<=n;++i) read(f[i]);
build(1,1,m);
g[1].resize(n+1),f=f/~g[1],f.resize(m+1);
solve(1,1,m,f);
return 0;
}
LG5158 多项式快速插值
给出 \(n\) 个点 \((X_i,Y_i)\)。求一个 \(n-1\) 次的多项式 \(f(x)\)，使得 \(f(X_i)= Y_i\mod{998244353}\)。
\(1⩽n⩽100000\)
题解
思路是优化拉格朗日插值。先看看原式
\[f(x)=\sum_{i = 1}^{n} y_i \prod_{i \not = j} \frac{x - x_j}{x_i - x_j}
\]
分离常数
\[f(x)=\sum_{i = 1}^{n}{ y_i\over \prod_{i \not = j}{x_i - x_j}} \prod_{i \not = j}(x - x_j)
\]
来考虑一下\({ y_i\over \prod_{i \not = j}{x_i - x_j}}\)这个东西，上面是个常数，那么只考虑下面。如果我们设多项式\(g(x)=\prod_{i=1}^n (x-x_i)\)，那么下面那个东西就是\({g(x_i)\over (x-x_i)}\)
这分子分母全为\(0\)怎么求？根据洛必达法则，如果
\[\lim_{x\to a}f(x)=0,\lim_{x\to a}g(x)=0
\]
则有
\[\lim_{x\to a}{f(x)\over g(x)}=\lim_{x\to a}{f'(x)\over g'(x)}
\]
那么我们代入之后可以发现\({g(x_i)\over (x-x_i)}=g'(x_i)\)
先分治NTT算出\(g\)，然后对\(g'\)多点求值把每个点处的值算出来就好了
那么接下来我们就考虑分治，设\(f_{l,r}\)表示\((x_l,y_l)\sim(x_r,y_r)\)这些点算出来的答案，则有
\[\begin{aligned}
f_{l,r}&=\sum_{i = l}^{r}{ y_i\over g'(x_i)} \prod_{j=l,i \not = j}^r(x - x_j)\\&=\prod_{j=mid+1}^r(x - x_j)\sum_{i = l}^{mid}{ y_i\over g'(x_i)} \prod_{j=l,i \not = j}^{mid}(x - x_j)+\prod_{j=l}^{mid}(x - x_j)\sum_{i = mid+1}^{r}{ y_i\over g'(x_i)} \prod_{j=mid+1,i \not = j}^{r}(x - x_j)\\&=\prod_{j=mid+1}^r(x - x_j)f_{l,mid}+\prod_{j=l}^{mid}(x - x_j)f_{mid+1,r}\\
\end{aligned}
\]
时间复杂度为\(O(n\log^2n)\)。这题卡常，需要小范围多项式乘法暴力才能过。
typedef vector<int> polynomial;
void num_trans(polynomial&a,int dir){
int lim=a.size();
static vector<int> rev,w[2];
if(rev.size()!=lim){
rev.resize(lim);
int len=log2(lim);
for(int i=0;i<lim;++i) rev[i]=rev[i>>1]>>1|(i&1)<<(len-1);
for(int dir=0;dir<2;++dir){
static co int g[2]={3,332748118};
w[dir].resize(lim);
w[dir][0]=1,w[dir][1]=fpow(g[dir],(mod-1)/lim);
for(int i=2;i<lim;++i) w[dir][i]=mul(w[dir][i-1],w[dir][1]);
}
}
for(int i=0;i<lim;++i)if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int step=1;step<lim;step<<=1){
int quot=lim/(step<<1);
for(int i=0;i<lim;i+=step<<1){
int j=i+step;
for(int k=0;k<step;++k){
int t=mul(w[dir][quot*k],a[j+k]);
a[j+k]=add(a[i+k],mod-t),a[i+k]=add(a[i+k],t);
}
}
}
if(dir){
int ilim=fpow(lim,mod-2);
for(int i=0;i<lim;++i) a[i]=mul(a[i],ilim);
}
}
polynomial poly_inv(polynomial a,int n){
polynomial b(1,fpow(a[0],mod-2));
if(n==1) return b;
int lim=2;
for(;lim<n;lim<<=1){
polynomial a1(a.begin(),a.begin()+lim);
a1.resize(lim<<1),num_trans(a1,0);
b.resize(lim<<1),num_trans(b,0);
for(int i=0;i<lim<<1;++i) b[i]=mul(add(2,mod-mul(a1[i],b[i])),b[i]);
num_trans(b,1),b.resize(lim);
}
a.resize(lim<<1),num_trans(a,0);
b.resize(lim<<1),num_trans(b,0);
for(int i=0;i<lim<<1;++i) b[i]=mul(add(2,mod-mul(a[i],b[i])),b[i]);
num_trans(b,1),b.resize(n);
return b;
}
polynomial poly_div(polynomial f,polynomial g){
int n=f.size()-1,m=g.size()-1;
reverse(g.begin(),g.end()),g.resize(n-m+1),g=poly_inv(g,n-m+1);
reverse(f.begin(),f.end()),f.resize(n-m+1);
int lim=1<<int(ceil(log2(((n-m)<<1)+1)));
f.resize(lim),num_trans(f,0);
g.resize(lim),num_trans(g,0);
for(int i=0;i<lim;++i) f[i]=mul(f[i],g[i]);
num_trans(f,1),f.resize(n-m+1);
return reverse(f.begin(),f.end()),f;
}
polynomial poly_mod(polynomial f,polynomial g){
int n=f.size()-1,m=g.size()-1;
polynomial q=poly_div(f,g);
int lim=1<<int(ceil(log2(n+1)));
q.resize(lim),num_trans(q,0);
g.resize(lim),num_trans(g,0);
for(int i=0;i<lim;++i) q[i]=mul(q[i],g[i]);
num_trans(q,1);
for(int i=0;i<m;++i) f[i]=add(f[i],mod-q[i]);
return f.resize(m),f;
}
polynomial poly_mul(polynomial a,polynomial b){
int n=a.size()-1,m=b.size()-1;
if((LL)n*m<=5000){ // edit 1
polynomial c(n+m+1);
for(int i=0;i<=n;++i)
for(int j=0;j<=m;++j) c[i+j]=add(c[i+j],mul(a[i],b[j]));
return c;
}
int lim=1<<int(ceil(log2(n+m+1)));
a.resize(lim),num_trans(a,0);
b.resize(lim),num_trans(b,0);
for(int i=0;i<lim;++i) a[i]=mul(a[i],b[i]);
num_trans(a,1),a.resize(n+m+1);
return a;
}
polynomial poly_add(polynomial a,polynomial b){
if(a.size()<b.size()) swap(a,b);
for(int i=0;i<b.size();++i) a[i]=add(a[i],b[i]);
return a;
}
co int N=100000+1;
#define lc (x<<1)
#define rc (x<<1|1)
namespace eval{
int n,m,a[N],ans[N];
polynomial f,s[N<<2];
void build(int x,int l,int r){
if(l==r){
s[x].resize(2);
s[x][0]=mod-a[l],s[x][1]=1;
return;
}
int mid=(l+r)>>1;
build(lc,l,mid),build(rc,mid+1,r);
s[x]=poly_mul(s[lc],s[rc]);
}
void solve(int x,int l,int r,polynomial f){
if(l==r){
ans[l]=f[0];
return;
}
int mid=(l+r)>>1;
solve(lc,l,mid,poly_mod(f,s[lc])),solve(rc,mid+1,r,poly_mod(f,s[rc]));
}
void main(){
build(1,1,m);
if(n>=m) f=poly_mod(f,s[1]);
solve(1,1,m,f);
}
}
int X[N],Y[N];
polynomial g[N<<2],f[N<<2];
void build(int x,int l,int r){
if(l==r){
g[x].resize(2);
g[x][0]=mod-X[l],g[x][1]=1;
return;
}
int mid=(l+r)>>1;
build(lc,l,mid),build(rc,mid+1,r);
g[x]=poly_mul(g[lc],g[rc]);
}
void solve(int x,int l,int r){
if(l==r){
f[x].assign(1,mul(Y[l],fpow(eval::ans[l],mod-2)));
return;
}
int mid=(l+r)>>1;
solve(lc,l,mid),solve(rc,mid+1,r);
f[x]=poly_add(poly_mul(f[lc],g[rc]),poly_mul(f[rc],g[lc]));
}
int main(){
// freopen("LG5158.in","r",stdin),freopen("LG5158.out","w",stdout);
int n=read<int>();
for(int i=1;i<=n;++i) read(X[i]),read(Y[i]);
build(1,1,n);
eval::n=n-1,eval::f.resize(n);
for(int i=0;i<=eval::n;++i) eval::f[i]=mul(g[1][i+1],i+1);
eval::m=n;
for(int i=1;i<=eval::m;++i) eval::a[i]=X[i];
eval::main();
solve(1,1,n);
for(int i=0;i<n;++i) printf("%d ",f[1][i]);
return 0;
}
静渊以有谋，疏通而知事。
