CodeForces 623 E.Transforming Sequence（dp+FFT+快速幂）

Description

用$[1,2^k-1]$之间的整数构造一个长度为$n$的序列$a_1,...,a_n$，令$b_i=a_1|a_2|...|a_i,1\le i\le n$，问使得$b$序列严格递增的方案数

Input

两个整数$n,k(1\le n\le 10^{18},1\le k\le 30000)$

Output

输出满足条件的$a$序列个数，结果模$10^9+7$

Sample Input

1 2

Sample Output

3

Solution

为保证$b$序列严格递增，每次$a_i$需要至少在一个$a_1$~$a_{i-1}$都是$0$的位上是$1$，故序列长度至多为$k$，即如$n>k$则无解，下面只考虑$n\le k$的情况

设$dp[i][j]$为前$i$个数占据$j$个二进制位的方案数，不考虑这些位在$k$个位中的位置，则答案应为$\sum\limits_{i=1}^kC_k^i\cdot dp[n][i]$

假设要求$dp[x+y][i]$，假设前$x$个数已经占据了$i$位中的$j$位，方案数$C_i^j\cdot dp[x][j]$，而后$y$个数要把剩下的$i-j$位占据，而在前$x$个数占据的$j$位上可以随意取值，方案数$dp[y][i-j]\cdot 2^{yj}$，故有$dp[x+y][i]=\sum\limits_{j=0}^idp[x][j]\cdot dp[y][i-j]\cdot 2^{yj}\cdot C_i^j$

在求$dp[n][i]$时，如果$n$是偶数，那么有$\frac{dp[n][i]}{i!}=\sum\limits_{j=0}^i\frac{dp[\frac{n}{2}][j]\cdot 2^{\frac{n}{2}j}}{j!}\cdot \frac{dp[\frac{n}{2}][i-j]}{(i-j)!}$

如果$n$是奇数，可以先用上面的方法求出$dp[n-1]$的值，然后$\frac{dp[n][i]}{i!}=\sum\limits_{j=0}^i\frac{dp[n-1][j]\cdot 2^j}{j!}\cdot \frac{1}{(i-j)!}$

故可以快速幂求$dp[n]$，每次两个序列卷积得到一个新的序列，$FFT$即可，时间复杂度$O(klog_2^2k)$

Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
#define maxfft 65536+5
#define mod 1000000007
const double pi=acos(-1.0);
struct cp
{
    double a,b;
    cp operator +(const cp &o)const {return (cp){a+o.a,b+o.b};}
    cp operator -(const cp &o)const {return (cp){a-o.a,b-o.b};}
    cp operator *(const cp &o)const {return (cp){a*o.a-b*o.b,b*o.a+a*o.b};}
    cp operator *(const double &o)const {return (cp){a*o,b*o};}
    cp operator !() const{return (cp){a,-b};}
}w[maxfft];
int pos[maxfft];
void fft_init(int len)
{
    int j=0;
    while((1<<j)<len)j++;
    j--;
    for(int i=0;i<len;i++)
        pos[i]=pos[i>>1]>>1|((i&1)<<j);
}
void fft(cp *x,int len,int sta)
{
    for(int i=0;i<len;i++)
        if(i<pos[i])swap(x[i],x[pos[i]]);
    w[0]=(cp){1,0};
    for(unsigned i=2;i<=len;i<<=1)
    {
        cp g=(cp){cos(2*pi/i),sin(2*pi/i)*sta};
        for(int j=i>>1;j>=0;j-=2)w[j]=w[j>>1];
        for(int j=1;j<i>>1;j+=2)w[j]=w[j-1]*g;
        for(int j=0;j<len;j+=i)
        {
            cp *a=x+j,*b=a+(i>>1);
            for(int l=0;l<i>>1;l++)
            {
                cp o=b[l]*w[l];
                b[l]=a[l]-o;
                a[l]=a[l]+o;
            }
        }
    }
    if(sta==-1)for(int i=0;i<len;i++)x[i].a/=len,x[i].b/=len;
}
cp x[maxfft],y[maxfft],z[maxfft];
void FFT(int *a,int *b,int n,int m,int *c)
{
    int len=1;
    while(len<n+m)len<<=1;
    fft_init(len);
    for(int i=0;i<len;i++)
    {
        int aa=i<n?a[i]:0,bb=i<m?b[i]:0;
        x[i]=(cp){(aa>>15),(aa&32767)},y[i]=(cp){(bb>>15),(bb&32767)};
    }
    fft(x,len,1),fft(y,len,1);
    for(int i=0;i<len;i++)
    {
        int j=len-1&len-i;
        z[i]=((x[i]+!x[j])*(y[i]-!y[j])+(x[i]-!x[j])*(y[i]+!y[j]))*(cp){0,-0.25};
    }
    fft(z,len,-1);
    for(int i=0;i<n+m-1;i++)
    {
        ll ta=(ll)(z[i].a+0.5)%mod;
        ta=(ta<<15)%mod;
        c[i]=ta;
    }
    for(int i=0;i<len;i++)
    {
        int j=len-1&len-i;
        z[i]=(x[i]-!x[j])*(y[i]-!y[j])*(cp){-0.25,0}+(x[i]+!x[j])*(y[i]+!y[j])*(cp){0,0.25};
    }
    fft(z,len,-1);
    for(int i=0;i<n+m-1;i++)
    {
        ll ta=(ll)(z[i].a+0.5)%mod,tb=(ll)(z[i].b+0.5)%mod;
        ta=(ta+(tb<<30))%mod;
        c[i]=(c[i]+ta)%mod;
    }
}
int Pow(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=(ll)ans*a%mod;
        a=(ll)a*a%mod;
        b>>=1;
    }
    return ans;
}
#define maxn 30005
int fact[maxn],inv[maxn];
void init(int n=30000)
{
    fact[0]=1;
    for(int i=1;i<=n;i++)fact[i]=(ll)i*fact[i-1]%mod;
    inv[1]=1;
    for(int i=2;i<=n;i++)inv[i]=mod-(ll)(mod/i)*inv[mod%i]%mod;
    inv[0]=1;
    for(int i=1;i<=n;i++)inv[i]=(ll)inv[i]*inv[i-1]%mod;
}
int C(int n,int m)
{
    return (ll)fact[n]*inv[m]%mod*inv[n-m]%mod;
}
int f[maxfft],f1[maxfft];
void Solve(int n,int k)
{
    if(n==1)
    {
        f[0]=0;
        for(int i=1;i<=k;i++)f[i]=1;
        return ;
    }
    Solve(n/2,k);
    for(int i=1;i<=k;i++)
    {
        f[i]=(ll)f[i]*inv[i]%mod;
        f1[i]=(ll)f[i]*Pow(2,(ll)(n/2)*i%(mod-1))%mod;
    }
    FFT(f,f1,k+1,k+1,f);
    for(int i=1;i<=k;i++)f[i]=(ll)f[i]*fact[i]%mod;
    if(n&1)
    {
        for(int i=1;i<=k;i++)
        {
            f[i]=(ll)f[i]*inv[i]%mod*Pow(2,i)%mod;
            f1[i]=inv[i];
        }
        FFT(f,f1,k+1,k+1,f);
        for(int i=1;i<=k;i++)f[i]=(ll)f[i]*fact[i]%mod;
    }
}
int main()
{
    init();
    int n,k;
    while(~scanf("%d%d",&n,&k))
    {
        if(n>k)
        {
            printf("0\n");
            continue;
        }
        Solve(n,k);
        int ans=0;
        for(int i=1;i<=k;i++)ans=(ans+(ll)f[i]*C(k,i)%mod)%mod;
        printf("%d\n",ans);
    } 
    return 0;
}