Blog of RuSun

\begin {array}{c} \mathfrak {One Problem Is Difficult} \\\\ \mathfrak {Because You Don't Know} \\\\ \mathfrak {Why It Is Diffucult} \end {array}

万能欧几里得算法模板

对于一条直线 $\frac {a x + b} c$ ,如果穿过整数直线,与 $x$ 平行则执行 $U$ ,与 $y$ 平行则执行 $R$ 。那么可以在 $\log(a + c)$ 时间内完成 $n$ 内的操作。

例题 I

求 $\sum\limits_{i=0}^{n}\lfloor \frac{ai+b}{c} \rfloor, \sum\limits_{i=0}^{n}{\lfloor \frac{ai+b}{c} \rfloor}^2, \sum\limits_{i=0}^{n}i\lfloor \frac{ai+b}{c} \rfloor$ 。

对于第一个问题,对于 $U$ 操作为 ++t ,对于 $R$ 操作为 res += t 。那么可以刻画出这个问题。

考虑结合,两组操作如何合并。记 $f, g, h$ 为三个问题的答案,$u, r$ 表示两种操作的操作的次数。考虑操作 $1$ 对操作 $2$ 的增加的贡献,记 $s _ i$ 表示在 $i$ 处的点值,那么有

$$
\begin {aligned}
f & = \sum _ {i = 0} ^ {r - 1} s _ i \\
g & = \sum _ {i = 0} ^ {r - 1} s _ i ^ 2 \\
h & = \sum _ {i = 0} ^ {r - 1} i s _ i
\end {aligned}
$$

$$
\begin {aligned}
f & \gets f _ 1 + \sum _ {i = 0} ^ {r _ 2 - 1} (s _ i + u _ 1) = f _ 1 + f _ 2 + u _ 1 r _ 2 \\
g & \gets g _ 1 + \sum _ {i = 0} ^ {r _ 2 - 1} (s _ i + u _ 1) ^ 2 = g _ 1 + \sum _ {i = 0} ^ {r _ 2 - 1} (s _ i ^ 2 + 2s _ i u _ 1 + u _ 1 ^ 2) = g _ 1 + g _ 2 + 2 f _ 2 u _ 1 + u _ 1 ^ 2 r _ 2 \\
h & \gets h _ 1 + \sum _ {i = 0} ^ {r _ 2 - 1} (i + r _ 1) (s _ i + u _ 1) = h _ 1 + h _ 2 + \frac {r _ 2 (r _ 2 - 1)} 2u _ 1 + r _ 1 f _ 2 + r _ 2 r _ 1 u _ 1
\end {aligned}
$$

另外注意对于开头的处理。

查看代码
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
#include <cstdio>
#include <algorithm>
using namespace std;
template <class Type>
void read (Type &x)
{
char c;
bool flag = false;
while ((c = getchar()) < '0' || c > '9')
c == '-' && (flag = true);
x = c - '0';
while ((c = getchar()) >= '0' && c <= '9')
x = (x << 1) + (x << 3) + c - '0';
flag && (x = ~x + 1);
}
template <class Type, class ...rest>
void read (Type &x, rest &...y) { read(x), read(y...); }
template <class Type>
void write (Type x)
{
x < 0 && (putchar('-'), x = ~x + 1);
x > 9 && (write(x / 10), 0);
putchar('0' + x % 10);
}
typedef long long LL;
const int mod = 998244353;
struct Node
{
int f, g, h, u, r;
Node (int _f = 0, int _g = 0, int _h = 0, int _u = 0, int _r = 0) : f(_f), g(_g), h(_h), u(_u), r(_r) { }
friend Node operator * (Node a, Node b)
{ return Node((a.f + b.f + (LL)a.u * b.r) % mod, (a.g + b.g + 2ll * b.f * a.u + (LL)a.u * a.u % mod * b.r) % mod, (a.h + b.h + (LL)b.f * a.r + (b.r - 1ll) * b.r / 2 % mod * a.u + (LL)a.r * b.r % mod * a.u) % mod, (a.u + b.u) % mod, (a.r + b.r) % mod); }
};
const Node I, U(0, 0, 0, 1, 0), R(0, 0, 0, 0, 1);
Node binpow (Node b, int k)
{
Node res = I;
for (; k; k >>= 1, b = b * b)
if (k & 1) res = res * b;
return res;
}
Node ugcd(int a, int b, int c, int n, Node U, Node R)
{
b %= c;
if (a >= c) return ugcd(a % c, b, c, n, U, binpow(U, a / c) * R);
LL m = ((LL)a * n + b) / c;
if (m == 0) return binpow(R, n);
swap(U, R);
return binpow(U, (c - b - 1) / a) * R * ugcd(c, c - b - 1, a, m - 1, U, R) * binpow(U, n - (c * m - b - 1) / a);
}
int main ()
{
int T; read(T);
for (int n, a, b, c; T; --T)
{
read(n, a, b, c);
Node res = binpow(U, b / c) * R * ugcd(a, b, c, n, U, R);
write(res.f), putchar(' '), write(res.g), putchar(' '), write(res.h), puts("");
}
return 0;
}

例题II

求 $\sum\limits_{x=1}^LA^xB^{\left\lfloor\frac{Px+R}{Q}\right\rfloor}$ 。其中 $A, B$ 为矩阵。

此题中,$U$ 操作为在后面乘 $B$ ,$R$ 操作为在前面乘 $A$ ,并将答案贡献上。那么两组操作合并,考虑第一组对第二组的影响,即将 $A ^ {k _ A}$ 乘在前面,$B ^ {k _ B}$ 乘在后面。

查看代码
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
#include <cstdio>
#include <algorithm>
using namespace std;
template <class Type>
void read (Type &x)
{
char c;
bool flag = false;
while ((c = getchar()) < '0' || c > '9')
c == '-' && (flag = true);
x = c - '0';
while ((c = getchar()) >= '0' && c <= '9')
x = (x << 1) + (x << 3) + c - '0';
flag && (x = ~x + 1);
}
template <class Type, class ...rest>
void read (Type &x, rest &...y) { read(x), read(y...); }
template <class Type>
void write (Type x)
{
x < 0 && (putchar('-'), x = ~x + 1);
x > 9 && (write(x / 10), 0);
putchar('0' + x % 10);
}
typedef __int128 L;
typedef long long LL;
const int N = 30, mod = 998244353;
int m;
struct Matrix
{
int w[N][N];
auto &operator [] (int k) { return w[k]; }
void print ()
{
for (int i = 0; i < m; ++i, puts(""))
for (int j = 0; j < m; ++j, putchar(' '))
write(w[i][j]);
}
Matrix (bool op = false)
{
for (int i = 0; i < m; ++i)
for (int j = 0; j < m; ++j)
w[i][j] = op && i == j;
}
friend Matrix operator * (Matrix a, Matrix b)
{
Matrix res = Matrix();
for (int i = 0; i < m; ++i)
for (int j = 0; j < m; ++j)
for (int k = 0; k < m; ++k)
res[i][j] = (res[i][j] + (LL)a[i][k] * b[k][j]) % mod;
return res;
}
friend Matrix operator + (Matrix a, Matrix b)
{
Matrix res = Matrix();
for (int i = 0; i < m; ++i)
for (int j = 0; j < m; ++j)
res[i][j] = (a[i][j] + b[i][j]) % mod;
return res;
}
} A, B;
struct Node
{
Matrix f, u, r;
Node (Matrix _f = Matrix(0), Matrix _u = Matrix(1), Matrix _r = Matrix(1)) : f(_f), u(_u), r(_r) { }
friend Node operator * (Node a, Node b)
{ return Node(a.f + a.r * b.f * a.u, a.u * b.u, a.r * b.r); }
};
Node binpow (Node b, LL k)
{
Node res = Node();
for (; k; k >>= 1, b = b * b)
if (k & 1) res = res * b;
return res;
}
Node ugcd(LL a, LL b, LL c, LL n, Node U, Node R)
{
b %= c;
if (a >= c) return ugcd(a % c, b, c, n, U, binpow(U, a / c) * R);
LL m = ((L)a * n + b) / c;
if (m == 0) return binpow(R, n);
swap(U, R);
return binpow(U, (c - b - 1) / a) * R * ugcd(c, c - b - 1, a, m - 1, U, R) * binpow(U, n - ((L)c * m - b - 1) / a);
}
int main ()
{
LL a, b, c, n;
read(a, c, b, n, m);
for (int i = 0; i < m; ++i)
for (int j = 0; j < m; ++j)
read(A[i][j]);
for (int i = 0; i < m; ++i)
for (int j = 0; j < m; ++j)
read(B[i][j]);
Node U(Matrix(), B, Matrix(1)), R(A, Matrix(1), A);
(binpow(U, b / c) * ugcd(a, b, c, n, U, R)).f.print();
return 0;
}