「万能欧几里得」学习笔记

考虑平面直角坐标系上的直线 \(y = \frac{ax+b}c\)，作出所有 \(x = k\) \((k \in \mathbb Z)\) 和所有 \(y = h\) \((h \in \mathbb Z)\)。
从左往右扫描这条直线，当遇到该直线与某水平线的交点，则记录一个 \(\texttt U\)；遇到该直线与某垂直线的交点，则记录一个 \(\texttt R\)。规定整点时水平线的优先级更高。
我们记录的 \(\texttt U, \texttt R\) 实际上可以是任意幺半群信息，当然相同的字母对应相同的元素。

对于参数 \(n,a,b,c\) 和 \(\texttt U,\texttt R\) 执行递归过程 \(f(n,a,b,c,\texttt U,\texttt R)\)：表示在 \(y = \frac{ax+b}c\) \((x \in (0, n])\) 上考虑这个问题。
根据欧几里得过程，考虑讨论：

若 \(a \ge c\)，考虑如何将 \(a\) 取模 \(c\)。
事实上这就是 \(f(n,a \bmod c,b,c,\texttt U,\texttt U^{\lfloor a/c\rfloor}\texttt R)\)。

若 \(a < c\)，考虑如何交换 \(a,c\)。
记 \(m = \left\lfloor\frac{an+b}c\right\rfloor\) 即 \(\texttt U\) 的个数。
如果 \(m = 0\)，以 \(\texttt R^n\) 终止递归。
否则，简单计算得，可以通过在开头补上 \(\texttt R^{\lfloor(c-b-1)/a\rfloor}\texttt U\)，结尾补上 \(\texttt R^{n-\lfloor(cm-b-1)/a\rfloor}\) 的方式，递归到 \(f(m-1,c,b,a,\texttt R,\texttt U)\)。

以下是我擅自从叉姐代码改出来的的一个板子：

#include <algorithm>

using ll = long long;

using namespace std;

namespace UniversalEuclidean {
  template<class MonoidT>
  MonoidT qpow(MonoidT a, ll b) {
    MonoidT ret = MonoidT::identity();
    for (; b; b >>= 1) {
      if (b & 1)
        ret = ret * a;
      a = a * a;
    }
    return ret;
  }
  template<class MonoidT>
  MonoidT sum(ll n, ll a, ll b, ll c) {
    MonoidT r = MonoidT::R();
    MonoidT u = MonoidT::U();
    MonoidT prefix = qpow(u, b / c) * r;
    MonoidT suffix = MonoidT::identity();
    b %= c;
    while (true) {
      if (a >= c)
        r = qpow(u, a / c) * r,
        a %= c;
      else {
        ll m = (a * n + b) / c;
        if (!m)
          return prefix * qpow(r, n) * suffix;
        prefix = prefix * qpow(r, (c - b - 1) / a) * u,
        suffix = qpow(r, n - (c * m - b - 1) / a) * suffix;
        b = (c - b - 1) % a, n = m - 1;
        swap(a, c), swap(u, r);
      }
    }
  }
}

这里是一个 LibreOJ 板子的 AC 代码：

#include <cstdio>
#include <algorithm>

using ll = long long;

using namespace std;

const int mod = 998244353;
inline int norm(int x) {
  return x >= mod ? x - mod : x;
}
inline int reduce(int x) {
  return x < 0 ? x + mod : x;
}
inline int neg(int x) {
  return x ? mod - x : 0;
}
inline void add(int &x, int y) {
  if ((x += y - mod) < 0)
    x += mod;
}
inline void sub(int &x, int y) {
  if ((x -= y) < 0)
    x += mod;
}
inline void fam(int &x, int y, int z) {
  x = (x + (ll)y * z) % mod;
}
inline int qpow(int a, int b) {
  int ret = 1;
  for (; b; b >>= 1) {
    if (b & 1)
      ret = (ll)ret * a % mod;
    a = (ll)a * a % mod;
  }
  return ret;
}

struct Matrix {
  static const int N = 20;
  int a[N][N];
  inline Matrix() {
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        a[i][j] = 0;
  }
  inline int *operator[](int x) {
    return a[x];
  }
  inline const int *operator[](int x) const {
    return a[x];
  }
  inline static Matrix identity() {
    Matrix ret;
    for (int i = 0; i < N; ++i)
      ret[i][i] = 1;
    return ret;
  }
  inline Matrix operator+(const Matrix &o) const {
    Matrix ret;
    for (int i = 0; i < N; ++i)
      for (int j = 0; j < N; ++j)
        ret[i][j] = norm(a[i][j] + o[i][j]);
    return ret;
  }
  inline Matrix operator*(const Matrix &o) const {
    Matrix ret;
    for (int k = 0; k < N; ++k)
      for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
          fam(ret[i][j], a[i][k], o[k][j]);
    return ret;
  }
} A, B, ans;

struct Monoid {
  Matrix AProd, BProd, sum;
  inline Monoid(Matrix a = Matrix(), Matrix b = Matrix(), Matrix c = Matrix()): AProd(a), BProd(b), sum(c) {}
  static Monoid identity() {
    return Monoid(Matrix::identity(), Matrix::identity());
  }
  inline static Monoid R() {
    return Monoid(A, Matrix::identity(), A);
  }
  inline static Monoid U() {
    return Monoid(Matrix::identity(), B);
  }
  inline Monoid operator*(const Monoid &o) const {
    return Monoid(AProd * o.AProd, BProd * o.BProd, sum + AProd * o.sum * BProd);
  }
};

namespace UniversalEuclidean {
  template<class MonoidT>
  MonoidT qpow(MonoidT a, ll b) {
    MonoidT ret = MonoidT::identity();
    for (; b; b >>= 1) {
      if (b & 1)
        ret = ret * a;
      a = a * a;
    }
    return ret;
  }
  template<class MonoidT>
  MonoidT sum(ll n, ll a, ll b, ll c) {
    MonoidT r = MonoidT::R();
    MonoidT u = MonoidT::U();
    MonoidT prefix = qpow(u, b / c);
    MonoidT suffix = MonoidT::identity();
    b %= c;
    while (true) {
      if (a >= c)
        r = qpow(u, a / c) * r,
        a %= c;
      else {
        ll m = ((__int128)a * n + b) / c;
        if (!m)
          return prefix * qpow(r, n) * suffix;
        prefix = prefix * qpow(r, (c - b - 1) / a) * u,
        suffix = qpow(r, n - ((__int128)c * m - b - 1) / a) * suffix;
        b = (c - b - 1) % a, n = m - 1;
        swap(a, c), swap(u, r);
      }
    }
  }
}

ll a, b, c, n;
int m;

int main() {
  scanf("%lld%lld%lld%lld%d", &a, &c, &b, &n, &m);
  for (int i = 0; i < m; ++i)
    for (int j = 0; j < m; ++j)
      scanf("%d", A[i] + j);
  for (int i = 0; i < m; ++i)
    for (int j = 0; j < m; ++j)
      scanf("%d", B[i] + j);
  ans = UniversalEuclidean::sum<Monoid>(n, a, b, c).sum;
  for (int i = 0; i < m; ++i)
    for (int j = 0; j < m; ++j)
      printf("%d%c", ans[i][j], " \n"[j == m - 1]);
}