zoukankan      html  css  js  c++  java
  • 多点求值与快速插值

    多点求值

    给出 $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;
    }
  • 相关阅读:
    Zabbix 3.2.1 安装 Graphtree3.0.4 或 Graphtree3.2.x
    jquery-1
    AngularJS (1)
    css-2 (Text Font)
    css
    Java经验
    js经验
    mysql经验
    MySQL 函数
    jquery 效果
  • 原文地址:https://www.cnblogs.com/xjqxjq/p/12241568.html
Copyright © 2011-2022 走看看