CodeForces - 623E（倍增+ntt）

吃花椒的妙酱

已于 2022-08-11 23:24:13 修改

阅读量475

点赞数

文章标签：算法图论

于 2022-08-11 23:17:14 首次发布

本文链接：https://blog.csdn.net/m0_53688600/article/details/126294722

版权

题目大意： 一个序列的变换定义为对于每个ai变成or_sum[i],or_sum[i]表示前i个数的或值和(ai的或和)，问有多少个长度为n的序列满足变换后是个严格递增序列(n<=1e18,k<=3e4)
思路：
设 $a i^{'}$ 表示变换后位置i上的数
想一下，不难发现变换后的ai’对应的二进制中1的数量比 $a_{i-1}$ 中的要多，
当n>k时显然ans=0。

将ai看成01串后，1的位置并不重要了，如果能求出f(i)，f(i)表示 $a_n'$ 中有i个1的合法序列数，那么 $ans=\sum_{i=n}^kf(i)*C(k,i)$
考虑dp,设dp[i][j]表示前i个数组成的合法序列数且 $a_i’$ 中有j个1，我们要求的就是f(i) = dp[n][i]。
考虑从i-1推到i，有转移方程 $dp[i][j]=\sum_{t=0}^{j-1}dp[i-1][t]*2^{t}*C(j,t)$
意义是 $a_{i-1}'$ 中钦定t个位置是1然后 $a_i'$ 中对应的这t个位置可选1可不选。
直接暴力转移是 $O(nk^2)$ 的。
展开一下组合数,移项得:
$\frac{dp[i][j]}{j!}=\sum_{t=0}^{j-1}\frac{dp[i-1][t]*2^{t}}{t!}*\frac{1}{(j-t)!}$
观察到后面是个卷积形式，用fft优化后是O(nklogk)，再加上n次fft的大肠数，差不多2e10以上的量级了，还是过不了。
上面的式子还可以继续拓展（自己想不到qaq，更多是从数学意义上推的）
$dp[n+m][j]=\sum_{t=0}^{j}dp[n][t]*dp[m][j-t]*2^{mt}*C(j,t)$
意义是长度为n+m的序列，前部分 $a_n$ 有t个1，后半部分有m个数，每个数这个t个1可选可不选，其他意义同上。
展开组合数，移项得：
$\frac{dp[n+m][j]}{j!}=\sum_{t=0}^{j}\frac{dp[n][t]*2^{mt}}{t!}*\frac{dp[m][j-t]}{(j-t)!}$
然后我们把m换成n得：
$\frac{dp[2n][j]}{j!}=\sum_{t=0}^{j}\frac{dp[n][t]*2^{nt}}{t!}*\frac{dp[n][j-t]}{(j-t)!}$
然后就可以愉快的倍增搞了，当然长度为奇数的情况下，需要当成m = n+1，单独卷一次。
这样只要卷logn次，复杂度O(klogklogn),由于模数是1e9+7，还要做任意模数。常数十分大，博主用的3模数ntt，常数直接起飞，6100ms过的233。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define IOS ios::sync_with_stdio(false),cin.tie(nullptr) 
#define _for(i,a,b) for(int i=(a) ;i<=(b) ;i++)
#define _rep(i,a,b) for(int i=(a) ;i>=(b) ;i--)
#define mst(v,s) memset(v,s,sizeof(v))
#define pii pair<int ,int >
#define pb(v) push_back(v)
#define all(v) v.begin(),v.end()
#define int long long
#define inf 0x3f3f3f3f
#define INF 0x7f7f7f7f7f7f7f7f
#define endl "\n"
#define fi first
#define se second
#define ls p<<1
#define rs p<<1|1
#define lson p<<1,l,mid
#define rson p<<1|1,mid+1,r
#define AC return 0
#define ldb long double
const int N=1e6+5;
const double eps=1e-8;
const int maxn = 6e6+10;
int Pow(int x,int d,int p){
    int tans = 1;
    if(d == 0)return 1%p;
    int a = Pow(x,d/2,p);
    tans = 1ll*a*a%p;
    if(d%2)tans = 1ll*tans*x%p;
    return tans%p;
}
typedef vector<int> Poly;//多项式定义
int F1[maxn],F2[maxn];
int rev[maxn];
void NTT(int * A,int lim,int opt,int p) {
    int i, j, k, m, gn, g, tmp;
    for(int i = 0; i < lim; ++i)rev[i] = (i & 1)*(lim >> 1) + (rev[i >> 1] >> 1);
    for(i = 0; i < lim; ++i)if (rev[i] < i) swap(A[i], A[rev[i]]);
    for(m = 2; m <= lim; m <<= 1) {
        k = m >> 1;
        gn = Pow(3,(p - 1) / m,p);
        for(i = 0; i < lim; i += m) {
            g = 1;
            for (j = 0; j < k; ++j, g = 1ll * g * gn % p) {
                tmp = 1ll * A[i + j + k] * g % p;
                A[i + j + k] = (A[i + j] - tmp + p) % p;
                A[i + j] = (A[i + j] + tmp) % p;
            }
        }
    }
    if(opt == -1){
        reverse(A+1,A+lim);
        int inv = Pow(lim,p-2,p);
        for(i = 0; i < lim; ++i) A[i] = 1ll * A[i] * inv % p;
    }
}
Poly mul(const Poly & A,const Poly & B,int p){
    int n = A.size(),m = B.size();
    int siz = n + m - 1;
    Poly C(siz);
    int fsiz = 1;
    while(fsiz <= siz)fsiz <<=1;
    for(int i = 0;i < fsiz;i++)F1[i] = F2[i] = 0;
    for(int i = 0;i < n;i++)F1[i] = A[i];
    for(int i = 0;i < m;i++)F2[i] = B[i];
    NTT(F1,fsiz,1,p);
    NTT(F2,fsiz,1,p);
    for(int i = 0;i < fsiz;i++)F1[i] = 1ll*F1[i]*F2[i]%p;
    NTT(F1,fsiz,-1,p);
    for(int i = 0;i < siz;i++){
        C[i] = F1[i];
    }
    return C;
}
ll ksc(ll x,ll y,ll p) {return (x*y-(ll)((long double)x/p*y+1e-8)*p+p)%p;}
int mm[3]={469762049,998244353,1004535809};
int n,m,k;
const int mod = 1e9+7;
Poly ans[3];
Poly cal(Poly &F,Poly &G,int mod){//任意模数多项式乘法
    int n = F.size(),m=G.size();
    Poly H(n+m-1,0);
    for(int i=0 ;i<3 ;i++) ans[i] = mul(F,G,mm[i]);
    int M = mm[0]*mm[1];
    int M0 = ksc(mm[1],Pow(mm[1],mm[0]-2,mm[0]),M);
    int M1 = ksc(mm[0],Pow(mm[0],mm[1]-2,mm[1]),M);
    for(int i =0 ;i<n+m-1 ;i++){
        ll A=(ksc(M0,ans[0][i],M)+ksc(M1,ans[1][i],M))%M;
        ll k=((ans[2][i]-A)%mm[2]+mm[2])%mm[2]*Pow(M%mm[2],mm[2]-2,mm[2])%mm[2];
        H[i] = ((M%mod)*(k%mod)%mod+A%mod)%mod;
    }
    return H;
}
int f[N],ni_f[N];
void ini(){
    f[0]=1;
    int mx = 3e4;
    _for(i,1,mx) f[i] =  (f[i-1] * i)%mod;
    ni_f[mx] = Pow(f[mx],mod-2,mod);
    _rep(i,mx-1,0) ni_f[i] = ni_f[i+1] * (i+1)%mod;
}
ll C(int n,int m){
    if( m==n || m==0 ) return 1;
    if( n < m ) return 0;
    return f[n] * ni_f[n-m]%mod * ni_f[m]%mod;
}
void fz(Poly &F){
    F.resize(k+1,0);
}
void solve(int len,Poly &F,Poly &G,int p){
    if( len==1 ){
        F[0]=1;
        _for(i,1,k) F[i] = 1;
        return;
    }
    fz(F),fz(G);
    solve(len>>1,F,G,p);
    F[0] = G[0] = 0;
    _for(i,1,k){
        F[i] = F[i] * ni_f[i]%p,
        G[i] = F[i] * Pow(2,(len>>1)*i%(p-1),p)%p;
    }
    fz(F),fz(G);
    F = cal(F,G,p);
    if( len&1 ){
        _for(i,1,k) {
            F[i] = F[i] * Pow(2,i,p)%p,
            G[i] = ni_f[i];
        }
        fz(F),fz(G);
        F = cal(F,G,p);
    }
    _for(i,1,k) {
        F[i] = F[i] * f[i]%p;
    }
}
signed main(){
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif  
    IOS;
    ini();
    cin>>n>>k;
    if( n > k ){
        cout<<0<<endl;
        return 0;
    }
    Poly F(k+1,0),G(k+1,0);
    solve(n,F,G,mod);
    int res = 0;
    for(int i=n ;i<=k ;i++){
        res = (res + C(k,i)*F[i]%mod)%mod;
    }
    cout<<res<<endl;
    AC; 
}