• 51nod_1172_Partial_Sums_V2 拆系数FFT


    题目链接

    定义一个操作:

    [ A_k(i)=sum_{j=0}^{i}A_{k-1}(j) ]

    已知(A_0(x)),求(A_k(x) mod 10^9+7)。A的次数界为(nleq 50,000,kleq 10^9)

    解析:
    (A_k(i))只会对(jgeq i)(A_m(i))产生贡献,不妨假设(A_0(i)=[i=0]),来看一看这个多项式变换k次后的(A_k(i))是多少。

    [egin{align} k=0 & & A= &{1,0,0,0,0,0cdots} onumber\ k=1 & & A= &{1,1,1,1,1,1cdots} onumber\ k=2 & & A= &{1,2,3,4,5,6cdots} onumber\ k=3 & & A= &{1,3,6,10,15,21cdots} onumber\ end{align} ]

    可以发现,(A_0(0))(A_k(j))的贡献倍率为(C_{k+j-1}^{k-1})
    同理,(A_0(i))(A_k(j))的贡献倍率就是(C_{k+(j-i)-1}^{k-1})
    因此:

    [ A_k(j)=sum_{i=0}^{j}{k+(j-i)-1choose k-1}A_0(i)\ 令B(i)={i+k-1choose k-1},则\ A_k(j)=sum_{i=0}^{j}B(j-i)A_0(i)\ ]

    只要构建出B数组,就可以卷积求出(A_k(j))
    (mod 10^9+7)意义下递推B数组。

    [ B(0)=1\ B(i)=frac{B(i-1)*(i+k-1)}{i} ]

    复杂度(O(nlog n))。若是线性预处理1~n的逆元的话复杂度为(O(n))

    [ps:线性求[1,n]在素数p模意义下的逆元:\inv(1)=1,inv(i)=(p-p/i)cdot inv(p mod i) mod p ]

    由于模数较大,直接FFT每一项的值数量级为(10^9*10^9*n=10^{23}),double 和long double 的精度都不足以满足。
    一般有两种姿势。
    第一种是多模数NTT,一般最后合并时要用到高精度中国剩余定理或者一些特殊技巧,而且NTT次数较多,效率较低。
    另一种是拆系数FFT,即将两个多项式A,B拆成四个。

    [ 定义A_p(i)和A_q(i): A(i)=A_p(i)+A_q(i)cdot Base\ 即: A_q(i)=leftlfloorfrac{A(i)}{Base} ight floor,A_p(i)=A(i)mod Base\ B同理\ Base 一般取sqrt{mod}最优 ]

    这样的话

    [egin{align} C(i)= & sum_{j=0}^{i}A(j)B(i-j) onumber\ = & sum_{j=0}^{i}(A_p(j)+Basecdot A_q(j))(B_p(i-j)+Basecdot B_q(i-j)) onumber\ = & sum_{j=0}^{i}A_p(j)B_p(i-j)+Basesum_{j=0}^{i}A_p(j)B_q(i-j) onumber\ & +Basesum_{j=0}^{i}A_q(j)B_p(i-j)+Base^2sum_{j=0}^{i} A_q(j)B_q(i-j) onumber\ end{align} ]

    所以变成了4遍卷积。每一项的阈值数量级为(10^4),卷积结果的阈值数量级为(10^{13})double精度是可以承受的。
    做完卷积之后乘以对应的(Base)再加起来就是(C(i))了。
    复杂度(O(nlog n)),常数略大。

    #include<stdio.h>
    #include<math.h>
    #include<iostream>
    using namespace std;
    typedef long long LL;
    typedef double LDB;
    const int mod=1000000007;
    const LDB PI=2*acos(-1.0);
    inline int Mod1(int x){return x>=mod?x-mod:x;}
    const int Block=31623;
    int n,K,N,L;
    int A[50010],B[50010],Inv[50010];
    struct Complex{
        LDB a,b;
        Complex(LDB x=0,LDB y=0){a=x,b=y;}
        Complex operator + (const Complex& x)const{return Complex(a+x.a,b+x.b);}
        Complex operator * (const Complex& x)const{return Complex(a*x.a-b*x.b,a*x.b+b*x.a);}
        Complex operator - (const Complex& x)const{return Complex(a-x.a,b-x.b);}
        void operator *= (const LDB& x){a*=x,b*=x;}
        void print(){
            printf("(%.3f,%.3f) ",(double)a,(double)b);
        }
    }s1[131082],s2[131082],t1[131082],t2[131082],wn[2][131082];
    int Rev[131082];
    void FFT(Complex* f,int Ty){
        int i,j,l,c,t,k;
        for(i=0;i<N;++i) if(Rev[i]<i) swap(f[Rev[i]],f[i]);
        Complex w,x;
        for(l=2,c=1,t=(N>>1);l<=N;c=l,l<<=1,t>>=1){
            for(i=0;i<N;i+=l){
                for(j=0,k=0;j<c;++j,k+=t){
                    x=f[i+j+c]*wn[Ty][k];
                    f[i+j+c]=f[i+j]-x;
                    f[i+j]=f[i+j]+x;
                }
            }
        }
        if(Ty){
            LDB inv=(LDB)1/N;
            for(int i=0;i<N;++i) f[i]*=inv;
        }
    }
    void Work(){
        for(int i=0;i<n;++i){
            s1[i].a=A[i]/Block;
            s2[i].a=A[i]%Block;
            t1[i].a=B[i]/Block;
            t2[i].a=B[i]%Block;
        }
        for(N=1,L=0;N<n;N<<=1,++L); N<<=1;
        for(int i=0;i<N;++i){
            Rev[i]=((Rev[i>>1]>>1)|((i&1)<<L));
            wn[0][i]=Complex(cos(PI*i/N),sin(PI*i/N));
            wn[1][i]=Complex(cos(-PI*i/N),sin(-PI*i/N));
        }
        FFT(s1,0); FFT(s2,0);
        FFT(t1,0); FFT(t2,0);
        Complex a,b,c,d;
        for(int i=0;i<N;++i){
            a=s1[i],b=s2[i],c=t1[i],d=t2[i];
            s1[i]=a*c;
            s2[i]=b*c;
            t1[i]=a*d;
            t2[i]=b*d;
        }
        FFT(s1,1); FFT(s2,1);
        FFT(t1,1); FFT(t2,1);
        int q,w,e,r,ans;
        for(int i=0;i<n;++i){
            ans=0;
            q=(LL)(s1[i].a+0.5)%mod;
            w=(LL)(s2[i].a+0.5)%mod;
            e=(LL)(t1[i].a+0.5)%mod;
            r=(LL)(t2[i].a+0.5)%mod;
            ans=Mod1(ans+r);
            ans=Mod1(ans+(LL)Block*Mod1(w+e)%mod);
            ans=Mod1(ans+(LL)Block*Block*q%mod);
            printf("%d
    ",ans);
        }
    }   
    int main(){
        scanf("%d%d",&n,&K);
        for(int i=0;i<n;++i) scanf("%d",&A[i]);
        if(K==0){
            for(int i=0;i<n;++i) printf("%d
    ",A[i]);
            return 0;
        }
        --K;
        Inv[1]=1; for(int i=2;i<n;++i) Inv[i]=(LL)(mod-mod/i)*Inv[mod%i]%mod;
        B[0]=1; for(int i=1;i<n;++i) B[i]=(LL)B[i-1]*Inv[i]%mod*(K+i)%mod;
        Work();
        getchar(); getchar();
        return 0;
    }
    
  • 相关阅读:
    TFTP服务器的使用
    SecureCRT Ver 8.1.4 整合汉化绿色版一体包
    SecureCRT Ver 8.1.4 整合汉化绿色版一体包
    深度好贴,mark一下!
    DataSnap——利用TParams进行多表事务更新
    临时表经典使用范例
    PHP对文件的操作方法
    中国银行支付接口(ecshop版)
    常用支付接口实例php版
    PHPCMS网站二次开发配置要点
  • 原文地址:https://www.cnblogs.com/kito/p/7015960.html
走看看 - 开发者的网上家园