【codeforces 623E】dp+FFT+快速幂

题目大意：用$[1,2^k-1]$之间的证书构造一个长度为$n$的序列$a_i$，令$b_i=a_1\ or\ a_2\ or\ ...\ or a_i$，问使得b序列严格递增的方案数，答案对$10^9+7$取模。

数据范围，$n≤10^{18}$，$k≤30000$。

考虑用dp来解决这一题，我们用$f[i][j]$来表示前$i$个数中，使用了$j$个二进制位（注意！并不是前$j$个），那么答案显然为$\sum_{i=0}^{k} \binom{n}{i} \times f[n][i]$。

考虑如何用前面求得的数值来更新$f[x+y][i]$，不妨设$j∈[1,i]$。

不难推出，用了$x$个数，在$i$个二进制位中选用了$j$个二进制位的方案数为$\binom{i}{j} \times f[x][j]$。

然后，用掉$y$个数，并选用余下$i-j$个二进制位的方案数为$f[y][i-j]$。

考虑到前面$x$个数已经选择了$j$个二进制位，那么剩下的$y$个数，在这$j$个位置上，均可以随便填0或1，方案数为$(2^j)^y$。

通过上文分析，得$f[x+y][i]=\sum_{j=1}^{i} f[x][j] \times \binom{i}{j} \times f[y][i-j] \times (2^j)^y$。

通过简单整理，$=i!\sum_{j=1}^{i} \frac{f[x][j]\times (2^j)^y}{j!} \times \frac{f[y][i-j]}{(i-j)!}$。

然后，我们就可以通过NTT，进行dp式子的转移。

不过此题的模数非常恶心，所以需要用任意模数FFT。

考虑到$n$范围非常大，所以$x$和$y$的选择必须要有技巧，我们可以用类似快速幂的算法，加速转移，详情可见代码。

时间复杂度为$O(k\ log\ k\ log\ n)$。

 #include<bits/stdc++.h>

 #define L long long

 #define MOD 1000000007

 #define H 16

 #define M 1<<H

 #define hh 32768

 #define PI acos(-1)

 using namespace std;

 int nn;

 int k; L n;

 L pow_mod(L x,L k){

     L ans=;

     while(k){

         if(k&) ans=ans*x%MOD;

         k>>=; x=x*x%MOD;

     }

     return ans;

 }

 struct cp{

     double i,r;

     cp(){i=r=;}

     cp(double rr,double ii){i=ii; r=rr;}

     friend cp operator +(cp a,cp b){return cp(a.r+b.r,a.i+b.i);}

     friend cp operator -(cp a,cp b){return cp(a.r-b.r,a.i-b.i);}

     friend cp operator *(cp a,cp b){return cp(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}

     L D(){L hhh=(r+0.499); return hhh%MOD;}

 };

 cp w[M][H];

 void init(){

     for(int i=,j=;j<H;j++,i<<=){

         for(int k=;k<i;k++)

         w[k][j]=cp(cos(*PI*k/i),sin(*PI*k/i));

     }

 }

 void change(cp a[],int n){

     for(int i=,j=;i<n-;i++){

         if(i<j) swap(a[i],a[j]);

         int k=n>>;

         while(j>=k) j-=k,k>>=;

         j+=k;

     }

 }

 void FFT(cp a[],int n,int on){

     change(a,n);

     cp wn,u,t;

     for(int h=,i=;h<=n;h<<=,i++){

         for(int j=;j<n;j+=h){

             for(int k=j;k<j+(h>>);k++){

                 wn=w[k-j][i]; if(on==-) wn.i=-wn.i;

                 u=a[k]; t=a[k+(h>>)]*wn;

                 a[k]=u+t; a[k+(h>>)]=u-t;

             }

         }

     }

     if(on==-)

         for(int i=;i<n;i++) a[i].r=a[i].r/n;

 }

 struct poly{

     L a[M];

     poly(){memset(a,,sizeof(a));}

     friend poly operator *(poly a,poly b){

         poly c;

         static cp A[M],B[M],C[M],D[M],E[M],F[M],G[M];

         memset(A,,sizeof(A)); memset(B,,sizeof(B));

         memset(C,,sizeof(C)); memset(D,,sizeof(D));

         for(int i=;i<nn;i++) A[i].r=a.a[i]%hh,B[i].r=a.a[i]/hh;

         for(int i=;i<nn;i++) C[i].r=b.a[i]%hh,D[i].r=b.a[i]/hh;

         FFT(A,nn,); FFT(B,nn,); FFT(C,nn,); FFT(D,nn,);

         for(int i=;i<nn;i++){

             E[i]=A[i]*C[i];

             F[i]=A[i]*D[i]+B[i]*C[i];

             G[i]=B[i]*D[i];

         }

         FFT(E,nn,-); FFT(F,nn,-); FFT(G,nn,-);

         for(int i=;i<nn;i++)

             c.a[i]=(E[i].D()+F[i].D()*hh%MOD+G[i].D()*hh%MOD*hh%MOD)%MOD;

         for(int i=k+;i<nn;i++) c.a[i]=;

         return c;

     }

 };

 L fac[M]={},invfac[M]={};

 L C(int n,int m){return fac[n]*invfac[m]%MOD*invfac[n-m]%MOD;}

 poly ans,f,f1,f2;

 int main(){

     init();

     cin>>n>>k; n--;

     fac[]=; for(int i=;i<=k;i++) fac[i]=fac[i-]*i%MOD;

     invfac[k]=pow_mod(fac[k],MOD-);

     for(int i=k-;~i;i--) invfac[i]=invfac[i+]*(i+)%MOD;

     for(nn=;nn<=(k*);nn<<=);

     L now=;

     for(int i=;i<=k;i++) f.a[i]=;

     ans=f;

     while(n){

         if(n&){

             f1=ans; f2=f;

             for(int i=;i<=k;i++)

             f1.a[i]=f1.a[i]*invfac[i]%MOD*pow_mod(pow_mod(,i),now)%MOD;

             for(int i=;i<=k;i++)

             f2.a[i]=f2.a[i]*invfac[i]%MOD;

             ans=f1*f2;

             for(int i=;i<=k;i++)

             ans.a[i]=ans.a[i]*fac[i]%MOD;

         }

         f1=f; f2=f;

         for(int i=;i<=k;i++)

         f1.a[i]=f1.a[i]*invfac[i]%MOD*pow_mod(pow_mod(,i),now)%MOD;

         for(int i=;i<=k;i++)

         f2.a[i]=f2.a[i]*invfac[i]%MOD;

         f=f1*f2;

         for(int i=;i<=k;i++)

         f.a[i]=f.a[i]*fac[i]%MOD;

         n>>=; now<<=;

     }

     L sum=;

     for(int i=;i<=k;i++)

     sum=(sum+ans.a[i]*C(k,i))%MOD;

     cout<<sum<<endl;

 }

巴特西

【codeforces 623E】dp+FFT+快速幂

最新文章

热门文章