首页 > 代码库 > 51nod 1348 乘积之和

51nod 1348 乘积之和

用(r-l+2)维向量f[l,r]表示区间[l,r]内选i个数(0<=i<=r-l+1)相乘的所有方案之和,可以发现f[l,r]=f[l,m]*f[m+1,r],题目模数100003较小,每次卷积后答案上界大约为1e16,用ntt在两个1e9左右的模数下计算后CRT合并即可,复杂度为O(nlog2n),要注意常数优化

#include<cstdio>#include<cstring>#include<algorithm>typedef long long i64;int _(){    int x=0,f=1,c=getchar();    while(c<48)c==-&&(f=-1),c=getchar();    while(c>47)x=x*10+c-48,c=getchar();    return x*f;}int pow(int a,int n,int p){    int v=1;    for(;n;n>>=1){        if(n&1)v=i64(v)*a%p;        a=i64(a)*a%p;    }    return v;}i64 mul(i64 a,i64 b,i64 p){    i64 s=0;    a%=p;b%=p;    while(b){        if(b&1)(s+=a)%=p;        (a<<=1)%=p;        b>>=1;    }    return s;}int N,X,r[262144];template<const int p,const int g>void ntt(int*a,int t){    for(int i=0;i<N;++i)if(i>r[i])std::swap(a[i],a[r[i]]);    for(int i=1;i<N;i<<=1){        int w=pow(g,(t*(p-1)/(i*2)+p-1),p);        for(int j=0;j<N;j+=i<<1){            int e=1,*b=a+j,*c=b+i;            for(int k=0;k<i;++k,e=i64(e)*w%p){                int x=b[k],y=c[k]*i64(e)%p;                b[k]=(x+y)%p;                c[k]=(x-y+p)%p;            }        }    }    if(t==-1){        i64 I=pow(N,p-2,p);        for(int i=0;i<N;++i)a[i]=a[i]*I%p;    }}int n,q,v0[50007],mem[65536*20+7],*mp=mem,vs[4][65536+7];const int p1=998244353,g1=3,p2=950009857,g2=7;const i64 m1=i64(p1)*pow(p1,p2-2,p2),m2=i64(p2)*pow(p2,p1-2,p1),ps=i64(p1)*p2;int*calc(int L,int R){    int*pos=mp;mp+=R-L+1;    if(L==R){        *pos=v0[L];        return pos;    }    int M=L+R>>1;    int*lp=calc(L,M)-1;    int*rp=calc(M+1,R)-1;    for(N=2,X=0;N<R-L+2;N<<=1,++X);    if(R-L+1<=16){        for(int i=0;i<2;++i)memset(vs[i],0,N*sizeof(int)),vs[i][0]=1;        for(int i=1;i<=M-L+1;++i)vs[0][i]=lp[i];        for(int i=1;i<=R-M;++i)vs[1][i]=rp[i];        for(int i=1;i<=R-L+1;++i){            for(int j=0;j<=i;++j)pos[i-1]=(pos[i-1]+i64(vs[0][j])*vs[1][i-j])%100003;        }        return pos;    }    for(int i=1;i<N;++i)r[i]=r[i>>1]>>1|(i&1)<<X;    for(int i=0;i<4;++i)memset(vs[i],0,N*sizeof(int)),vs[i][0]=1;    for(int i=1;i<=M-L+1;++i)vs[0][i]=vs[1][i]=lp[i];    for(int i=1;i<=R-M;++i)vs[2][i]=vs[3][i]=rp[i];    ntt<p1,g1>(vs[0],1);    ntt<p1,g1>(vs[2],1);    ntt<p2,g2>(vs[1],1);    ntt<p2,g2>(vs[3],1);    for(int i=0;i<N;++i)vs[0][i]=i64(vs[0][i])*vs[2][i]%p1;    for(int i=0;i<N;++i)vs[1][i]=i64(vs[1][i])*vs[3][i]%p2;    ntt<p1,g1>(vs[0],-1);    ntt<p2,g2>(vs[1],-1);    for(int i=1;i<=R-L+1;++i)pos[i-1]=(mul(m2,vs[0][i],ps)+mul(m1,vs[1][i],ps))%ps%100003;    return pos;}int main(){    n=_();q=_();    for(int i=1;i<=n;++i)v0[i]=_();    int*ans=calc(1,n)-1;    while(q--)printf("%d\n",ans[_()]);    return 0;}

 

51nod 1348 乘积之和