省选前学不会用到的新知识点多有趣（
1. 普通的类欧几里得
对于 (p<0) 或 (pge q) 或 (r<0) 或 (rge q) 的情况，可以直接把整除的部分提出来，设 (p'=pmod q,r'=rmod q)。
对于其他情况，我们可以改计算 (y=frac{px+r}{q}) 的下面的整点个数变为上面的整点个数。设 (m=lfloorfrac{pn+r}{q} floor)
于是就递归到了更小的情况，时间复杂度 (O(logmax{p,q}))。
然后你暴推式子，终于过了 luogu 模板题的时候，你翻到了 LOJ 的模板题。
（你兴奋地骂了句cnm就把这题跳了）
2. 万能欧几里得
实际上有一个有技巧的推柿子方法。
就是对于一条直线 (y=frac{px+r}{q})，当它经过一个横坐标为整数的点时记一个字符 ( ext{R})，经过一个纵坐标为整数的点时记一个字符 ( ext{U})，特别的，经过整点的时候先记 ( ext{U}) 再记 ( ext{R})，得到了一个字符串，每个字符表示可能被计算贡献的一部分。我们用一个结构体Node表示一个字符串的贡献。以上面的第一道题为例，设 (cnt1,cnt2) 为 ( ext{U}) 和 ( ext{R}) 的个数，(sum) 表示当前的答案。
我们定义两个字符串的拼接为两个Node的乘法。这个乘法必须满足结合律，但不需要满足交换律（实际上也不太可能）。
（不太理解的话可以画画图）
于是 ( ext{U}=(1,0,0), ext{R}=(0,1,0))。我们用函数 ( exttt{calc(p,q,r,n,U,R)}) 计算 (y=frac{px+r}{q}(xin (0,n])) 形成的字符串，往上/右的贡献分别为 ( ext{U,R}) 时的答案。规定 (0le r<q)。
若 (pge q) 则返回 ( exttt{calc(p%q,q,r,n,U,U^(p/q)*R)})，因为此时两个 ( ext{R}) 之间以及第一个 ( ext{R}) 之前至少有 (lfloorfrac{p}{q} floor) 个 ( ext{U})，所以可以进行合并。
设 (m=lfloorfrac{pn+r}{q} floor)。若 (m=0) 则返回 ( ext{R}^n)，因为此时没有 ( ext{U})。
其他情况，我们应当把直线对 (y=x) 做对称，再把 ( ext{U,R}) 交换，变为 (y=frac{qx-r}{p}(xin(0,m]))，答案应当是一样的。但此时到整点时先记 ( ext{R}) 再记 ( ext{U})，所以可以把直线向下平移得到 (y=frac{qx-r-1}{p})，就变为了先记 ( ext{U}) 再记 ( ext{R})。
但是 (r'=-r-1) 又不满足 (0le r<q) 的限制了，所以可以再把直线往左平移得到 (y=frac{qx+q-r-1}{p})，再往下平移得到 (y=frac{qx+((q-r-1)mod p)}{p})，这就满足了限制。特判 ((0,1]) 和 ((m,frac{pn+r}{q}]) 的贡献。于是返回 ( exttt{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(logfrac{q}{p})=O(log q-log p))。所以总时间复杂度为 (O(Flogmax{p,q}))，其中 (O(F)) 表示 Node 乘法的时间复杂度。
3. [LOJ6440] 万能欧几里得
题目描述：求
其中 (A,B) 为 (n imes n) 的矩阵。
数据范围：(nle 20,p,q,r,l,lfloorfrac{pl}{q} floorle 10^{18},0le 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], "
"[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();
}