【模板】万能欧几里得算法

不等式与整值函数

当 \(a, b, c\in \mathbb N, c\neq 0\) 时，有

\[\left\lfloor\frac{a}{c}\right\rfloor=\frac{a-a\bmod c}{c} \]

这是定义吧。

\[\left\lfloor\frac{a}{c}\right\rfloor=\left\lceil\frac{a+1}{c}\right\rceil-1 \]

\[\left\lceil\frac{a}{c}\right\rceil=\left\lfloor\frac{a-1}{c}\right\rfloor+1 \]

他们说是鸽巢原理的推论。

\[a<b\iff a+1\leq b \]

注意我已经说了 \(a, b\) 是整数。所以我们发现：

\[a<\frac{b}{c}\iff a<\left\lceil \frac{b}{c}\right\rceil\iff a\leq \left\lfloor \frac{b-1}{c}\right\rfloor \]

\[a>\frac{b}{c}\iff a>\left\lfloor \frac{b}{c}\right\rfloor \iff a\geq \left\lceil \frac{b+1}{c}\right\rceil \]

注意这三条都是 iff 关系。

万能欧几里得算法

怎么做

\(solve(p, q, r, l, U, R)\) 表示对直线 \(y=\dfrac{px+r}{q}\ \left(x\in(0, l], 0\leq r<q\right)\) 上，每次碰到横线执行 \(U\) 操作，碰到竖线执行 \(R\) 操作，规定整点上先 \(U\) 再 \(R\)。

第一部分：\(p\geq q\) 时，将问题分解为 \(p\bmod q\) 的子问题。

我们感性理解一下怎么做：令 \(k=\left\lfloor p/q\right\rfloor\)，则

\[y=\frac{(p-kq)x+kqx+r}{q}=\frac{(p-kq)x+r}{q}+kx \]

意思是说每次 \(R\) 之前都有 \(k\) 次额外的 \(U\)（\(+kx\) 部分），所以这个问题 \(=solve(p\bmod q, q, r, l, U, U^{\left\lfloor p/q\right\rfloor}R)\)。

第二部分：\(p<q\) 时，目标是交换 \(p, q\) 与 \(U, R\) 的地位。

这一部分的核心是一个交换推导。因为第 \(a\) 个 \(R\) 之前有 \(\left\lfloor (pa+r)/q\right\rfloor\) 个 \(U\)，不妨设第 \(b\) 个 \(U\) 将在第 \(a\) 个 \(R\) 之后，有：

\[b>\left\lfloor \frac{pa+r}{q}\right\rfloor\implies b>\frac{pa+r}{q}\implies qb>pa+r\implies a<\frac{qb-r}{p}\implies a\leq \left\lfloor \frac{qb-r-1}{p}\right\rfloor \]

也就是说我们成功的交换了 \(U, R\) 的地位，现在第 \(a\leq \left\lfloor (qb-r-1)/p\right\rfloor\) 个 \(R\) 将在第 \(b\) 个 \(U\) 之前了，意思是第 \(b\) 个 \(U\) 之前有 \(\left\lfloor (qb-r-1)/p\right\rfloor\) 个 \(R\)。也就是说我们可以改成求直线 \(y=(qb-r-1)/p\) 的信息了。记 \(R'=U,U'=R\)，现在就是有第 \(b\) 个 \(R'\) 之前有 \(\left\lfloor (qb-r-1)/p\right\rfloor\) 个 \(U'\)。

（这一步有另外一个理解是说你将坐标轴关于 \(y=x\) 对称变为求直线 \(y=(qx-r)/p\) 的信息，但是因为原来要求整点上先 \(U\) 后 \(R\) 现在无法满足，于是将直线向下移动一点点改为 \(y=(qx-r-1)/p\)。至于为什么 \(-1\) 看上面的证明。）

但是我们会遇到一个惊人的负数 which we don't like 并且不符合 \(0\leq r\) 的条件。因为原来有 \(0\leq r<q\)，考虑加一个 \(q\)，写成 \(y=(q(x-1)+q-r-1)/p\)，至于 \(x=1\) 单独计算，它前面有 \(\left\lfloor(q-r-1)/p\right\rfloor\) 个 \(U'\)。

此外还有边界问题，原来相当于是说有 \(l\) 个 \(R\) 操作，现在变成有 \(cntU\) 个 \(R'=U\)，仔细考察，\(cntU=(pl+r)/q\)，然后最后一个 \(R'\) 后面的操作也要算进去，就是 \(l-\left\lfloor (q\cdot cntU-r-1)/p\right\rfloor\) 个 \(U'\)。

也即我们返回

\[R^{\left\lfloor(q-r-1)/p\right\rfloor}U\cdot solve(q, q-r-1, p, cntU-1, R, U)\cdot R^{l-\left\lfloor (q\cdot cntU-r-1)/p\right\rfloor} \]

即可。\(cntU=(pl+r)/q\) 可能需要巨大整数运算。

注意这里有一个细节：\(q-r-1\geq p\) 时，不符合递归的要求，需要修正，修正方法是直接将它 \(\bmod p\)，不需要其它额外贡献，因为根据 \(\text{solve}\) 的定义它只需要管直线就行了。但在一开始进入递归时 \(r\geq q\) 则不仅需要取模还需要前插 \(U^{\left\lfloor r/q\right\rfloor}\)，这是因为 \(\text{solve}\) 只管直线不管下面的点。

怎么写

\[\text{solve}(p, q, r, l, U, R) =\begin{cases} \varnothing, \text{when}\ l=0; \\ \text{solve}(p\bmod q, q, r, l, U, U^{\left\lfloor p/q\right\rfloor}R), \text{when}\ p\geq q;\\ \text{let}\ m=\left\lfloor\dfrac{pl+r}{q}\right\rfloor\ \text{in}\begin{cases} R^l, \text{when}\ m=0;\\ R^{\left\lfloor(q-r-1)/p\right\rfloor}U\cdot \text{solve}(q, (q-r-1)\bmod p, p, m-1, R, U)\cdot R^{l-\left\lfloor(qm-r-1)/p\right\rfloor}, \text{otherwise}. \end{cases}, \text{otherwise}. \end{cases} \]

这里 \(m\) 就表示 \(cntU\)。然后 \(U, R\) 的操作，就是和矩阵有关的维护信息的那一套理论了。

洛谷 P5170 【模板】类欧几里得算法

#include <bits/stdc++.h>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define endl "\n"
#define debug(...) void(0)
#endif
typedef long long LL;
template <class T>
using must_int = enable_if_t<is_integral<T>::value, int>;
template <unsigned umod>
struct modint {
  static constexpr int mod = umod;
  unsigned v;
  modint() : v(0) {}
  template <class T, must_int<T> = 0>
  modint(T x) {
    x %= mod;
    v = x < 0 ? x + mod : x;
  }
  modint operator+() const { return *this; }
  modint operator-() const { return modint() - *this; }
  friend int raw(const modint &self) { return self.v; }
  friend ostream &operator<<(ostream &os, const modint &self) {
    return os << raw(self);
  }
  modint &operator+=(const modint &rhs) {
    v += rhs.v;
    if (v >= umod) v -= umod;
    return *this;
  }
  modint &operator-=(const modint &rhs) {
    v -= rhs.v;
    if (v >= umod) v += umod;
    return *this;
  }
  modint &operator*=(const modint &rhs) {
    v = 1ull * v * rhs.v % umod;
    return *this;
  }
  modint &operator/=(const modint &rhs) {
    assert(rhs.v);
    return *this *= qpow(rhs, mod - 2);
  }
  template <class T, must_int<T> = 0>
  friend modint qpow(modint a, T b) {
    modint r = 1;
    for (; b; b >>= 1, a *= a)
      if (b & 1) r *= a;
    return r;
  }
  friend modint operator+(modint lhs, const modint &rhs) { return lhs += rhs; }
  friend modint operator-(modint lhs, const modint &rhs) { return lhs -= rhs; }
  friend modint operator*(modint lhs, const modint &rhs) { return lhs *= rhs; }
  friend modint operator/(modint lhs, const modint &rhs) { return lhs /= rhs; }
  bool operator==(const modint &rhs) const { return v == rhs.v; }
  bool operator!=(const modint &rhs) const { return v != rhs.v; }
};
typedef modint<998244353> mint;
struct node {
  mint cu, cr, sr, su, suu, sur;
  friend node operator*(const node &S, const node &T) {
    return node{
        S.cu + T.cu, S.cr + T.cr, S.sr + T.sr + S.cr * T.cr,
        S.su + T.su + S.cu * T.cr,
        S.suu + T.suu + T.cr * S.cu * S.cu + 2 * S.cu * T.su,
        // (S.cu + ?)^2 = sum S.cu^2 + sum 2S.cu? + sum?^2
        S.sur + T.sur + T.cr * S.cr * S.cu + S.cu * T.sr + S.cr * T.su
        // (S.cr + ?)(S.cu + %) = sum S.crS.cu +sum?%+?sumS.cu+%sum S.cr
    };
  }
};
node qpow(node a, LL b) {
  node r{};
  for (; b; b >>= 1, a = a * a)
    if (b & 1) r = r * a;
  return r;
}
node solve(LL p, LL q, LL r, LL l, node U, node R) {
  // q > r
  if (!l) return {};
  if (p >= q) return solve(p % q, q, r, l, U, qpow(U, p / q) * R);
  auto div = [](__int128 p, LL q, LL r, LL l) { return (LL)((p * l + r) / q); };
  LL cntu = div(p, q, r, l);
  if (!cntu) return qpow(R, l);
  return qpow(R, (q - r - 1) / p) * U *
         solve(q, p, (q - r - 1) % p, cntu - 1, R, U) *
         qpow(R, l - div(q, p, -r - 1, cntu));
}
int main() {
#ifndef LOCAL
  cin.tie(nullptr)->sync_with_stdio(false);
#endif
  int t;
  cin >> t;
  while (t--) {
    int a, b, c, l;
    cin >> l >> a >> b >> c;
    auto ret = qpow({1, 0, 0, 0, 0, 0}, b / c) *
               solve(a, c, b % c, l, {1, 0, 0, 0, 0, 0}, {0, 1, 1, 0, 0, 0});
    cout << ret.su + b / c << " "
         << ret.suu + 1ll * (b / c) * (b / c) << " " << ret.sur
         << endl;
  }
  return 0;
}