CodeForces 1152F2 Neko Rules the Catniverse (Large Version)

题目链接：http://codeforces.com/problemset/problem/1152/F2

题目大意

　　见http://codeforces.com/problemset/problem/1152/F1，此题 n 最大能到 10⁹。

分析

　　在 F1 的基础上，我们发现 dp[i + 1] 数组的每个值均可以通过 dp[i] 数组的有限几个数求得，而 dp[i + 2] 数组的每个值也均可以通过 dp[i + 1] 数组相同位置的有限几个数求得，于是我们可以考虑用矩阵快速幂来求 dp[n]。

　　由于 dp 数组的第二维和第三维数值并不是特别大，我们可以把 dp 数组的后两个维度压制成一个维度，于是 dp[i][j][sta] 就变成 dp[i][cur]，其中 cur = j * (1 << m) + sta。

　　假设 dp[i][cur] 对某一个 dp[i + 1][newcur] 有贡献，根据上一题的帖子，有两种可能：

1. 不访问：dp[i + 1][newcur] += dp[i][cur]。
2. 访问：dp[i + 1][newcur] += dp[i][cur] * (1 + 后m个星球被访问过的星球个数)。

　　定义 base 矩阵是我们所要构造的矩阵，设 base[cur][newcur] 代表某一个状态 cur 对一个新状态 newcur 所做贡献的系数，于是有：

1. 不访问：base[cur][newcur] = 1。
2. 访问：base[cur][newcur] =  (1 + 后m个星球被访问过的星球个数)。
3. 其他(cur 压根不可能变成 newcur)：base[cur][newcur] = 0。

　　于是 $dp[i + 1][newcur] = sum_{cur = 0}^{maxCur} dp[i][cur] * base[cur][newcur]$。

　　根据这个式子可以构造 dp 矩阵如下（以 k = 1， m = 1 为例）：

$$
dp(i) = egin{bmatrix}
dp[i][0] & dp[i][1] & dp[i][2] & dp[i][3]
end{bmatrix} 
$$

　　相应 base 系数矩阵如下：

$$
base = egin{bmatrix}
base[0][0] & base[0][1] & base[0][2] & base[0][3] \
base[1][0] & base[1][1] & base[1][2] & base[1][3] \
base[2][0] & base[2][1] & base[2][2] & base[2][3] \
base[3][0] & base[3][1] & base[3][2] & base[3][3]
end{bmatrix}
$$

　　base 矩阵填上数后为：

$$
base = egin{bmatrix}
1 & 1 & 0 & 0 \
0 & 0 & 0 & 0 \
0 & 0 & 1 & 1 \
1 & 2 & 0 & 0
end{bmatrix}
$$

　　于是有：

$$
egin{align*}
dp(i) &= dp(i - 1) * base \
dp(n) &= dp(0) * base^{n}
end{align*}
$$

　　不知道关于矩阵怎么构造的说的明不明白，有不足的地方还请指出。

　　PS：如果不压制维度，那么 base 矩阵就要弄四维的了，四维矩阵乘法反正我是不会。

代码如下

　　时间复杂度：$O(log n*(k*2^m)^3)$

#include <bits/stdc++.h>
using namespace std;
 
#define INIT() ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
#define Rep(i,n) for (int i = 0; i < (n); ++i)
#define For(i,s,t) for (int i = (s); i <= (t); ++i)
#define rFor(i,t,s) for (int i = (t); i >= (s); --i)
#define ForLL(i, s, t) for (LL i = LL(s); i <= LL(t); ++i)
#define rForLL(i, t, s) for (LL i = LL(t); i >= LL(s); --i)
#define foreach(i,c) for (__typeof(c.begin()) i = c.begin(); i != c.end(); ++i)
#define rforeach(i,c) for (__typeof(c.rbegin()) i = c.rbegin(); i != c.rend(); ++i)
 
#define pr(x) cout << #x << " = " << x << "  "
#define prln(x) cout << #x << " = " << x << endl
 
#define LOWBIT(x) ((x)&(-x))
 
#define ALL(x) x.begin(),x.end()
#define INS(x) inserter(x,x.begin())
 
#define ms0(a) memset(a,0,sizeof(a))
#define msI(a) memset(a,inf,sizeof(a))
#define msM(a) memset(a,-1,sizeof(a))

#define MP make_pair
#define PB push_back
#define ft first
#define sd second
 
template<typename T1, typename T2>
istream &operator>>(istream &in, pair<T1, T2> &p) {
    in >> p.first >> p.second;
    return in;
}
 
template<typename T>
istream &operator>>(istream &in, vector<T> &v) {
    for (auto &x: v)
        in >> x;
    return in;
}
 
template<typename T1, typename T2>
ostream &operator<<(ostream &out, const std::pair<T1, T2> &p) {
    out << "[" << p.first << ", " << p.second << "]" << "
";
    return out;
}

inline int gc(){
    static const int BUF = 1e7;
    static char buf[BUF], *bg = buf + BUF, *ed = bg;
    
    if(bg == ed) fread(bg = buf, 1, BUF, stdin);
    return *bg++;
} 

inline int ri(){
    int x = 0, f = 1, c = gc();
    for(; c<48||c>57; f = c=='-'?-1:f, c=gc());
    for(; c>47&&c<58; x = x*10 + c - 48, c=gc());
    return x*f;
}
 
typedef long long LL;
typedef unsigned long long uLL;
typedef pair< double, double > PDD;
typedef pair< int, int > PII;
typedef pair< string, int > PSI;
typedef set< int > SI;
typedef vector< int > VI;
typedef map< int, int > MII;
typedef pair< LL, LL > PLL;
typedef vector< LL > VL;
typedef vector< VL > VVL;
const double EPS = 1e-10;
const LL inf = 0x7fffffff;
const LL infLL = 0x7fffffffffffffffLL;
const LL mod = 1e9 + 7;
const int maxN = 1e5 + 7;
const LL ONE = 1;
const LL evenBits = 0xaaaaaaaaaaaaaaaa;
const LL oddBits = 0x5555555555555555;

void add_mod(LL &a, LL b) {
    a = (a + b) % mod;
}

struct Matrix{
    int row, col;
    LL MOD;
    VVL mat;
    
    Matrix(int r, int c, LL p = mod) : row(r), col(c), MOD(p) { 
        mat.assign(r, VL(c, 0));
    }
    Matrix(const Matrix &x, LL p = mod) : MOD(p){
        mat = x.mat;
        row = x.row;
        col = x.col;
    }
    Matrix(const VVL &A, LL p = mod) : MOD(p){
        mat = A;
        row = A.size();
        col = A[0].size();
    }
    
    // x * 单位阵 
    inline void E(int x = 1) {
        assert(row == col);
        Rep(i, row) mat[i][i] = x;
    }
    
    inline VL& operator[] (int x) {
        assert(x >= 0 && x < row);
        return mat[x];
    }
    
    inline Matrix operator= (const VVL &x) {
        row = x.size();
        col = x[0].size();
        mat = x;
        return *this;
    }
    
    inline Matrix operator+ (const Matrix &x) {
        assert(row == x.row && col == x.col);
        Matrix ret(row, col);
        Rep(i, row) {
            Rep(j, col) {
                ret.mat[i][j] = mat[i][j] + x.mat[i][j];
                ret.mat[i][j] %= MOD;
            }
        }
        return ret;
    }
    
    inline Matrix operator* (const Matrix &x) {
        assert(col == x.row);
        Matrix ret(row, x.col);
        Rep(k, x.col) {
            Rep(i, row) {
                if(mat[i][k] == 0) continue;
                Rep(j, x.col) {
                    ret.mat[i][j] += mat[i][k] * x.mat[k][j];
                    ret.mat[i][j] %= MOD;
                }
            }
        }
        return ret;
    }
    
    inline Matrix operator*= (const Matrix &x) { return *this = *this * x; }
    inline Matrix operator+= (const Matrix &x) { return *this = *this + x; }
    
    inline void print() {
        Rep(i, row) {
            Rep(j, col) {
                cout << mat[i][j] << " ";
            }
            cout << endl;
        }
    }
};

// 矩阵快速幂，计算x^y 
inline Matrix mat_pow_mod(Matrix x, LL y) {
    Matrix ret(x.row, x.col);
    ret.E();
    while(y){
        if(y & 1) ret *= x;
        x *= x;
        y >>= 1;
    }
    return ret;
}

int n, k, m, maxSta; 
LL ans;

int toId(int j, int sta) {
    return j * maxSta + sta;
}

int main(){
    INIT(); 
    cin >> n >> k >> m;
    maxSta = 1 << m; 
    Matrix dp(1, (k + 1) * maxSta);
    dp.mat[0][0] = 1; 
    Matrix base((k + 1) * maxSta, (k + 1) * maxSta);
    
    // 构造base矩阵 
    Rep(j, k + 1) {
        Rep(sta, maxSta) {
            int newsta = (sta << 1) % maxSta;
            int cur = toId(j, sta), newcur = toId(j, newsta);
            
            // 不选 i + 1 个星球从cur->newcur的状态变化 
            base.mat[cur][newcur] = 1;
            // 选 i + 1 个星球从cur->newcur的状态变化 
            if (j < k) {
                newcur = toId(j + 1, newsta | 1);
                base.mat[cur][newcur] = __builtin_popcount(sta) + 1;
            }
        }
    }
    
    base = mat_pow_mod(base, n);
    dp *= base;
    
    Rep(sta, maxSta) add_mod(ans, dp.mat[0][toId(k, sta)]);
    cout << ans << endl;
    return 0;
}