qoj9479 And DNA

statement

很巧妙的 trick，称之为数值倍增法，在 CF1844G 中亦有提及。

下文用 \(\oplus\) 指代按位与，用 \(x_i\) 指代 \(a_i\) 的最低位。

这里，\(\oplus\) 操作是位独立的，因此我们不妨按位从低到高构造这个序列。

我们试着求出对于 \(M\) 的方案数 \(T(M)\)。

---

如果 \(M\) 是偶数。意味着对于每个 \(a_i\)，

• 如果 \(x_i=0\)，那么 \(x_{i-1}\oplus x_{i+1}=0\)，也即两边至少有一个 \(0\)。此时最低位贡献 \(0\)。
• 如果 \(x_i=1\)，那么 \(x_{i-1}\oplus x_{i+1}=1\)，也即两边都是 \(1\)。此时最低位贡献 \(2\)。

不难证明只有 \(x_i\) 全为 \(0\) 或全为 \(1\) 可以满足条件。

那么直接拆掉最后一位往下递归，有 \(T(2K)=T(K)+T(K-1)\)。

---

如果 \(M\) 是奇数。意味着对于每个 \(x_i\)，

• 如果 \(x_i=0\)，那么 \(x_{i-1}\oplus x_{i+1}=1\)，也即两边都是 \(1\)。
• 如果 \(x_i=1\)，那么 \(x_{i-1}\oplus x_{i+1}=0\)，也即两边至少有一个 \(0\)。

不管哪种情况最低位均贡献 \(1\)。

同样不难证明，条件满足当且仅当 \(x\) 能被拆成若干个 \(011\) 或者 \(01\) 拼成的串。

注意第一个段可以循环位移若干次。例如 \(01101\) 可以位移成 \(10101\)。

考虑其生成函数 \(\dfrac{2x^2+3x^3}{1-x^2-x^3}\)，矩阵快速幂求值即可求得方案数 \(P\)，这一项仅跟 \(n\) 有关。

那么有 \(T(2K+1)=P\times T(K)\)。

---

那么直接求值即可。

#include <iostream>
#include <algorithm>

const int N = 40, O = 998244353;
#define rep(i,a,b) for(int i(a);i<(b);++i)
auto incr = [](auto& x, auto&& y) { x += y, x -= (x >= O) * O; };
typedef long long i64;

int n, m;

static constexpr int M = 3;
struct matrix {
    int mat[M][M];
    inline int* operator[] (int k) { return mat[k]; }
    inline const int* operator[] (int k) const {return mat[k]; }
};

inline matrix operator * (const matrix &a, const matrix& b) {
    matrix c{0};
    rep(i, 0, M) rep(j, 0, M) rep(k, 0, M)
        incr(c[i][k], 1ll * a[i][j] * b[j][k] %O);
    return c;
}

matrix F {
    0, 1, 0, 
    1, 0, 1, 
    1, 0, 0,
};

matrix pF[N];
inline void init() {
    pF[0] = F;
    rep(i, 1, N) {
        pF[i] = pF[i-1] * pF[i-1];
    }
}

inline i64 calc(int n) {
    if(n <= 0 || n == 2 || n == 3) return std::max(n, 1);
    if(n == 1) return 0;
    matrix init {
        3, 2, 0,
        0, 0, 0, 
        0, 0, 0, 
    };
    i64 p = n - 3;
    for(int i = 0; i < N; ++i) {
        if(p >> i & 1) {
            init = init * pF[i];
        }
    }
    return init[0][0];
}

i64 P = 0;

i64 dfs(int n) {
    if(n == 0) return 1;
    if(n & 1) {
        return dfs(n / 2) * P %O;
    } else {
        return (dfs(n / 2) + dfs(n / 2 - 1)) %O;
    }
}

signed main() {
    std::ios::sync_with_stdio(0), std::cin.tie(0), std::cout.tie(0);
    std::cin >> n >> m;
    init(); P = calc(n); 
    std::cerr << P << "\n";
    std::cout << dfs(m) << "\n";
}