【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)$。
1 #include<bits/stdc++.h> 2 #define L long long 3 #define MOD 1000000007 4 #define H 16 5 #define M 1<<H 6 #define hh 32768 7 #define PI acos(-1) 8 using namespace std; 9 int nn; 10 int k; L n; 11 12 L pow_mod(L x,L k){ 13 L ans=1; 14 while(k){ 15 if(k&1) ans=ans*x%MOD; 16 k>>=1; x=x*x%MOD; 17 } 18 return ans; 19 } 20 21 struct cp{ 22 double i,r; 23 cp(){i=r=0;} 24 cp(double rr,double ii){i=ii; r=rr;} 25 friend cp operator +(cp a,cp b){return cp(a.r+b.r,a.i+b.i);} 26 friend cp operator -(cp a,cp b){return cp(a.r-b.r,a.i-b.i);} 27 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);} 28 L D(){L hhh=(r+0.499); return hhh%MOD;} 29 }; 30 31 cp w[M][H]; 32 void init(){ 33 for(int i=2,j=0;j<H;j++,i<<=1){ 34 for(int k=0;k<i;k++) 35 w[k][j]=cp(cos(2*PI*k/i),sin(2*PI*k/i)); 36 } 37 } 38 39 void change(cp a[],int n){ 40 for(int i=0,j=0;i<n-1;i++){ 41 if(i<j) swap(a[i],a[j]); 42 int k=n>>1; 43 while(j>=k) j-=k,k>>=1; 44 j+=k; 45 } 46 } 47 void FFT(cp a[],int n,int on){ 48 change(a,n); 49 cp wn,u,t; 50 for(int h=2,i=0;h<=n;h<<=1,i++){ 51 for(int j=0;j<n;j+=h){ 52 for(int k=j;k<j+(h>>1);k++){ 53 wn=w[k-j][i]; if(on==-1) wn.i=-wn.i; 54 u=a[k]; t=a[k+(h>>1)]*wn; 55 a[k]=u+t; a[k+(h>>1)]=u-t; 56 } 57 } 58 } 59 if(on==-1) 60 for(int i=0;i<n;i++) a[i].r=a[i].r/n; 61 } 62 63 struct poly{ 64 L a[M]; 65 poly(){memset(a,0,sizeof(a));} 66 friend poly operator *(poly a,poly b){ 67 poly c; 68 static cp A[M],B[M],C[M],D[M],E[M],F[M],G[M]; 69 memset(A,0,sizeof(A)); memset(B,0,sizeof(B)); 70 memset(C,0,sizeof(C)); memset(D,0,sizeof(D)); 71 for(int i=0;i<nn;i++) A[i].r=a.a[i]%hh,B[i].r=a.a[i]/hh; 72 for(int i=0;i<nn;i++) C[i].r=b.a[i]%hh,D[i].r=b.a[i]/hh; 73 FFT(A,nn,1); FFT(B,nn,1); FFT(C,nn,1); FFT(D,nn,1); 74 for(int i=0;i<nn;i++){ 75 E[i]=A[i]*C[i]; 76 F[i]=A[i]*D[i]+B[i]*C[i]; 77 G[i]=B[i]*D[i]; 78 } 79 FFT(E,nn,-1); FFT(F,nn,-1); FFT(G,nn,-1); 80 for(int i=0;i<nn;i++) 81 c.a[i]=(E[i].D()+F[i].D()*hh%MOD+G[i].D()*hh%MOD*hh%MOD)%MOD; 82 for(int i=k+1;i<nn;i++) c.a[i]=0; 83 return c; 84 } 85 }; 86 L fac[M]={0},invfac[M]={0}; 87 L C(int n,int m){return fac[n]*invfac[m]%MOD*invfac[n-m]%MOD;} 88 poly ans,f,f1,f2; 89 int main(){ 90 init(); 91 cin>>n>>k; n--; 92 fac[0]=1; for(int i=1;i<=k;i++) fac[i]=fac[i-1]*i%MOD; 93 invfac[k]=pow_mod(fac[k],MOD-2); 94 for(int i=k-1;~i;i--) invfac[i]=invfac[i+1]*(i+1)%MOD; 95 for(nn=1;nn<=(k*2);nn<<=1); 96 L now=1; 97 for(int i=1;i<=k;i++) f.a[i]=1; 98 ans=f; 99 while(n){ 100 if(n&1){ 101 f1=ans; f2=f; 102 for(int i=1;i<=k;i++) 103 f1.a[i]=f1.a[i]*invfac[i]%MOD*pow_mod(pow_mod(2,i),now)%MOD; 104 for(int i=1;i<=k;i++) 105 f2.a[i]=f2.a[i]*invfac[i]%MOD; 106 ans=f1*f2; 107 for(int i=1;i<=k;i++) 108 ans.a[i]=ans.a[i]*fac[i]%MOD; 109 } 110 f1=f; f2=f; 111 for(int i=1;i<=k;i++) 112 f1.a[i]=f1.a[i]*invfac[i]%MOD*pow_mod(pow_mod(2,i),now)%MOD; 113 for(int i=1;i<=k;i++) 114 f2.a[i]=f2.a[i]*invfac[i]%MOD; 115 f=f1*f2; 116 for(int i=1;i<=k;i++) 117 f.a[i]=f.a[i]*fac[i]%MOD; 118 n>>=1; now<<=1; 119 } 120 L sum=0; 121 for(int i=1;i<=k;i++) 122 sum=(sum+ans.a[i]*C(k,i))%MOD; 123 cout<<sum<<endl; 124 }
