类欧几里得学习笔记
省选前学不会用到的新知识点多有趣（
1. 普通的类欧几里得
对于 \(p<0\) 或 \(p\ge q\) 或 \(r<0\) 或 \(r\ge q\) 的情况，可以直接把整除的部分提出来，设 \(p'=p\bmod q,r'=r\bmod q\)。
对于其他情况，我们可以改计算 \(y=\frac{px+r}{q}\) 的下面的整点个数变为上面的整点个数。设 \(m=\lfloor\frac{pn+r}{q}\rfloor\)
于是就递归到了更小的情况，时间复杂度 \(O(\log\max\{p,q\})\)。
然后你暴推式子，终于过了 luogu 模板题的时候，你翻到了 LOJ 的模板题。
（你兴奋地骂了句cnm就把这题跳了）
2. 万能欧几里得
实际上有一个有技巧的推柿子方法。
就是对于一条直线 \(y=\frac{px+r}{q}\)，当它经过一个横坐标为整数的点时记一个字符 \(\text{R}\)，经过一个纵坐标为整数的点时记一个字符 \(\text{U}\)，特别的，经过整点的时候先记 \(\text{U}\) 再记 \(\text{R}\)，得到了一个字符串，每个字符表示可能被计算贡献的一部分。我们用一个结构体Node表示一个字符串的贡献。以上面的第一道题为例，设 \(cnt1,cnt2\) 为 \(\text{U}\) 和 \(\text{R}\) 的个数，\(sum\) 表示当前的答案。
我们定义两个字符串的拼接为两个Node的乘法。这个乘法必须满足结合律，但不需要满足交换律（实际上也不太可能）。
（不太理解的话可以画画图）
于是 \(\text{U}=(1,0,0),\text{R}=(0,1,0)\)。我们用函数 \(\texttt{calc(p,q,r,n,U,R)}\) 计算 \(y=\frac{px+r}{q}(x\in (0,n])\) 形成的字符串，往上/右的贡献分别为 \(\text{U,R}\) 时的答案。规定 \(0\le r<q\)。
若 \(p\ge q\) 则返回 \(\texttt{calc(p%q,q,r,n,U,U^(p/q)*R)}\)，因为此时两个 \(\text{R}\) 之间以及第一个 \(\text{R}\) 之前至少有 \(\lfloor\frac{p}{q}\rfloor\) 个 \(\text{U}\)，所以可以进行合并。
设 \(m=\lfloor\frac{pn+r}{q}\rfloor\)。若 \(m=0\) 则返回 \(\text{R}^n\)，因为此时没有 \(\text{U}\)。
其他情况，我们应当把直线对 \(y=x\) 做对称，再把 \(\text{U,R}\) 交换，变为 \(y=\frac{qx-r}{p}(x\in(0,m])\)，答案应当是一样的。但此时到整点时先记 \(\text{R}\) 再记 \(\text{U}\)，所以可以把直线向下平移得到 \(y=\frac{qx-r-1}{p}\)，就变为了先记 \(\text{U}\) 再记 \(\text{R}\)。
但是 \(r'=-r-1\) 又不满足 \(0\le r<q\) 的限制了，所以可以再把直线往左平移得到 \(y=\frac{qx+q-r-1}{p}\)，再往下平移得到 \(y=\frac{qx+((q-r-1)\bmod p)}{p}\)，这就满足了限制。特判 \((0,1]\) 和 \((m,\frac{pn+r}{q}]\) 的贡献。于是返回 \(\texttt{U^{(q-r-1)/p}*R*calc(q,p,(q-r-1)%p,m-1,U,R)*U^{n-(m*q-r-1)/p}}\)。
（其实往上平移和往左平移是一样的，只是为了方便理解，这一部分可能有点绕，建议自己画画图理解一下）
关于时间复杂度，这里用到了Node 的快速幂，但其实时间复杂度不超过 \(O(\log\frac{q}{p})=O(\log q-\log p)\)。所以总时间复杂度为 \(O(F\log\max\{p,q\})\)，其中 \(O(F)\) 表示 Node 乘法的时间复杂度。
3. [LOJ6440] 万能欧几里得
题目描述：求
其中 \(A,B\) 为 \(n\times n\) 的矩阵。
数据范围：\(n\le 20,p,q,r,l,\lfloor\frac{pl}{q}\rfloor\le 10^{18},0\le A_{i,j},B_{i,j}<998244353\)。
这个用普通方法是做不了的，需要用上面的高级方法。
设Node结构体为矩阵三元组 \((X,Y,ans)\)，表示 \(A^{cnt1},B^{cnt2}\) 和答案。于是乘法就是
然后把上面的做法直接套用就完事了。
#include<bits/stdc++.h>
#define Rint register int
using namespace std;
typedef long long LL;
typedef __int128 LLL;
const int mod = 998244353;
template<typename T>
inline void read(T &x){
int ch = getchar(); x = 0;
for(;ch < '0' || ch > '9';ch = getchar());
for(;ch >= '0' && ch <= '9';ch = getchar()) x = x * 10 + ch - '0';
}
inline void qmo(int &x){x += (x >> 31) & mod;}
LL p, q, r, l, n;
struct Mat {
int a[20][20];
Mat(){memset(a, 0, sizeof a);}
Mat operator = (const Mat &o){memcpy(a, o.a, sizeof a); return *this;}
Mat operator + (const Mat &o) const {
Mat res;
for(Rint i = 0;i < n;++ i)
for(Rint j = 0;j < n;++ j)
qmo(res.a[i][j] = a[i][j] + o.a[i][j] - mod);
return res;
}
Mat operator * (const Mat &o) const {
Mat res;
for(Rint i = 0;i < n;++ i)
for(Rint k = 0;k < n;++ k)
for(Rint j = 0;j < n;++ j)
qmo(res.a[i][j] += (LL) a[i][k] * o.a[k][j] % mod - mod);
return res;
}
void print(){
for(Rint i = 0;i < n;++ i)
for(Rint j = 0;j < n;++ j) printf("%d%c", a[i][j], " \n"[j == n - 1]);
}
} A, B, I;
struct Node {
Mat x, y, ans;
Node(const Mat &a = I, const Mat &b = I, const Mat &c = Mat()): x(a), y(b), ans(c){}
Node operator * (const Node &o) const {return Node(x * o.x, y * o.y, ans + x * o.ans * y);}
void print(){ans.print();}
} aa, bb;
Mat ksm(Mat a, LL b){
Mat res = I;
for(;b;b >>= 1, a = a * a) if(b & 1) res = res * a;
return res;
}
Node ksm(Node a, LL b){
Node res;
for(;b;b >>= 1, a = a * a) if(b & 1) res = res * a;
return res;
}
Node calc(LL p, LL q, LL r, LL n, const Node &a, const Node &b){
LL m = ((LLL) p * n + r) / q;
if(!m) return ksm(b, n);
if(p >= q) return calc(p % q, q, r, n, a, ksm(a, p / q) * b);
return ksm(b, (q - r - 1) / p) * a * calc(q, p, (q - r - 1) % p, m - 1, b, a) * ksm(b, n - ((LLL) m * q - r - 1) / p);
}
int main(){
read(p); read(q); read(r); read(l); read(n);
for(Rint i = 0;i < n;++ i) I.a[i][i] = 1;
for(Rint i = 0;i < n;++ i)
for(Rint j = 0;j < n;++ j)
read(A.a[i][j]);
for(Rint i = 0;i < n;++ i)
for(Rint j = 0;j < n;++ j)
read(B.a[i][j]);
aa.x = I; aa.y = B; aa.ans = Mat();
bb.x = A; bb.y = I; bb.ans = A;
(Node(I, ksm(B, r / q), Mat()) * calc(p, q, r % q, l, aa, bb)).print();
}
