万能欧几里得 学习笔记
模拟赛又双叒叕被卡科技了，学个技能防身/yun（tzc 附体）
考虑这样的问题：将平面内所有 \(x=c\) 或 \(y=c\)（其中 \(c\in\Z\)）的直线标记出来。考虑一条直线 \(y=\dfrac{px+r}q\)（其中 \(p,r\in\N,q\in\N_+\)，而作为 \(r\)，模 \(q\) 的同余类显然是等价的，所以还可以保证 \(r<q\)），令这条直线上的一个动点在横坐标区间 \((0,L]\)（\(L\in\N\)）内从左往右移动，每次遇到竖线就写下一个字符 R，遇到横线就写 U（如果同时遇到，即遇到整点，规定先 U 后 R），最终得到一个 RU 串。再考虑一个幺半群（你可以自己定义）\(\mathbb G\)，里面有两个给定的特殊的元素 \(R,U\in\mathbb G\)，希望快速求出这个 RU 串的值。事实证明，很多类似 \(\sum\limits_{i=1}^Lf\!\left(i,\left\lfloor\dfrac{pi+r}{q}\right\rfloor\right)\) 的和式求值问题都可以归约到这个问题。
思路是类似欧几里得算法的递归：若 \(p\geq q\) 则令 \(p\gets p\bmod q\)，否则交换 \(p,q\)。\(p\geq q\) 的情况（相对）比较简单。注意到当横坐标增大 \(1\) 个单位时，纵坐标增大 \(\dfrac pq\) 个单位，在这个矩形内显然恰好会有一个 R，在该 R 之前会有 \(\left\lfloor\dfrac pq\right\rfloor\sim\left\lfloor\dfrac pq\right\rfloor+1\) 个 U。于是我们考虑将这至少的 \(\left\lfloor\dfrac pq\right\rfloor\) 个 U 砍掉，那么就相当于白送了这么多个 U 然后把纵坐标的步长改为 \(\left\{\dfrac pq\right\}\)（表示小数部分）。这也就相当于令 \(R\gets L^{\left\lfloor\frac pq\right\rfloor}R\) 并令 \(p\gets p\bmod q\)（想一想，为什么）。
接下来考虑第二部分，交换 \(p,q\)、\(R,U\) 的地位（细节比较毒瘤）。如果直接关于 \(y=x\) 翻转坐标系，那么会出现三个问题：
- 由于 \(r\geq 0\)，我们递归到的问题的直线在第一象限中必须要与纵轴有交点，而直接翻转会变成和横轴有交点。
- 翻转之后 \(L\) 不是整数了。
- 原本遇到整点先
U后R，翻转之后就变成先R后U了。
我们一个个来解决这些问题。先看第一个，考虑找到第一个遇到的 U（也就是交于直线 \(y=1\) 处），先算上在这之前的贡献，然后砍掉纵坐标区间 \((0,1]\)，再翻转上面的区域，截距就为正了。前面的贡献是多少？此处我们算出第 \(b\) 个 U 之前有多少个 R，一劳永逸。
答案显然是最大的整数 \(x\) 使 \(\left\lfloor\dfrac{px+r}q\right\rfloor<b\)。由于 \(b\) 是整数，即 \(\dfrac{px+r}q<b\)，\(x<\dfrac{qb-r}p\)，即 \(x\leq\dfrac{qb-r-1}p\)。所以答案就是 \(\left\lfloor\dfrac{qb-r-1}p\right\rfloor\)。
继续上面的讨论。带入 \(b=1\)，前面的贡献显然是 \(R^{\left\lfloor\frac{q-r-1}p\right\rfloor}U\)（由于 \(p<q\)，显然不可能在此处遇到整点，故无需担心后面还有一个 \(R\)）。而砍掉下面翻转上面后，新的直线解析式的 \(p,q\) 交换，截距是多少？就是原直线在 \(y=1\) 的交点横坐标，容易解出为 \(\dfrac{q-r}p\)，故新的 \(r\) 为 \(q-r\)。
考虑第二个问题。解决方案是找到最后一个 U，对其下方的区域进行递归，上面的贡献特殊考虑。显然一共有 \(L_0=\left\lfloor\dfrac{pL+r}{q}\right\rfloor\) 个 U（此时其实注意到解决第一个问题时没有考虑无 U 的情况，现在可以通过 \(L_0=0\) 判断，如果确实无 U 那显然返回 \(R^L\)），并且之前砍掉了一个 U，所以新的 \(L\) 等于 \(L_0-1\)。上面的贡献是多少？显然只有 R，总共有 \(L\) 个 R，减去最后一个 U 前面的 R 数量 \(\left\lfloor\dfrac{qL_0-r-1}p\right\rfloor\)（通过带入 \(b=L_0\) 得到），所以贡献是 \(R^{L-\left\lfloor\frac{qL_0-r-1}{p}\right\rfloor}\)。
最后一个问题：如何处理翻转之后整点处顺序变了？方法非常巧妙：只要令 \(r\gets r-1\) 即可。这样直线往下移了一点点，原来的整点 \(px\equiv 0\pmod q\) 处确实变成先 R 后 U 了，而原来的 \(px\equiv 1\pmod q\) 变成了现在的整点，但人家本来就是先 U 后 R。所以新的 \(r\) 其实等于旧的 \(q-r-1\)，最后为了保证 \(r\) 在模意义范围内，还要 \(\bmod p\)。
到此为止这个算法就讲完啦，分析一下复杂度。跟 gcd 一样，每两次令 \(p\gets p\bmod q\)，分 \(q\leq\dfrac p2\) 与否两种情况讨论发现都有 \(p\) 至少减半，所以递归次数是 \(\mathrm O(\log(p+q))\)。但里面有快速幂，这样就是 2log 吗？错。首先 \(L\) 不增，看作常数。那么每次快速幂的指数都是 \(\mathrm O\!\left(\dfrac qp\right)\)，其 log 便为 \(\mathrm O(\log q-\log p)\)。下一次快速幂时，\(q\) 变为了上次的 \(p\)，\(p\) 变为了更小的 \(q'\)，于是总复杂度加起来就是 \(\log q-\log p+\log p-\log q'+\log q'-\cdots=\mathrm O(\log q)\)，还是一个 log 的（很巧妙啊很巧妙）。所以总复杂度就是 \(\mathrm O(T(\log(p+q)+\log L))\)，其中 \(T\) 是幺半群 \(\mathbb G\) 一次乘法的复杂度。
code
opr euclid(int p, int q, int r, int L, opr U, opr R) {
if(p >= q) {
return euclid(p % q, q, r, L, U, (U ^ p / q) * R);
} else {
int L0 = (p * L + r) / q;
if(!L0) return R ^ L;
opr res = (R ^ ((q - r - 1) / p)) * U;
res = res * euclid(q, p, (q - r - 1) % p, L0 - 1, R, U);
res = res * (R ^ (L - (L0 * q - r - 1) / p));
return res;
}
}
例题：P5170 - 【模板】类欧几里得算法
用万欧水类欧，好！
我们只需要自定义幺半群 \(\mathbb G\)，其它的交给万欧板子就好。
这个地方的标准做法应该是开一个结构体记录答案和转移要用的 \(\Delta\)，只需要 6 个变量，但是我做的时候没想到，使用了最无脑的做法——记录操作，使用了 11 个变量（好劣啊）……那就讲一下我的方法吧，反正标准方法在下一道例题里会提到。
考虑记录一个结构体，里面描述一个操作集合。设当前横坐标下取整为 \(I\)，纵坐标下取整为 \(V\)，三问的答案分别为 \(A,B,C\)，考虑记录 \(A,B,C\) 分别加上关于 \(I,V\) 的一个什么东西，然后 \(I,V\) 分别加上多少。那么是加上什么东西呢？初始情况（即 \(R,U\) 内部）显然是直接加上 \(V,V^2,IV,1,1\)， 但这样合并不满足封闭性。经尝试，如果我们记录 \(A,B,C\) 分别加上的是 \(V\) 的一次函数、\(V\) 的二次函数、\(I,V\) 的多项式，那么是封闭的。于是我们设操作集合为（按顺序地）：
- 令 \(A\gets A+aV+b\)。
- 令 \(B\gets B+cV^2+dV+e\)。
- 令 \(C\gets C+fVI+gV+hI+i\)。
- 令 \(V\gets V+j\)。
- 令 \(I\gets I+k\)。
记录 \(a\sim k\) 这 11 个量。合并就不说了！！因为既简单又复杂（矛 盾 文 学 奖）。\(U\) 是 \(j=1\)，其余量为 \(0\)；\(R\) 是 \(a = c = f = g = k = 1\)，其余量为 \(0\)。注意到题面里 \(i\) 从 \(0\) 开始，所以前两问需要加上 \(i=0\) 的贡献。以及万欧算法需要 \(r<q\)，而初始没有保证，所以要在前面乘上 \(U^{\left\lfloor\frac rq\right\rfloor}\)，然后令 \(r\gets r\bmod q\)。最后计算答案时，由于初始时 \(I=V=0\)，所以只要取常数项 \(b,e,i\)。
code
struct opr {
int a, b;
int c, d, e;
int f, g, h, i;
int j, k;
opr(int o = 0) {
memset(this, 0, sizeof(*this));
if(o == 1) j = 1;
else if(o == 2) a = c = f = g = k = 1;
}
void mod() {
a %= P, b %= P, c %= P, d %= P, e %= P, f %= P, g %= P, h %= P, i %= P, j %= P, k %= P;
}
friend opr operator*(opr x, opr y) {
opr z;
z.a = x.a + y.a;
z.b = x.b + y.a * x.j + y.b;
z.c = x.c + y.c;
z.d = x.d + 2 * y.c * x.j + y.d;
z.e = x.e + y.c * x.j % P * x.j + y.d * x.j + y.e;
z.f = x.f + y.f;
z.g = x.g + y.f * x.k + y.g;
z.h = x.h + y.f * x.j + y.h;
z.i = x.i + y.f * x.j % P * x.k + y.g * x.j + y.h * x.k + y.i;
z.j = x.j + y.j;
z.k = x.k + y.k;
z.mod();
return z;
}
friend opr operator^(opr x, int y) {
opr res;
while(y) {
if(y & 1) res = res * x;
x = x * x;
y >>= 1;
} return res;
}
};
opr euclid(int p, int q, int r, int L, opr U, opr R) {
if(p >= q) {
return euclid(p % q, q, r, L, U, (U ^ p / q) * R);
} else {
int L0 = (p * L + r) / q;
if(!L0) return R ^ L;
opr res = (R ^ ((q - r - 1) / p)) * U;
res = res * euclid(q, p, (q - r - 1) % p, L0 - 1, R, U);
res = res * (R ^ (L - (L0 * q - r - 1) / p));
return res;
}
}
void mian() {
int n = read(), p = read(), r = read(), q = read();
opr res = (opr(1) ^ r / q) * euclid(p, q, r % q, n, opr(1), opr(2));
prt(add(res.b, r / q % P)), pc(' '), prt(add(res.e, (r / q) * (r / q) % P)), pc(' '), prt(res.i), pc('\n');
}
例题：loj6440 - 万能欧几里得
loj 上的模板题。（终于有属于我们万欧自己的模板题啦！！！）
考虑在结构体里维护答案：即维护假装是从一开始就执行该 RU 串得到的答案 \(ans\)。那么考虑两个 RU 串合并，为了算出合并后的 \(ans\) 还需要什么信息。那显然新的 \(ans\) 等于左边串的 \(ans\) 加上经过左边的 RU 作用后的右边串的 \(ans\)。幸运的是，后者是可求的，设左边有 \(i\) 个 \(R\)，\(j\) 个 \(U\)，那么 \(\sum\limits_{i=1}^kA^{f(i)}B^{g(i)}\) 经过这些作用之后就是 \(A^i\!\left(\sum\limits_{d=1}^kA^{f(d)}B^{g(d)}\right)B^j\)。于是我们只要再维护 RU 串执行完毕后的 \(A^i\) 和 \(B^j\) 的值即可，总共维护三个矩阵。
\(R\) 是 \(A^i=A,B^j=I,ans=A\)，\(U\) 是 \(A^i=I,B^j=B,ans=O\)，依然要在之前乘上 \(U^{\left\lfloor\frac rq\right\rfloor}\) 然后令 \(r\gets r\bmod q\)。时间复杂度 \(\mathrm O\!\left(n^3\log\right)\)。有两个细节：一是 ll 乘以 ll 需要 i128（虽然运算结果不会爆，但是过程会爆）；以及 1e9 超过了 998244353 需要手动取模/yun
code
constexpr int N = 30;
int n;
struct matrix {
int a[N][N];
matrix(int x = 0) {
memset(a, 0, sizeof(a));
REP(i, 1, n) a[i][i] = x;
}
int *operator[](int x) { return a[x]; }
friend matrix operator*(matrix f, matrix g) {
matrix h;
REP(i, 1, n) REP(j, 1, n) REP(k, 1, n) addto(h[i][j], f[i][k] * g[k][j] % P);
return h;
}
friend matrix operator+(matrix f, matrix g) {
matrix h;
REP(i, 1, n) REP(j, 1, n) h[i][j] = add(f[i][j], g[i][j]);
return h;
}
void read() {
REP(i, 1, n) REP(j, 1, n) a[i][j] = ::read();
}
void prt() {
REP(i, 1, n) REP(j, 1, n) ::prt(a[i][j]), pc(" \n"[j == n]);
}
} A, B;
struct opr {
matrix ans, Ai, Bj;
opr(int o = 0) {
ans = matrix(), Ai = Bj = matrix(1);
if(o == 1) Bj = B;
else if(o == 2) Ai = ans = A;
}
friend opr operator*(opr x, opr y) {
opr z;
z.ans = x.ans + x.Ai * y.ans * x.Bj;
z.Ai = x.Ai * y.Ai;
z.Bj = x.Bj * y.Bj;
return z;
}
friend opr operator^(opr x, int y) {
opr res;
while(y) {
if(y & 1) res = res * x;
x = x * x;
y >>= 1;
} return res;
}
};
opr euclid(int p, int q, int r, int L, opr U, opr R) {
// debug(mt(p, q, r, L));
if(p >= q) {
return euclid(p % q, q, r, L, U, (U ^ p / q) * R);
} else {
int L0 = ((i128)p * L + r) / q;
if(!L0) return R ^ L;
opr res = (R ^ ((q - r - 1) / p)) * U;
res = res * euclid(q, p, (q - r - 1) % p, L0 - 1, R, U);
res = res * (R ^ (L - ((i128)L0 * q - r - 1) / p));
return res;
}
}
void mian() {
int p = read(), q = read(), r = read(), L = read(); n = read();
A.read(), B.read();
opr res = (opr(1) ^ r / q) * euclid(p, q, r % q, L, opr(1), opr(2));
res.ans.prt();
}
例题：loj138 - 类欧几里得算法
万欧，好😍；类欧，不好😫。
这道题的类欧做法好像要推不少式子，然后还要插值，最终复杂度是 \(k^4\)。而用万欧做脑子都不用动，复杂度也一样，常数应该也更小！我们万欧真是太厉害啦😋。
考虑类似洛谷模板题那样维护三个操作，分别表示 \(x\) 加多少，\(f(x)\) 加多少，和答案加多少。答案加的东西怎么定义才能封闭呢？不难发现，定义成 \(x\) 和 \(f(x)\) 的多项式就好，共 \((k_1+1)(k_2+1)\) 项。于是我们维护操作：
- \(ans\gets ans+\sum\limits_{i=0}^{k_1}\sum\limits_{j=0}^{k_2}a_{i,j}\)。
- \(x\gets x+I\)。
- \(f(x)\gets f(x)+V\)。
\(U\) 的 \(V=1\)，\(R\) 不仅 \(I=1\)，还由于 \((x+1)^{k_1}\) 的二项式展开，有 \(a_{i,k_2}=\dbinom{k_1}{i}\)。如何合并？只讲对系数矩阵 \(a\) 的合并。肯定要累加上左元的矩阵，接下来设 \(I,V\) 为左元的对应项，\(a\) 为右元的 \(a\)，那么 \((x+I)^i(f(x)+V)^j\) 对合并后的 \(a\) 的贡献你可以把这两个多项式二项式展开，然后对应项系数乘起来再乘以 \(a_{i,j}\) 累加到对应位置。这样复杂度就是四方的（\(\mathrm O\!\left(Tk^4\log\right)\)），而且跑不满。
code
constexpr int N = 12;
comb_init(N);
int k1, k2;
struct opr {
int a[N][N], I, V;
int *operator[](int x) { return a[x]; }
opr(int o = 0) {
memset(this, 0, sizeof(*this));
if(o == 1) V = 1;
else if(o == 2) {
I = 1;
REP(i, 0, k1) a[i][k2] = comb[k1][i];
}
}
friend opr operator*(opr f, opr g) {
opr h;
h.I = add(f.I, g.I);
h.V = add(f.V, g.V);
REP(i, 0, k1) REP(j, 0, k2) h[i][j] = f[i][j];
REP(i, 0, k1) REP(j, 0, k2) {
int l[N], r[N];
int now = 1;
REP(k, 0, i) l[i - k] = comb[i][k] * now % P * g[i][j] % P, now = now * f.I % P;
now = 1;
REP(k, 0, j) r[j - k] = comb[j][k] * now % P, now = now * f.V % P;
REP(k, 0, i) REP(o, 0, j) addto(h[k][o], l[k] * r[o] % P);
}
return h;
}
friend opr operator^(opr x, int y) {
opr res;
while(y) {
if(y & 1) res = res * x;
x = x * x;
y >>= 1;
} return res;
}
};
opr euclid(int p, int q, int r, int L, opr U, opr R) {
// debug(mt(p, q, r, L));
if(p >= q) {
return euclid(p % q, q, r, L, U, (U ^ p / q) * R);
} else {
int L0 = (p * L + r) / q;
if(!L0) return R ^ L;
opr res = (R ^ ((q - r - 1) / p)) * U;
res = res * euclid(q, p, (q - r - 1) % p, L0 - 1, R, U);
res = res * (R ^ (L - (L0 * q - r - 1) / p));
return res;
}
}
void mian() {
int L = read(), p = read(), r = read(), q = read(); k1 = read(), k2 = read();
opr res = (opr(1) ^ r / q) * euclid(p, q, r % q, L, opr(1), opr(2));
int ans = res[0][0];
addto(ans, qpow(0, k1) * qpow(r / q, k2));
prt(ans), pc('\n');
}
