多点求值与快速插值
多点求值
给出 $n$ 次多项式 $A(x)$ ,求出 $m$ 个 $x_i$ 对应的 $A(x_i)$
考虑分治,设 $L(x)=\prod_{i=1}^{\frac{n}{2}}(x-x_i)$ , $R(x)=\prod_{i=\frac{n}{2}+1}^n(x-x_i)$
对于 $i \in [1,\frac{n}{2}],F(x_i)=(F \mod L)(x_i)$ , 对于 $i \in (\frac{n}{2},n],F(x_i)=(F \mod R)(x_i)$
就是对于左半部分来说, $F(x)=L(x)q(x)+r(x)$ ,把 $x_i$ 带入后只剩下了 $r(x)$ 的部分,右半部分也同理
于是我们在 $O(nlog^2n)$ 完成了多点求值
快速插值
考虑拉格朗日插值: $F(x)=\sum_{i=1}^n\frac{\prod_{j \ne i}(x-x_j)}{\prod_{j \ne i}(x_i-x_j)}y_i$
如果暴力做的话只能做到 $O(n^2)$ ,考虑怎样优化
下面的部分是个常数,设 $M(x)=\prod_{i=1}^n(x-x_i)$ ,考虑对于每个 $i$ , 下面的式子相当于 $\frac{M(x)}{x-x_i}$ 带入 $x_i$ 之后的值,但是出现了 $\frac{0}{0}$ , 根据洛必达法则,它的值相当于 $M'(x)$ 带入 $x_i$ 后的值,于是可以分治+ $Ntt$ 后求导,再多点求值求出每个 $x_i$ 对应的值
于是 $F(x)=\sum_{i=1}^n\prod_{j \ne i}(x-x_j)v_i$ , 这一部分就可以递归处理
设 $L(x)=\prod_{i=1}^{\frac{n}{2}}(x-x_i)$ , $R(x)=\prod_{i=\frac{n}{2}+1}^n(x-x_i)$,
于是 $F(x)=R(x)\sum_{i=1}^{\frac{n}{2}}v_i\prod_{j=1}^{\frac{n}{2}}[j \ne i](x-x_j)+L(x)\sum_{i=\frac{n}{2}+1}^{n}v_i\prod_{j=\frac{n}{2}+1}^{n}[j \ne i](x-x_j)$
于是我们在 $O(nlog^2n)$ 完成了快速插值
代码
(代码本机跑是挂的但交洛谷是OK的不知道为什么,哪位好心人帮我解答感谢您嘞)
#include <bits/stdc++.h> #define _(d) while(d(isdigit(c=getchar()))) using namespace std; int Rd(){ char c;_(!);int x=c^48;_()x=(x<<3)+(x<<1)+(c^48);return x; } int vv[20]; void Wr(int x){ if (!x){putchar(48);return;} int v=0;while(x) vv[++v]=x%10,x/=10; for (int i=v;i;i--) putchar(vv[i]^48); } const int N=4e5+5,P=998244353; int A[N],B[N],t,p,re[N],C[N],D[N],E[N],G[2][19][N],c[N],*d[N],e[N],*f[N],g[N],*h[N]; int X(int x){return x>=P?x-P:x;} int K(int x,int y){ int z=1; for (;y;y>>=1,x=1ll*x*x%P) if (y&1) z=1ll*z*x%P; return z; } void init(){ for (int i=0;i<19;i++){ G[0][i][0]=G[1][i][0]=1; G[0][i][1]=K(3,(P-1)/(1<<(i+1))); G[1][i][1]=K((P+1)/3,(P-1)/(1<<(i+1))); for (int j=2;j<(1<<i);j++) G[0][i][j]=1ll*G[0][i][j-1]*G[0][i][1]%P, G[1][i][j]=1ll*G[1][i][j-1]*G[1][i][1]%P; } } void pre(int l){ for (t=1,p=0;t<l;t<<=1,p++); for (int i=0;i<t;i++) re[i]=(re[i>>1]>>1)|((i&1)<<(p-1)); } void Ntt(int *a,int o){ for (int i=0;i<t;i++) if (i<re[i]) swap(a[i],a[re[i]]); for (int v=0,i=1;i<t;i<<=1,v++){ for (int x,y,j=0;j<t;j+=(i<<1)) for (int k=0;k<i;k++) x=a[j+k],y=1ll*G[o][v][k]*a[i+j+k]%P, a[j+k]=X(x+y),a[i+j+k]=X(x-y+P); } if (o) for (int i=0,v=K(t,P-2);i<t;i++) a[i]=1ll*a[i]*v%P; } void dao(int *a,int l){ for (int i=1;i<l;i++) a[i-1]=1ll*i*a[i]%P; a[l-1]=0; } void inv(int *a,int *b,int l){ if (l==1){ b[0]=K(a[0],P-2); return; } inv(a,b,(l+1)>>1); for (int i=0;i<l;i++) A[i]=a[i],B[i]=b[i]; pre(l<<1);Ntt(A,0);Ntt(B,0); for (int i=0;i<t;i++) A[i]=1ll*A[i]*B[i]%P*B[i]%P; Ntt(A,1); for (int i=0;i<l;i++) b[i]=X(X(b[i]<<1)+P-A[i]); for (int i=0;i<t;i++) A[i]=B[i]=0; } void dvs(int *f,int *g,int *q,int *d,int n,int m){ reverse(f,f+n+1);reverse(g,g+m+1); for (int i=min(n-m,m);~i;i--) E[i]=g[i]; inv(E,q,n-m+1); for (int i=0;i<=n-m;i++) A[i]=q[i],B[i]=f[i]; pre((n-m+1)<<1);Ntt(A,0);Ntt(B,0); for (int i=0;i<t;i++) A[i]=1ll*A[i]*B[i]%P; Ntt(A,1); for (int i=0;i<=n-m;i++) q[i]=A[i]; for (int i=0;i<t;i++) A[i]=B[i]=0; reverse(q,q+n-m+1); reverse(f,f+n+1);reverse(g,g+m+1); for (int i=0;i<=m;i++) A[i]=g[i]; for (int i=0;i<=n-m;i++) B[i]=q[i]; pre(n+2);Ntt(A,0);Ntt(B,0); for (int i=0;i<t;i++) A[i]=1ll*A[i]*B[i]%P; Ntt(A,1); for (int i=0;i<m;i++) d[i]=X(f[i]-A[i]+P); for (int i=0;i<t;i++) A[i]=B[i]=E[i]=0; } #define Ls k<<1 #define Rs k<<1|1 #define mid ((l+r)>>1) void build(int k,int l,int r){ if (l==r){ c[k]=1;d[k]=new int[2]; d[k][0]=X(P-C[l]); d[k][1]=1; return; } build(Ls,l,mid);build(Rs,mid+1,r); c[k]=c[Ls]+c[Rs]; d[k]=new int[c[k]+1]; for (int i=0;i<=c[Ls];i++) A[i]=d[Ls][i]; for (int i=0;i<=c[Rs];i++) B[i]=d[Rs][i]; pre(c[k]+2);Ntt(A,0);Ntt(B,0); for (int i=0;i<t;i++) A[i]=1ll*A[i]*B[i]%P; Ntt(A,1); for (int i=0;i<=c[k];i++) d[k][i]=A[i]; for (int i=0;i<t;i++) A[i]=B[i]=0; } void push(int k,int l){ e[k]=l;f[k]=new int[l+1]; } void qry(int k,int l,int r){ if (l==r){D[l]=f[k][0];return;} push(Ls,c[Ls]-1);push(Rs,c[Rs]-1); dvs(f[k],d[Ls],C,f[Ls],e[k],c[Ls]); for (int i=0;i<=e[k]-c[Ls];i++) C[i]=0; dvs(f[k],d[Rs],C,f[Rs],e[k],c[Rs]); for (int i=0;i<=e[k]-c[Rs];i++) C[i]=0; qry(Ls,l,mid);qry(Rs,mid+1,r); } void multip(int *a,int n,int *b,int m){ push(1,m-1); for (int i=0;i<m;i++) f[1][i]=a[i]; qry(1,1,m); } void work(int k,int x,int y,int l,int r){ for (int i=0;i<=g[x];i++) A[i]=h[x][i]; for (int i=0;i<=c[y];i++) B[i]=d[y][i]; pre(r-l+2);Ntt(A,0);Ntt(B,0); for (int i=0;i<t;i++) A[i]=1ll*A[i]*B[i]%P; Ntt(A,1); for (int i=0;i<=g[k];i++) h[k][i]=X(h[k][i]+A[i]); for (int i=0;i<t;i++) A[i]=B[i]=0; } void calcf(int k,int l,int r){ if (l==r){ h[k]=new int[1]; h[k][0]=D[l]; return; } calcf(Ls,l,mid);calcf(Rs,mid+1,r); g[k]=r-l;h[k]=new int[r-l+1]; work(k,Ls,Rs,l,r);work(k,Rs,Ls,l,r); } void insv(int *a,int *b,int n){ for (int i=1;i<=n;i++) C[i]=a[i]; build(1,1,n); dao(d[1],n+1); for (int i=0;i<=n;i++) C[i]=0; multip(d[1],n-1,a,n); for (int i=1;i<=n;i++) D[i]=1ll*b[i]*K(D[i],P-2)%P; calcf(1,1,n); for (int i=0;i<n;i++) Wr(h[1][i]),putchar(i<n-1?' ':'\n'); } int n,a[N],b[N]; int main(){ init();n=Rd(); for (int i=1;i<=n;i++) a[i]=Rd(),b[i]=Rd(); insv(a,b,n); return 0; }