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

发表于 | 分类于 |

考虑平面直角坐标系上的直线 \(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)\)。

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#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 代码： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#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]);
}

Gitalking ...