「CF623E」 Transforming Sequence

March 14, 2019

1323views

1737 words

没想到非 Chinese Round 也会出这么毒瘤的题目......

题目大意：对于正整数序列 $A$，定义序列 $B$： $B_1=A_1,B_i=B_{i-1}\ or\ A_i,i\in[2,n]B $ 其中 or 为位或运算。每一个序列合法，满足对于$ \forall i\in[1,n]$，有 $A_i\in[1,2^k])$ 而且对于 $\forall i\in[2,n]$，有$B_i>B_{i-1}$。给出 $n(n \leq 10^{18}),k(k \leq 3* 10^4)$ ，求合法序列数。

对 $10^9 + 7$ 取模。

首先我们观察到每次至少要多一个变成 $1$ 的二进制位否则不可能满足要求。

所以我们按照这个，自然地考虑动态规划：

设 $F_{i,j}$ 表示已经选择了 $i$ 个数，已经有 $j$ 位变成了 $1$ 的方案数。首先观察到如果这样设的话会超时，所以我们考虑如何把转移写成卷积的形式。我们考虑如何用已知的值更新 $f_{x+y,i}$。

我们设 $j \in [1,i]$。这个 $x$ 和 $y$ 的枚举实质就是枚举那部分钦定那部分不钦定，所以前 $x$ 个数中选择 $j$ 位的方案显然是 $(^i_j) \times f_{x,j}$（组合数是考虑在 $i$ 个二进制位中任选 $j$ 个）。后 $y$ 个数中需要选择 $i-j$ 位，方案数为 $f_{y,i-j}$，然后我们考虑这些前面 $x$ 个数钦定了 $j$ 个位置，那么剩下的 $y$ 个数上在这 $j$ 个位置上可以随便填（因为或运算保证已经有一个是 $1$ 了），所以需要再乘上一个 $(2^{j})^{y}$。

转移方程大概如下：

$$f_{x+y,i} = \sum_{j=1}^i f_{x,j} \times (^i_j) \times f_{y,i-j} \times (2^{j})^y$$

通过小学就学习过的基本变换，不难得到：

$$f_{x+y,i} = i! \times \sum_{j=1}^i \frac{f_{x,j} \times (2^{j})^y }{j!} \times \frac{f_{y,i-j}}{(i-j)!}$$

现在我们就可以用卷积优化转移了，不过还需要处理两个小问题。

$10^9 + 7$ 不是个 NTT 模数哎怎么办呢？

有个科技叫做任意模数 NTT（MTT），建议大家去点一下。

我们考虑两个多项式做乘法，不能用 FFT 的原因是数字过大会爆精度。不如考虑把每个多项式写成 $A_{1} \times pp + A_{2}$，取 $pp = \sqrt{p} $。这样两个多项式 $A,B$ 相乘就转化为了求 $(A_{1} \times pp + A_{2}) \times (B_{1} \times pp + B_{2})$。

展开可得：$A_{1} \times B_{1} \times pp^2 + (A_1 \times B_2 + A_2 \times B_1) \times pp + A_2 \times B_2$。多做几遍 DFT 最后在外面取模就好了。不过貌似有国家队神仙提出了四遍 DFT 做法但是我不会。

这个 $n$ 这么大怎么选择 $x,y$ 呢？

考虑快速幂一样的把答案滚出来。详情见我的代码。

#include <algorithm>
#include <iostream>
#include <cstdio>
#include <queue>
#include <cmath>
#include <cstring>

#define fi first
#define se second
#define LL long long
#define P std::pair<int,int>
#define CLR(a,b) memset(a,b,sizeof(a))
#define FOR(i,a,b) for(int i = a;i <= b;++i)
#define ROF(i,a,b) for(int i = a;i >= b;--i)
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl
#define int LL
const int MAXN = (1<<17)+2;
const double PI = acos(-1.0);
int ha,haha;

inline int qpow(int a,int n=ha-2){
    int res = 1;
    while(n){
        if(n & 1) res = 1ll*res*a%ha;
        a = 1ll*a*a%ha;
        n >>= 1;
    }
    return res;
}

struct cp{
    double x,y;//x+iy
    cp(){x=0,y=0;}
    cp(long double x,long double y):x(x),y(y){}
    cp operator + (const cp &t) const {
        return cp(x+t.x,y+t.y);
    }
    cp operator - (const cp &t) const {
        return cp(x-t.x,y-t.y);
    }
    cp operator * (const cp &t) const {
        return cp(x*t.x-y*t.y,x*t.y+y*t.x);
    }
    LL D(){
        LL t = (x+0.499);
        return t%ha;
    }
};
int N,n,k;
int r[MAXN];

cp w[1<<17][100];

void init(){
    for(int i=2,j=0;j<17;j++,i<<=1){
        FOR(k,0,i-1) w[k][j] = cp(cos(2*PI*k/i),sin(2*PI*k/i));
    }
}

inline void FFT(cp *A,int limit,int opt){
    FOR(i,0,limit) if(i < r[i]) std::swap(A[i],A[r[i]]);
    for(int mid = 1,i=0;mid < limit;mid<<=1,++i){
        cp W(cos(PI/mid),sin(PI/mid));
        for(int j = 0;j < limit;j += (mid<<1)){
            cp w(1,0);
            for(int k = 0;k < mid;++k,w = w*W){
                w = ::w[k][i]; if(opt==-1) w.y=-w.y;
                cp x = A[j+k],y = w*A[j+mid+k];
                A[j+k] = x+y;A[j+mid+k] = x-y;
            }
        }
    }
    if(opt == -1){
        FOR(i,0,limit-1) A[i].x = A[i].x/limit;
    }
}

struct Poly{
    LL a[MAXN];

    Poly(){CLR(a,0);}
    
    friend Poly operator * (const Poly &a,const Poly &b){
        Poly ans;
        static cp A[MAXN],B[MAXN],C[MAXN],D[MAXN],E[MAXN],F[MAXN],G[MAXN];
        CLR(A,0);CLR(B,0);CLR(C,0);CLR(D,0);
        FOR(i,0,N-1) A[i].x = a.a[i]%haha,B[i].x = a.a[i]/haha;
        FOR(i,0,N-1) C[i].x = b.a[i]%haha,D[i].x = b.a[i]/haha;
        FFT(A,N,1);FFT(B,N,1);FFT(C,N,1);FFT(D,N,1);
        FOR(i,0,N-1){
            E[i] = A[i]*C[i];
            F[i] = A[i]*D[i]+B[i]*C[i];
            G[i] = B[i]*D[i];
        }
        FFT(E,N,-1);FFT(F,N,-1);FFT(G,N,-1);
        FOR(i,0,N-1) ans.a[i] = (E[i].D()%ha+F[i].D()%ha*haha%ha+G[i].D()%ha*haha%ha*haha%ha)%ha;
        FOR(i,k+1,N-1) ans.a[i] = 0;
        return ans;
    }
};

int inv[MAXN],fac[MAXN];

inline void pre(){
    fac[0] = 1;
    FOR(i,1,k) fac[i] = 1ll*fac[i-1]*i%ha;
    inv[k] = qpow(fac[k]);
    ROF(i,k-1,0) inv[i] = 1ll*inv[i+1]*(i+1)%ha;
    init();
}

Poly ans,f,f1,f2;

inline int C(int up,int down){
    return fac[up]*inv[down]%ha*inv[up-down]%ha;
}

signed main(){
    ha = 1e9+7;haha = 32768;
    scanf("%lld%lld%lld",&n,&k);--n;
    int len = 0;
    N = 1;pre();
    for(N=1;N<=(k*2);N<<=1,len++);
    FOR(i,0,N) r[i] = (r[i>>1]>>1)|((i&1)<<(len-1));
    FOR(i,1,k) f.a[i] = 1;
    ans = f;
    LL now = 1;
    while(n){
        if(n & 1){
            f1 = ans;f2 = f;
            FOR(i,1,k) f1.a[i] = f1.a[i]*inv[i]%ha*qpow(qpow(2,i),now)%ha;
            FOR(i,1,k) f2.a[i] = f2.a[i]*inv[i]%ha;
            ans = f1*f2;
           // FOR(i,1,k) if(ans.a[i]<0 && n == 29999) DEBUG(ans.a[i]);
            FOR(i,1,k) ans.a[i] = ans.a[i]*fac[i]%ha;
        }
        f1 = f2 = f;
        FOR(i,1,k) f1.a[i] = f1.a[i]*inv[i]%ha*qpow(qpow(2,i),now)%ha;
        FOR(i,1,k) f2.a[i] = f2.a[i]*inv[i]%ha;
        f = f1*f2;
        FOR(i,1,k) f.a[i] = f.a[i]*fac[i]%ha;
        n >>= 1;now <<= 1;
    }
    LL sum = 0;
    FOR(i,1,k) sum = (sum + ans.a[i]%ha*C(k,i)%ha)%ha;
    printf("%lld\n",sum);
    return 0;
}

「CF623E」 Transforming Sequence

RainAir • 2019 年 03 月 14 日

没想到非 Chinese Round 也会出这么毒瘤的题目......

对 $10^9 + 7$ 取模。

首先我们观察到每次至少要多一个变成 $1$ 的二进制位否则不可能满足要求。

所以我们按照这个，自然地考虑动态规划：

转移方程大概如下：

$$f_{x+y,i} = \sum_{j=1}^i f_{x,j} \times (^i_j) \times f_{y,i-j} \times (2^{j})^y$$

通过小学就学习过的基本变换，不难得到：

$$f_{x+y,i} = i! \times \sum_{j=1}^i \frac{f_{x,j} \times (2^{j})^y }{j!} \times \frac{f_{y,i-j}}{(i-j)!}$$

现在我们就可以用卷积优化转移了，不过还需要处理两个小问题。

$10^9 + 7$ 不是个 NTT 模数哎怎么办呢？

有个科技叫做任意模数 NTT（MTT），建议大家去点一下。

这个 $n$ 这么大怎么选择 $x,y$ 呢？

考虑快速幂一样的把答案滚出来。详情见我的代码。

#include <algorithm>
#include <iostream>
#include <cstdio>
#include <queue>
#include <cmath>
#include <cstring>

#define fi first
#define se second
#define LL long long
#define P std::pair<int,int>
#define CLR(a,b) memset(a,b,sizeof(a))
#define FOR(i,a,b) for(int i = a;i <= b;++i)
#define ROF(i,a,b) for(int i = a;i >= b;--i)
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl
#define int LL
const int MAXN = (1<<17)+2;
const double PI = acos(-1.0);
int ha,haha;

inline int qpow(int a,int n=ha-2){
    int res = 1;
    while(n){
        if(n & 1) res = 1ll*res*a%ha;
        a = 1ll*a*a%ha;
        n >>= 1;
    }
    return res;
}

struct cp{
    double x,y;//x+iy
    cp(){x=0,y=0;}
    cp(long double x,long double y):x(x),y(y){}
    cp operator + (const cp &t) const {
        return cp(x+t.x,y+t.y);
    }
    cp operator - (const cp &t) const {
        return cp(x-t.x,y-t.y);
    }
    cp operator * (const cp &t) const {
        return cp(x*t.x-y*t.y,x*t.y+y*t.x);
    }
    LL D(){
        LL t = (x+0.499);
        return t%ha;
    }
};
int N,n,k;
int r[MAXN];

cp w[1<<17][100];

void init(){
    for(int i=2,j=0;j<17;j++,i<<=1){
        FOR(k,0,i-1) w[k][j] = cp(cos(2*PI*k/i),sin(2*PI*k/i));
    }
}

inline void FFT(cp *A,int limit,int opt){
    FOR(i,0,limit) if(i < r[i]) std::swap(A[i],A[r[i]]);
    for(int mid = 1,i=0;mid < limit;mid<<=1,++i){
        cp W(cos(PI/mid),sin(PI/mid));
        for(int j = 0;j < limit;j += (mid<<1)){
            cp w(1,0);
            for(int k = 0;k < mid;++k,w = w*W){
                w = ::w[k][i]; if(opt==-1) w.y=-w.y;
                cp x = A[j+k],y = w*A[j+mid+k];
                A[j+k] = x+y;A[j+mid+k] = x-y;
            }
        }
    }
    if(opt == -1){
        FOR(i,0,limit-1) A[i].x = A[i].x/limit;
    }
}

struct Poly{
    LL a[MAXN];

    Poly(){CLR(a,0);}
    
    friend Poly operator * (const Poly &a,const Poly &b){
        Poly ans;
        static cp A[MAXN],B[MAXN],C[MAXN],D[MAXN],E[MAXN],F[MAXN],G[MAXN];
        CLR(A,0);CLR(B,0);CLR(C,0);CLR(D,0);
        FOR(i,0,N-1) A[i].x = a.a[i]%haha,B[i].x = a.a[i]/haha;
        FOR(i,0,N-1) C[i].x = b.a[i]%haha,D[i].x = b.a[i]/haha;
        FFT(A,N,1);FFT(B,N,1);FFT(C,N,1);FFT(D,N,1);
        FOR(i,0,N-1){
            E[i] = A[i]*C[i];
            F[i] = A[i]*D[i]+B[i]*C[i];
            G[i] = B[i]*D[i];
        }
        FFT(E,N,-1);FFT(F,N,-1);FFT(G,N,-1);
        FOR(i,0,N-1) ans.a[i] = (E[i].D()%ha+F[i].D()%ha*haha%ha+G[i].D()%ha*haha%ha*haha%ha)%ha;
        FOR(i,k+1,N-1) ans.a[i] = 0;
        return ans;
    }
};

int inv[MAXN],fac[MAXN];

inline void pre(){
    fac[0] = 1;
    FOR(i,1,k) fac[i] = 1ll*fac[i-1]*i%ha;
    inv[k] = qpow(fac[k]);
    ROF(i,k-1,0) inv[i] = 1ll*inv[i+1]*(i+1)%ha;
    init();
}

Poly ans,f,f1,f2;

inline int C(int up,int down){
    return fac[up]*inv[down]%ha*inv[up-down]%ha;
}

signed main(){
    ha = 1e9+7;haha = 32768;
    scanf("%lld%lld%lld",&n,&k);--n;
    int len = 0;
    N = 1;pre();
    for(N=1;N<=(k*2);N<<=1,len++);
    FOR(i,0,N) r[i] = (r[i>>1]>>1)|((i&1)<<(len-1));
    FOR(i,1,k) f.a[i] = 1;
    ans = f;
    LL now = 1;
    while(n){
        if(n & 1){
            f1 = ans;f2 = f;
            FOR(i,1,k) f1.a[i] = f1.a[i]*inv[i]%ha*qpow(qpow(2,i),now)%ha;
            FOR(i,1,k) f2.a[i] = f2.a[i]*inv[i]%ha;
            ans = f1*f2;
           // FOR(i,1,k) if(ans.a[i]<0 && n == 29999) DEBUG(ans.a[i]);
            FOR(i,1,k) ans.a[i] = ans.a[i]*fac[i]%ha;
        }
        f1 = f2 = f;
        FOR(i,1,k) f1.a[i] = f1.a[i]*inv[i]%ha*qpow(qpow(2,i),now)%ha;
        FOR(i,1,k) f2.a[i] = f2.a[i]*inv[i]%ha;
        f = f1*f2;
        FOR(i,1,k) f.a[i] = f.a[i]*fac[i]%ha;
        n >>= 1;now <<= 1;
    }
    LL sum = 0;
    FOR(i,1,k) sum = (sum + ans.a[i]%ha*C(k,i)%ha)%ha;
    printf("%lld\n",sum);
    return 0;
}

「CF623E」 Transforming Sequence

Leave a Comment Cancel reply

NOI2021 退役记

I'm back

SDOI2020 退役记

点双连通分量 & 圆方树学习笔记

CSP2019 游记

集合幂级数学习笔记

20联赛集训day2 题解

20联赛集训day5 题解

20联赛集训day3 题解

线性基学习笔记

「CF623E」 Transforming Sequence