Codeforces Round #536 (Div. 2) F. Lunar New Year and a Recursive Sequence

\[f_i=\left(\prod_{j=1}^{k}f_{i-j}^{b_j}\right) \bmod p \]

因为这个质数的原根是 \(3\)，\(f_i\) 就能被 \(3^{a_i}\) 表示
所以就是 $$\begin{aligned} 3^{{a_i}=&\left(\prod_{j=1}}3^{a_{i-j}b_j}\right) \bmod p\=&3^{\sum_{j=1}ka_{i-j}b_j} \bmod p \end{aligned}$$
那么就可以用矩阵预处理出递推矩阵，注意相乘过程中模数是 \(p-1\)
然后设 \(f_k=3^x\)，\(m=3^b\)，\(b\) 可以用 BSGS 求出，转移矩阵的那个元素为 \(a\)
即求 \(ax\equiv b \pmod {p-1}\) 的解，最后答案就是 \(3^{x}\)

#include <bits/stdc++.h>
#define ll long long

const int MOD = 998244353, INF = 0x3f3f3f3f;
const int N = 110;
int len, n, m;
int b[N];

struct Mat {
	int mat[N][N];
	Mat(int x = 0) {
		memset(mat, 0, sizeof(mat));
		for (int i = 1; i <= n; i++)
			mat[i][i] = x;
	}
	Mat operator * (const Mat &p) const {
		Mat c;
		for (int i = 1; i <= n; i++) 
			for (int j = 1; j <= n; j++) {
				ll temp = 0;
				for (int k = 1; k <= n; k++) {
					temp += 1LL * mat[i][k] * p.mat[k][j];
					if (temp >= (1LL << 59))
						temp %= (MOD - 1);
				}
				c.mat[i][j] = temp % (MOD - 1);
			}
		return c;
	}
};

Mat qp(Mat a, int b) {
	Mat c(1);
	for (; b; b >>= 1, a = a * a)
		if (b & 1)
			c = c * a;
	return c;
}

int qp(int a, int b, int p = MOD) {
	a %= p;
	if (!a) return 0;
	int res = 1;
	while (b > 0) {
		if (b & 1) res = 1LL * res * a % p;
		a = 1LL * a * a % p;
		b >>= 1;
	}
	return res % p;
}

int gcd(int a, int b) {
	while (b) {
		a %= b;
		std::swap(a, b);
	}
	return a;
}

struct Ha {
	static const int mod = 1e6 + 7;
	int head[mod + 5], cnt;
	struct E {
		int x, ans, ne;
	} e[mod << 1];
	void clear() {
		cnt = 1;
		memset(head, 0, sizeof(head));
	}
	Ha() { clear(); }
	void insert(int x, int v) {
		int u = x % mod;
		e[++cnt].x = x; e[cnt].ans = v; e[cnt].ne = head[u]; head[u] = cnt;
	}
	int operator[](const int &k) const {
		int u = k % mod;
		for (int i = head[u]; i; i = e[i].ne)
			if (e[i].x == k) return e[i].ans;
		return -1;
	}
} H;

int BSGS(int a, int b, int mod = MOD) {
	if (mod == 1) return 0;
	if (a >= mod) a %= mod;
	if (b >= mod) b %= mod;
	if (a == 0) return b == 0 ? 1 : -1;
	if (b == 1) return 0;
	int m = ceil(sqrt(mod + 0.5));
	H.clear();
	for (int B = 0; B < m; B++) {
		H.insert(b, B);
		b = 1LL * b * a % mod;
	}
	int base = qp(a, m, mod), cur = base;
	for (int A = 1; A <= m + 1; A++) {
		int B = H[cur];
		if (~B) return A * m - B;
		cur = 1LL * cur * base % mod;
	}
	return -1;
}

int exgcd(int a, int b, int &x, int &y) {
	if (!b) {
		x = 1, y = 0;
		return a;
	}
	int g = exgcd(b, a % b, x, y);
	int temp = x;
	x = y;
	y = temp - a / b * y;
	return g;
}

void equ(int a, int b, int c, int &x, int &y) {
	int g = exgcd(a, b, x, y);
	if (c % g) {
		x = y = -INF;
		return;
	}
	a /= g, b /= g, c /= g;
	x = (x % b + b) % b;
	x = 1LL * x * c % b;
	y = b;
}

int main() {
	scanf("%d", &n);
	for (int i = 1; i <= n; i++) 
		scanf("%d", b + i);
	Mat base;
	for (int i = 1; i <= n; i++)
		base.mat[1][i] = b[i];
	for (int i = 2; i <= n; i++)
		base.mat[i][i - 1] = 1;
	scanf("%d%d", &len, &m);
	base = qp(base, len - n);
	int t = base.mat[1][1], b = BSGS(3, m);
	if (b == -1) {
		puts("-1");
		return 0;
	}
	int x, y;
	equ(t, MOD - 1, b, x, y);
	printf("%d\n", x == -INF ? -1 : qp(3, x));
	return 0;
}