【学习笔记】万能欧几里得算法

#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;
}