多点求值与快速插值

Posted 2020-11-25 xjqxjq

多点求值

给出 $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}^nfrac{prod_{j e i}(x-x_j)}{prod_{j e 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}^nprod_{j e 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_iprod_{j=1}^{frac{n}{2}}[j e i](x-x_j)+L(x)sum_{i=frac{n}{2}+1}^{n}v_iprod_{j=frac{n}{2}+1}^{n}[j e 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?‘ ‘:‘
‘);
}
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;
}