万能欧几里得算法

问题

有一条直线 \(y=\frac{Px+K}{Q}\)，其中 \(P\ge 0\) 且 \(0\le K<Q\)。有两种操作 \(U,R\)，从左到右扫描这条直线在 \((0,L]\) 中的部分，若直线跨越了 \(y=k(k\in \mathbb{Z})\) 这条直线，则执行 \(U\) 操作；若直线跨越了 \(x=k(k\in \mathbb{Z})\) 这条直线，则执行 \(R\) 操作。若同时跨越则先执行 \(U\) 再执行 \(R\)。

操作满足结合律，且合并两个操作的时间复杂度是 \(O(1)\)。

万能欧几里得算法

我们设 \(\operatorname{solve}(P,Q,K,L,U,R)\) 为上述问题的答案。类似于欧几里得算法，我们讨论 \(P\) 和 \(Q\) 的大小关系：

• 若 \(P\ge Q\)，那么每个 \(R\) 操作前一定有至少 \(\lfloor \frac{P}{Q}\rfloor\) 个 \(U\) 操作。可以把 \(R\) 操作和这些 \(U\) 操作合并起来，递归到 \(\operatorname{solve}(P\bmod Q,Q,K,L,U,RU^{\lfloor \frac{P}{Q}\rfloor})\)。

• 若 \(P<Q\)，那么需要尝试递归到 \(\operatorname{solve}(Q,P,\dots)\) 之类的东西。原先的问题是，第 \(x\) 个 \(R\) 前面有 \(\lfloor \frac{Px+K}{Q}\rfloor\) 个 \(U\)，那么可以算出第 \(y\) 个 \(U\) 前有 \(\lfloor \frac{Qy-K-1}{P}\rfloor\) 个 \(R\)。

  但直接递归下去不满足 \(0\le K<Q\)。考虑把第一个 \(U\) 及以前的部分截去，这样 \(K\) 就会变成 \(Q-K-1\)。而第一个 \(U\) 前有 \(\lfloor \frac{Q-K-1}{P}\rfloor\) 个 \(R\)，于是递归到 \(R^{\lfloor \frac{Q-K-1}{P}\rfloor}U\operatorname{solve}(Q,P,Q-K-1,C,R,U)\)，其中 \(C=\lfloor\frac{PL+K}{Q}\rfloor\) 是 \(U\) 的个数。

  在最后一个 \(U\) 后面仍有一些 \(R\)，可以算出这样的 \(R\) 的个数是 \(L-\lfloor \frac{QC-K-1}{P}\rfloor\)，再乘上一些 \(R\) 即可。

尽管这个过程中需要做快速幂，但由于快速幂的指数可以估算为 \(\frac{Q}{P}\)，而 \(\log(\frac{Q}{P})=\log Q-\log P\)，递归到下一层后 \(P,Q\) 交换，所以这个 \(\log\) 可以被抵消掉。于是时间复杂度为 \(O(\log \max(P,Q))\)。

代码

ll Div(ll a, ll b, ll c, ll d) {
	return (__int128(a) * b + c) / d;
}
template<typename Info>
Info _Euclid(ll p, ll q, ll k, ll xlim, const Info& u, const Info& r) {
	if (!xlim) return Info();
	if (p >= q) return _Euclid(p % q, q, k, xlim, u, u * (p / q) + r);
	ll m = Div(xlim, p, k, q);
	if (!m) return r * xlim;
	ll cnt = xlim - Div(q, m, -k - 1, p);
	return r * ((q - k - 1) / p) + u + _Euclid(q, p, (q - k - 1) % p, m - 1, r, u) + r * cnt;
}
template<typename Info>
Info Euclid(ll p, ll q, ll k, ll xlim, const Info& u, const Info& r) {
	assert(0 <= p && 0 < q && 0 <= k && 0 <= xlim);
	return u * (k / q) + _Euclid(p, q, k % q, xlim, u, r);
}