【算法学习笔记】类欧几里得算法
一个基础的数论问题。
试求 \(\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor\) 的值，其中：\(a,b \ge 0\)，\(n,c >0\)
在Atcoder的AC库中有这样一个函数可以在 \(\mathcal{O}(lg(n + c + a + b))\) 的时间内解决问题。
函数代码 ↓
using ll = long long;
ll floor_sum(ll n, ll m, ll a, ll b) {
ll ans = 0;
if (a >= m) {
ans += (n - 1) * n * (a / m) / 2;
a %= m;
}
if (b >= m) {
ans += n * (b / m);
b %= m;
}
ll y_max = (a * n + b) / m, x_max = (y_max * m - b);
if (y_max == 0) return ans;
ans += (n - (x_max + a - 1) / a) * y_max;
ans += floor_sum(y_max, a, m, (a - x_max % a) % a);
return ans;
}
好奇它的证明过程，然后在 OI wiki 上找到相应文档，这个算法名为：类欧几里德算法
个人证明
另外补上个人证明：
当 \(a \ge c\) 或 \(b\ge c\) 时，
\[f(a,b,c,n) = \frac{n(n + 1)}{2}*\left\lfloor \frac{a}{c} \right\rfloor + (n + 1) * \left\lfloor \frac{b}{c} \right\rfloor + f(a\ mod\ c,b\ mod\ c,c,n). \]当 \(a < c\) 并且 \(b < c\) 时，
设 \(m = \left\lfloor \frac{an+b}{c} \right\rfloor\)，
\[f(a,b,c,n) = mn - f(c,c - b - 1,a,m-1). \]然后递归至 \(a = 0\) 即可.
例题
数形结合, 把式子稍微简单转换一下, 套用类欧几里得算法即可.
引入
设
其中 \(a,b,c,n\) 是常数。需要一个 \(O(\log n)\) 的算法。
这个式子和我们以前见过的式子都长得不太一样。带向下取整的式子容易让人想到数论分块，然而数论分块似乎不适用于这个求和。但是我们是可以做一些预处理的。
如果说 \(a\ge c\) 或者 \(b\ge c\)，意味着可以将 \(a,b\) 对 \(c\) 取模以简化问题：
那么问题转化为了 \(a<c,b<c\) 的情况。观察式子，你发现只有 \(i\) 这一个变量。因此要推就只能从 \(i\) 下手。在推求和式子中有一个常见的技巧，就是条件与贡献的放缩与转化。具体地说，在原式 \(\displaystyle f(a,b,c,n)=\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor\) 中，\(0\le i\le n\) 是条件，而 \(\left\lfloor \dfrac{ai+b}{c} \right\rfloor\) 是对总和的贡献。
要加快一个和式的计算过程，所有的方法都可以归约为 贡献合并计算。但你发现这个式子的贡献难以合并，怎么办？将贡献与条件做转化 得到另一个形式的和式。具体地，我们直接把原式的贡献变成条件：
现在多了一个变量 \(j\)，既然算 \(i\) 的贡献不方便，我们就想办法算 \(j\) 的贡献。因此想办法搞一个和 \(j\) 有关的贡献式。这里有另一个家喻户晓的变换方法，笔者概括为限制转移。具体来说，在上面的和式中 \(n\) 限制 \(i\) 的上界，而 \(i\) 限制 \(j\) 的上界。为了搞 \(j\)，就先把 j 放到贡献的式子里，于是我们交换一下 \(i,j\) 的求和算子，强制用 \(n\) 限制 \(j\) 的上界。
这样做的目的是让 \(j\) 摆脱 \(i\) 的限制，现在 \(i,j\) 都被 \(n\) 限制，而贡献式看上去是一个条件，但是我们仍把它叫作贡献式，再对贡献式做变换后就可以改变 \(i,j\) 的限制关系。于是我们做一些放缩的处理。首先把向下取整的符号拿掉
然后可以做一些变换
最后一步，向下取整得到：
这一步的重要意义在于，我们可以把变量 \(i\) 消掉了！具体地，令 \(m=\left\lfloor \frac{an+b}{c} \right\rfloor\)，那么原式化为
这是一个递归的式子。并且你发现 \(a,c\) 分子分母换了位置，又可以重复上述过程。先取模，再递归。这就是一个辗转相除的过程，这也是类欧几里德算法的得名。
容易发现时间复杂度为 \(\mathcal{O}(lg(n + m + a + b))\)。
同时关于 类欧几里德算法 有两个函数的拓展
扩展
理解了最基础的类欧几里德算法，我们再来思考以下两个变种求和式：
推导 g
我们先考虑 \(g\)，类似地，首先取模：
接下来考虑 \(a<c,b<c\) 的情况，令 \(m=\left\lfloor\frac{an+b}{c}\right\rfloor\)。之后的过程我会写得很简略，因为方法和上文略同：
这时我们设 \(t=\left\lfloor\frac{jc+c-b-1}{a}\right\rfloor\)，可以得到
推导 h
同样的，首先取模：
考虑 \(a<c,b<c\) 的情况，\(m=\left\lfloor\dfrac{an+b}{c}\right\rfloor, t=\left\lfloor\dfrac{jc+c-b-1}{a}\right\rfloor\).
我们发现这个平方不太好处理，于是可以这样把它拆成两部分：
这样做的意义在于，添加变量 \(j\) 的时侯就只会变成一个求和算子，不会出现 \(\sum\times \sum\) 的形式：
接下来考虑化简前一部分：
因此
在代码实现的时侯，因为 \(3\) 个函数各有交错递归，因此可以考虑三个一起整体递归，同步计算，否则有很多项会被多次计算。这样实现的复杂度是 \(O(\log n)\) 的。
模板代码实现
#define int long long using namespace std; const int P = 998244353; int i2 = 499122177, i6 = 166374059; struct data { data() { f = g = h = 0; } int f, g, h; }; // 三个函数打包 data calc(int n, int a, int b, int c) { int ac = a / c, bc = b / c, m = (a * n + b) / c, n1 = n + 1, n21 = n * 2 + 1; data d; if (a == 0) { // 迭代到最底层 d.f = bc * n1 % P; d.g = bc * n % P * n1 % P * i2 % P; d.h = bc * bc % P * n1 % P; return d; } if (a >= c || b >= c) { // 取模 d.f = n * n1 % P * i2 % P * ac % P + bc * n1 % P; d.g = ac * n % P * n1 % P * n21 % P * i6 % P + bc * n % P * n1 % P * i2 % P; d.h = ac * ac % P * n % P * n1 % P * n21 % P * i6 % P + bc * bc % P * n1 % P + ac * bc % P * n % P * n1 % P; d.f %= P, d.g %= P, d.h %= P; data e = calc(n, a % c, b % c, c); // 迭代 d.h += e.h + 2 * bc % P * e.f % P + 2 * ac % P * e.g % P; d.g += e.g, d.f += e.f; d.f %= P, d.g %= P, d.h %= P; return d; } data e = calc(m - 1, c, c - b - 1, a); d.f = n * m % P - e.f, d.f = (d.f % P + P) % P; d.g = m * n % P * n1 % P - e.h - e.f, d.g = (d.g * i2 % P + P) % P; d.h = n * m % P * (m + 1) % P - 2 * e.g - 2 * e.f - d.f; d.h = (d.h % P + P) % P; return d; } int T, n, a, b, c; signed main() { scanf("%lld", &T); while (T--) { scanf("%lld%lld%lld%lld", &n, &a, &b, &c); data ans = calc(n, a, b, c); printf("%lld %lld %lld\n", ans.f, ans.h, ans.g); } return 0; }
