【LOJ2325】【清华集训 2017】小 Y 和恐怖的奴隶主（期望DP，矩阵快速幂）

Description

https://loj.ac/problem/2325

---

Solution

首先了解一下这题的弱化版：抵制克苏恩
上面的博客中计算了一个期望$f$和一个概率$p$，其实只要计算$p$，答案为$\sum \frac{p_{t,i,j,k}}{i+j+k+1}$。
对于$p$，每次的转移都是一样的，可以用矩乘快速幂优化。
至于答案，我们给矩阵加上一行一列作为计数器即可。

---

Code

/************************************************
 * Au: Hany01
 * Date: Aug 20th, 2018
 * Prob: [LOJ2325][清华集训2017]小Y和恐怖的奴隶主
 * Email: hany01dxx@gmail.com & hany01@foxmail.com
 * Inst: Yali High School
************************************************/

#include<bits/stdc++.h>

using namespace std;

typedef unsigned long long LL;
typedef long double LD;
typedef pair<uint, uint> PII;
typedef unsigned int uint;
#define rep(i, j) for (register uint i = 0, i##_end_ = (j); i <= i##_end_; ++ i)
#define For(i, j, k) for (register uint i = (j), i##_end_ = (k); i <= i##_end_; ++ i)
#define Fordown(i, j, k) for (register uint i = (j), i##_end_ = (k); i >= i##_end_; -- i)
#define Set(a, b) memset(a, b, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define x first
#define y second
#define pb(a) push_back(a)
#define mp(a, b) make_pair(a, b)
#define SZ(a) ((uint)(a).size())
#define INF (0x3f3f3f3f)
#define INF1 (2139062143)
#define debug(...) fprintf(stderr, __VA_ARGS__)
#define y1 wozenmezhemecaia
#define MOD (998244353)
#define MOD1 (9223372036388749853ll)

template <typename T> inline bool chkmax(T &a, T b) { return a < b ? a = b, 1 : 0; }
template <typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }

namespace pb_ds {

    namespace io {
        const int MaxBuff = 1 << 16;
        const int Output = 1 << 22;
        char B[MaxBuff], *S = B, *T = B;
        #define getc() ((S == T) && (T = (S = B) + fread(B, 1, MaxBuff, stdin), S == T) ? 0 : *S++)
        char Out[Output], *iter = Out;
        inline void flush() {
            fwrite(Out, 1, iter - Out, stdout);
            iter = Out;
        }
    }

    template<class Type> inline Type read() {
        using namespace io;
        register char ch; register Type ans = 0; register bool neg = 0;
        while(ch = getc(), (ch < '0' || ch > '9') && ch != '-')     ;
        ch == '-' ? neg = 1 : ans = ch - '0';
        while(ch = getc(), '0' <= ch && ch <= '9') ans = ans * 10 + ch - '0';
        return neg ? -ans : ans;
    }

    template<class Type> inline void write(register Type x, register char ch = '\n') {
        using namespace io;
        if(!x) *iter++ = '0';
        else {
            if(x < 0) *iter++ = '-', x = -x;
            static int s[100]; register int t = 0;
            while(x) s[++t] = x % 10, x /= 10;
            while(t) *iter++ = '0' + s[t--];
        }
        *iter++ = ch;
    }
}
using namespace pb_ds;

uint m, K, id1[9], id2[9][9], id3[9][9][9], N, rec[9];
LL  n;

inline uint Pow(uint a, uint b) {
    static uint Ans;
    for (Ans = 1, a %= MOD; b; b >>= 1, a = (LL)a * a % MOD) if (b & 1) Ans = (LL)Ans * a % MOD;
    return Ans;
}
inline uint ad(uint x, uint y) { if ((x += y) >= MOD) return x - MOD; return x; }

inline void DFS(uint cur, uint res) {
	if (cur == m + 1) {
		if (m == 1) id1[rec[1]] = N ++;
		else if (m == 2) id2[rec[1]][rec[2]] = N ++;
		else id3[rec[1]][rec[2]][rec[3]] = N ++;
		return;
	}
	For(i, 0, res) rec[cur] = i, DFS(cur + 1, res - i);
}

struct Matrix {
    uint a[167][167], x, y;
    inline uint* operator [] (uint x) { return a[x]; }
}Ans, M[77];
inline Matrix operator * (Matrix A, Matrix B) {
    static Matrix C; C.x = A.x, C.y = B.y;
	static LL fuck;
    rep(i, C.x) rep(j, C.y) {
		fuck = 0;
        rep(k, A.y) {
			fuck += (LL)A[i][k] * B[k][j];
			if (fuck >= MOD1) fuck %= MOD;
		}
		C[i][j] = fuck % MOD;
    }
    return C;
}

int main()
{
#ifdef hany01
	freopen("loj2325.in", "r", stdin);
	freopen("loj2325.out", "w", stdout);
#endif

	static uint T = read<int>();
	m = read<int>(), K = read<int>(), DFS(1, K);
	if (m == 1) {
		For(i, 0, K) {
			register uint t = Pow(i + 1, MOD - 2);
			M[0][id1[i]][id1[i]] = t;
			if (i) M[0][id1[i - 1]][id1[i]] = (LL)t * i % MOD;
			M[0][N][id1[i]] = t;
		}
		M[0][N][N] = 1, M[0].x = M[0].y = N;
	} else if (m == 2) {
		For(i, 0, K) For(j, 0, K - i) {
			register uint t = Pow(i + j + 1, MOD - 2);
			M[0][id2[i][j]][id2[i][j]] = t;
			if (i) M[0][id2[i - 1][j]][id2[i][j]] = (LL)t * i % MOD;
			if (j) M[0][id2[i + 1][j - 1 + (i + j < K)]][id2[i][j]] = (LL)t * j % MOD;
			M[0][N][id2[i][j]] = t;
		}
		M[0][N][N] = 1, M[0].x = M[0].y = N;
	} else if (m == 3) {
		For(i, 0, K) For(j, 0, K - i) For(k, 0, K - i - j) {
			register uint t = Pow(i + j + k + 1, MOD - 2);
			M[0][id3[i][j][k]][id3[i][j][k]] = t;
			if (i) M[0][id3[i - 1][j][k]][id3[i][j][k]] = (LL)t * i % MOD;
			if (j) M[0][id3[i + 1][j - 1][k + (i + j + k < K)]][id3[i][j][k]] = (LL)t * j % MOD;
			if (k) M[0][id3[i][j + 1][k - 1 + (i + j + k < K)]][id3[i][j][k]] = (LL)t * k % MOD;
			M[0][N][id3[i][j][k]] = t;
		}
		M[0][N][N] = 1, M[0].x = M[0].y = N;
	}
	for (register uint t = 1; (1ll << t) <= 1e18; ++ t) M[t] = M[t - 1] * M[t - 1];

#ifdef hany01
	cerr << clock() * 1. / CLOCKS_PER_SEC << endl;
#endif

	while (T --) {
		n = read<LL>(), Ans.x = N, Ans.y = 0;
		rep(i, N) Ans[i][0] = 0;
		if (m == 1) Ans[id1[1]][0] = 1;
		else if (m == 2) Ans[id2[0][1]][0] = 1;
		else if (m == 3) Ans[id3[0][0][1]][0] = 1;
		for (register uint i = 0; i < 60; ++ i) if (n >> i & 1) Ans = M[i] * Ans;
		printf("%d\n", Ans[N][0]);
	}

#ifdef hany01
	cerr << clock() * 1. / CLOCKS_PER_SEC << endl;
#endif

	return 0;
}