没空写了，回头补下。
先放个板子。
struct Node {
Node operator*(Node b) {
// ...
}
};
Node pow(Node a, long long b) {
Node ans;
while (b) {
if (b & 1) ans = ans * a;
a = a * a;
b >>= 1;
}
return ans;
}
Node euclid(long long p, long long q, long long r, long long n, Node U, Node R) {
if (!n) return Node();
if (r >= q) return pow(U, r / q) * euclid(p, q, r % q, n, U, R);
if (p >= q) return euclid(p % q, q, r, n, U, pow(U, p / q) * R);
long long m = ((__int128_t) p * n + r) / q;
if (!m) return pow(R, n);
return pow(R, (q - r - 1) / p) * U * euclid(q, p, (q - r - 1) % p, m - 1, R, U)
* pow(R, n - ((__int128_t) q * m - r - 1) / p);
}
#include <bits/stdc++.h>
using namespace std;
const int P = 998244353;
struct Node {
long long x, y, sumx, sumy, sumy2, sumxy;
Node() : x(0), y(0), sumx(0), sumy(0), sumy2(0), sumxy(0) {}
Node operator*(Node b) {
Node a = *this, c;
c.x = (a.x + b.x) % P;
c.y = (a.y + b.y) % P;
c.sumx = (a.sumx + b.sumx + a.x * b.x) % P;
c.sumy = (a.sumy + b.sumy + a.y * b.x) % P;
c.sumy2 = (a.sumy2 + b.sumy2 + 2 * a.y * b.sumy + b.x * a.y % P * a.y) % P;
c.sumxy = (a.sumxy + b.sumxy + a.y * b.sumx + a.x * b.sumy + b.x * a.x % P * a.y) % P;
return c;
}
};
Node pow(Node a, int b) {
Node ans;
while (b) {
if (b & 1) ans = ans * a;
a = a * a;
b >>= 1;
}
return ans;
}
Node euclid(long long p, long long q, long long r, long long n, Node U, Node R) {
if (!n) return Node();
if (r >= q) return pow(U, r / q) * euclid(p, q, r % q, n, U, R);
if (p >= q) return euclid(p % q, q, r, n, U, pow(U, p / q) * R);
long long m = (1ll * p * n + r) / q;
if (!m) return pow(R, n);
return pow(R, (q - r - 1) / p) * U * euclid(q, p, (q - r - 1) % p, m - 1, R, U) * pow(R, n - (q * m - r - 1) / p);
}
int T;
int n, a, b, c;
int main() {
scanf("%d", &T);
Node U, R;
U.y = 1, R.x = 1, R.sumx = 1;
while (T--) {
scanf("%d%d%d%d", &n, &a, &b, &c);
auto ans = euclid(a, c, b, n, U, R);
ans.sumy = (ans.sumy + (b / c)) % P;
ans.sumy2 = (ans.sumy2 + 1ll * (b / c) * (b / c)) % P;
printf("%lld %lld %lld\n", ans.sumy, ans.sumy2, ans.sumxy);
}
return 0;
}
#include <bits/stdc++.h>
using namespace std;
const int P = 998244353;
const int MAXN = 22;
int n;
struct Matrix {
int a[MAXN][MAXN];
int* operator[](int b) { return a[b]; }
const int* operator[](const int b) const { return a[b]; }
Matrix() { memset(a, 0, sizeof a); }
Matrix(int n) { memset(a, 0, sizeof a); for (int i = 1; i <= n; i++) a[i][i] = 1; }
Matrix operator*(const Matrix &b) const {
Matrix c;
for (int k = 1; k <= n; k++) {
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
c[i][j] = (c[i][j] + 1ll * a[i][k] * b[k][j]) % P;
}
}
}
return c;
}
Matrix operator+(const Matrix &b) const {
Matrix c;
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
c[i][j] = (a[i][j] + b[i][j]) % P;
}
}
return c;
}
};
struct Node {
Matrix sum, x, y;
Node() : sum(), x(n), y(n) {}
Node operator*(Node b) {
Node a = *this, c;
c.sum = a.sum + a.x * b.sum * a.y;
c.x = a.x * b.x;
c.y = a.y * b.y;
return c;
}
};
Node pow(Node a, int b) {
Node ans;
while (b) {
if (b & 1) ans = ans * a;
a = a * a;
b >>= 1;
}
return ans;
}
Node euclid(long long p, long long q, long long r, long long n, Node U, Node R) {
if (!n) return Node();
if (r >= q) return pow(U, r / q) * euclid(p, q, r % q, n, U, R);
if (p >= q) return euclid(p % q, q, r, n, U, pow(U, p / q) * R);
long long m = ((__int128_t) p * n + r) / q;
if (!m) return pow(R, n);
return pow(R, (q - r - 1) / p) * U * euclid(q, p, (q - r - 1) % p, m - 1, R, U)
* pow(R, n - ((__int128_t) q * m - r - 1) / p);
}
long long p, q, r, l;
int main() {
scanf("%lld%lld%lld%lld%d", &p, &q, &r, &l, &n);
Node U, R;
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
scanf("%d", &R.x[i][j]);
R.sum[i][j] = R.x[i][j];
}
}
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
scanf("%d", &U.y[i][j]);
}
}
auto ans = euclid(p, q, r, l, U, R);
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
printf("%d ", ans.sum[i][j]);
}
printf("\n");
}
return 0;
}