类欧几里得算法
类欧几里得算法
最近读具体数学第二、三章有感，遂来挑战一下这题。这题做完以后感觉对求和技巧与规约的理解又更进了一步。
听说有个叫万能欧几里得算法的东西，是这类题目的通解做法，奈何本人水平太菜，有空再补。
类欧几里得算法要求以下三个式子的值：
\[\begin{aligned}
f(n,a,b,c)&=\sum_{i=0}^n \lfloor\frac{ai+b}{c}\rfloor\\ g(n,a,b,c)&=\sum_{i=0}^n \lfloor\frac{ai+b}{c}\rfloor^2\\ h(n,a,b,c)&=\sum_{i=0}^n i\lfloor\frac{ai+b}{c} \rfloor
\end{aligned}
\]
我们将每个式子设成函数的用处很快就能体现出来，因为类欧几里得算法之所以这么命名，不是因为其推导过程和欧几里得算法有类似之处，而是因为其递归方式和复杂度证明方式与欧几里得算法一致。
算 \(f\)
先考虑最简单的第一个式子如何求，这里说简单只是式子最简单，而其操作方式并不简单，如果明白了如何递归求解 \(f(n, a, b,c)\) 后面两个式子的推导思路就大差不差了。
根据具体数学的启示，我们先观察该式的边界情况，即 \(a=0\) 或 \(b=0\)，发现实际上只有 \(a=0\) 的答案是平凡的，\(b=0\) 也不太好处理。当 \(a=0\) 时，\(f(n,a,b,c) = (n+1)\lfloor \dfrac bc\rfloor\)。
既然 \(a=0\) 是好处理的，那么我们能否像欧几里得算法那样不断递归，将任何 \(a\) 都消成 \(0\) 呢？、
首先考虑 \(a\ge c\) 或者 \(b\ge c\) 的情况，受到具体数学 \(P77\) 对于类似的一个和式的处理的启发，我们能够将 \(a\) 与 \(b\) 拆分成 \(\lfloor \dfrac{a}{c}\rfloor c+ (a\bmod c)\) 与 \(\lfloor \dfrac{b}{c}\rfloor c+ (b\bmod c)\)，使得原式的范围缩小，从而得到规约的效果：
\[\begin{aligned}
f(n,a,b,c) &= \sum_{i=0}^n \lfloor\frac{(a\bmod c)i +\lfloor\dfrac{a}{c}\rfloor ci+(b\bmod c)+\lfloor\dfrac{b}{c}\rfloor c }{c}\rfloor\\
&=\sum_{i=0}^n \lfloor \frac{(a\bmod c)i+(b\bmod c)}{c}\rfloor+\sum_{i=0}^n i\lfloor\frac{a}{c} \rfloor+\sum_{i=0}^n \lfloor\frac{b}{c} \rfloor\\
&=f(n, a\bmod c, b\bmod c,c)+\frac{n(n+1)}{2}\lfloor\frac ac\rfloor+(n+1)\lfloor\frac bc\rfloor
\end{aligned}
\]
不断执行上面的操作只能让 \(a< c\) 且 \(b <c\)，接下来该怎么进一步缩小范围呢？
当 \(a<c\) 且 \(b<c\) 时，我们刚才换 \(a,b\) 的方法显然是行不通了，我们不妨换一种求和的思路，还记得具体数学 \(P73\) 的思想吗？我们不妨用 \(\sum_j [1\le j\le x]\) 来替换 \(\lfloor x\rfloor\)。
\[\begin{aligned}
f(n,a,b,c)&= \sum_{i=0}^n \sum_j [1\le j\le \frac{ai+b}{c}]\\
&\text{接下来换一个量进行枚举}\\
&=\sum_{j=1}^{\lfloor \frac{an+b}c \rfloor} \sum_{i}[i\le n][jc\le ai+b]\\
&\text{用 j+1替换j，使得下标符合f从0开始的形式，同时利用充要条件}n,m\in \Z ,n\le m\Leftrightarrow n<m+1\\
&=\sum_{j=0}^{\lfloor \frac{an+b}c \rfloor-1}\sum_i [i\le n][(j+1)c<ai+b+1]\\
&=\sum_{j=0}^{\lfloor \frac{an+b}c \rfloor-1}\sum_i [\frac{jc+c-b-1}{a}<i\le n]\\
&\text{由于} x<n\Leftrightarrow \lfloor x\rfloor <n \text{并结合具体数学 (3.12)}\\
&=\sum_{j=0}^{\lfloor \frac{an+b}c \rfloor-1}n-\lfloor\frac{jc+c-b-1}{a}\rfloor\\
&=\lfloor\frac{an+b}c\rfloor n-f(\lfloor\frac{an+b}c\rfloor-1, c,c-b-1,a)
\end{aligned}
\]
有两点需要注意：
- 第三步的代换十分巧妙，需要好好体会，因为这么做可以使区间变为左开右闭区间，这样统计区间内整数个数时减去的是一个下取整函数，而这个下取整函数的形式恰好可以规约
- \(b\) 对我们而言作用不大，我们主要想让这个类欧几里得过程最终能够到达 \(a=0\) 的边界情况，而我们这部分的操作相当于将 \(a,c\) 位置互换，从而使得 \(a\) 始终在减小，故而 \(a\) 最终总会到达 \(0\)，其实如果仔细看的话，只观察 \(a,c\) 的变化，就等同于欧几里得算法，根据主定理可以分析复杂度的确是 \(O(\log n)\)
总结来说，
\[f(n, a, b,c)=\begin{cases}
(n+1)\lfloor \dfrac bc\rfloor &a=0\\
f(n, a\bmod c, b\bmod c,c)+\dfrac{n(n+1)}{2}\lfloor\dfrac ac\rfloor+(n+1)\lfloor\dfrac bc\rfloor& a\ge c \ or\ b\ge c\\
\lfloor\dfrac{an+b}c\rfloor n-f(\lfloor\dfrac{an+b}c\rfloor-1, c,c-b-1,a) & otherwise\\
\end{cases}
\]
算 \(g\) 与 \(f\)
为什么 \(g\) 和 \(f\) 我并到一块了呢？看完后面的推导，我们就能够发现这两者是互推的。
先考虑 \(g\) 的求法，同样，先讨论 \(a=0\) 的平凡情况，此时 \(g(n, a, b,c)= (n+1)\lfloor\dfrac{b}{c}\rfloor ^2\)
当 \(a\ge c\) 或者 \(b\ge c\) 时，推导差不多，只需要将平方展开即可：
\[\begin{aligned}
g(n,a,b,c)&= \sum_{i=0}^n \left(\lfloor\frac{(a\bmod c)i +(b\bmod c)}{c}\rfloor+\lfloor\dfrac{a}{c}\rfloor i +\lfloor\dfrac{b}{c}\rfloor\right )^2\\
&= \sum_{i=0}^n \lfloor\frac{(a\bmod c)i +(b\bmod c)}{c}\rfloor^2 +\sum_{i=0}^n\lfloor\dfrac{a}{c}\rfloor^2i^2+\sum_{i=0}^n \lfloor\dfrac{b}{c}\rfloor^2\\& \quad+ 2\lfloor\dfrac{a}{c}\rfloor\sum_{i=0}^ni\lfloor\frac{(a\bmod c)i +(b\bmod c)}{c}\rfloor+2\lfloor\dfrac{a}{c}\rfloor\lfloor\dfrac{b}{c}\rfloor\sum_{i=0}^n
i \\ & \quad + 2\lfloor\dfrac{b}{c}\rfloor\sum_{i=0}^n\lfloor\frac{(a\bmod c)i +(b\bmod c)}{c}\rfloor\\
&=g(n,a\bmod c, b\bmod c, c)+\dfrac{n(n+1)(2n+1)}{6}\lfloor\dfrac{a}{c}\rfloor^2+(n+1)\lfloor\dfrac{b}{c}\rfloor^2\\&\quad+2\lfloor\dfrac{a}{c}\rfloor h(n, a\bmod c, b\bmod c, c)+n(n+1)\lfloor\dfrac{a}{c}\rfloor\lfloor\dfrac{b}{c}\rfloor +2\lfloor\dfrac{b}{c}\rfloor f(n, a\bmod c,b\bmod c, c)
\end{aligned}
\]
当 \(a<c\) 且 \(b<c\) 时，推导需要运用到具体数学 \(P31\) 页的 \((2.33)\) 变换和 \(P39\) 的展开和收缩变换：
\[\begin{aligned}
g(n,a,b,c)&=\sum_{i=0}^n \left(\sum_j [1\le j\le \frac{ai+b}{c}]\right)^2\\
&=\sum_{i=0}^n \left(\sum_j [1\le j\le \frac{ai+b}{c}]\right)\left(\sum_k [1\le k\le \frac{ai+b}{c}]\right)\\
&=\sum_{i=0}^n \sum_{j,k} [1\le j,k\le\frac{ai+b}{c}] \\
&\text{根据(2.33) }[1\le j,k\le n]+[1\le j=k\le n]=[1\le j\le k\le n]+[1\le k \le j\le n]\\
&= \sum_{i=0}^n\left(2\sum_{j,k} [1\le j\le k\le\frac{ai+b}{c} ]- \sum_{j,k}[1\le j=k\le \frac{ai+b}{c}]\right)\\
&\text{剩下的操作和求 f类似，先改变枚举顺序，再用k+1替换k}\\
&= 2\sum_{k=0}^{\lfloor \frac{an+b}c \rfloor - 1} \sum_{j=1}^{k+1}\sum_{i}[\frac{kc+c-b-1}{c}<i\le n]- \sum_{i=0}^ n \lfloor\frac{ai+b}{c}\rfloor\\
&=2\sum_{k=0}^{\lfloor \frac{an+b}c \rfloor - 1} (k+1)\left(n-\lfloor\frac{kc+c-b-1}{c}\rfloor\right)- \sum_{i=0}^ n \lfloor\frac{ai+b}{c}\rfloor\\
&= \lfloor \frac{an+b}c \rfloor\left(\lfloor \frac{an+b}c \rfloor+1\right)n-2h(\lfloor \frac{an+b}c \rfloor-1,c, c-b-1,a)\\&\quad-2f(\lfloor \frac{an+b}c \rfloor-1,c,c-b-1,a)-f(n,a,b,c)
\end{aligned}
\]
最后考虑 \(h\) 的求法，\(a=0\) 时，\(h(n,a,b,c) = \dfrac{n(n+1)}2 \lfloor\dfrac bc\rfloor\)
当 \(a\ge c\) 或 \(b\ge c\) 时，推导仍然比较简单，暴力展开就行：
\[\begin{aligned}
h(n,a,b,c)&=\sum_{i=0}^n i\left(\lfloor\frac{(a\bmod c)i +(b\bmod c)}{c}\rfloor+\lfloor\dfrac{a}{c}\rfloor i +\lfloor\dfrac{b}{c}\rfloor\right )\\
&=h(n, a\bmod c, b\bmod c,c)+\dfrac{n(n+1)(2n+1)}6 \lfloor\dfrac{a}{c}\rfloor+\dfrac{n(n+1)}2 \lfloor\dfrac{b}{c}\rfloor\\
\end{aligned}
\]
当 \(a<c\) 且 \(b<c\) 时，和之前一样：
\[\begin{aligned}
h(n,a,b,c)&=\sum_{i=0}^ni\sum_j[1\le j\le \dfrac{ai+b}c]\\
&= \sum_{j=0}^{\lfloor \frac{an+b}c\rfloor-1}\sum_i i[\frac{jc+c-b-1}{a}<i\le n]\\
&=\sum_{j=0}^{\lfloor \frac{an+b}c\rfloor-1}\left(n+1+\lfloor\frac{jc+c-b-1}{a}\rfloor\right)\left(n-\lfloor\frac{jc+c-b-1}{a}\rfloor\right)\\
&=\frac12\sum_{j=0}^{\lfloor \frac{an+b}c\rfloor-1}(n^2+n)-\frac12\sum_{j=0}^{\lfloor \frac{an+b}c\rfloor-1} \lfloor \frac{jc+c-b-1}a\rfloor^2-\frac12 \sum_{j=0}^{\lfloor \frac{an+b}c\rfloor-1}\lfloor \frac{jc+c-b-1}a\rfloor\\
&=\frac{n(n+1)}2 \lfloor \frac{an+b}c\rfloor-\frac12 g(\lfloor \frac{an+b}c\rfloor-1, c,c-b-1,a)\\&\quad-\frac12 f(\lfloor \frac{an+b}c\rfloor-1, c,c-b-1,a)
\end{aligned}
\]
推导工作全部进行完毕。
\[g(n,a,b,c)=\begin{cases}
(n+1)\lfloor\dfrac{b}{c}\rfloor ^2 &a=0\\
g(n,a\bmod c, b\bmod c, c)+\dfrac{n(n+1)(2n+1)}{6}\lfloor\dfrac{a}{c}\rfloor^2+(n+1)\lfloor\dfrac{b}{c}\rfloor^2\\+2\lfloor\dfrac{a}{c}\rfloor h(n, a\bmod c, b\bmod c, c)+n(n+1)\lfloor\dfrac{a}{c}\rfloor\lfloor\dfrac{b}{c}\rfloor \\+2\lfloor\dfrac{b}{c}\rfloor f(n, a\bmod c,b\bmod c, c) &a\ge c\ or \ b\ge c\\
\lfloor \dfrac{an+b}c \rfloor\left(\lfloor \dfrac{an+b}c \rfloor+1\right)n-2h(\lfloor \dfrac{an+b}c \rfloor-1,c, c-b-1,a)\\-2f(\lfloor \dfrac{an+b}c \rfloor-1,c,c-b-1,a)-f(n,a,b,c) &otherwise
\end{cases}
\]
\[h(n,a,b,c)=\begin{cases}
\dfrac{n(n+1)}2 \lfloor\dfrac bc\rfloor & a=0\\
h(n, a\bmod c, b\bmod c,c)+\dfrac{n(n+1)(2n+1)}6 \lfloor\dfrac{a}{c}\rfloor+\dfrac{n(n+1)}2 \lfloor\dfrac{b}{c}\rfloor &a\ge c\ or \ b\ge c\\
\dfrac{n(n+1)}2 \lfloor \dfrac{an+b}c\rfloor-\dfrac12 g(\lfloor \dfrac{an+b}c\rfloor-1, c,c-b-1,a)\\-\dfrac12 f(\lfloor \dfrac{an+b}c\rfloor-1, c,c-b-1,a) &otherwise
\end{cases}
\]
发现 \(g,h\) 是互推的，\(f\) 是独立的，即可写出程序。
然而常数过大，无法接受，但是，我们欣喜地发现，在同一种 \(a,b,c\) 的关系之下，\(f,g,h\) 所需的递归都是一致的，因此我们不妨同步算 \(f,g,h\)，一起计算的时间复杂度为 \(O(\log n)\)
附上代码供参考：
#include <iostream>
#include <cstdio>
using namespace std;
typedef long long ll;
const int mod = 998244353;
ll inv2 = 499122177, inv6 = 166374059;
inline ll mul(ll a, ll b) { return (a * b) % mod; }
inline ll add(ll a, ll b) { return (a + b) % mod; }
inline ll sqr(ll a) { return mul(a, a); }
/*
ll f(ll n, ll a, ll b, ll c); // sum (ai + b) / c
ll g(ll n, ll a, ll b, ll c); // sum ((ai + b) / c) ^ 2
ll h(ll n, ll a, ll b, ll c); // sum i((ai + b) / c)
ll f(ll n, ll a, ll b, ll c) {
if (!a) {
return mul(n + 1, b / c);
} else if (a >= c || b >= c) {
ll co1 = f(n, a % c, b % c, c), co2 = mul(mul(mul(a / c, n), n + 1), inv2);
ll co3 = mul(b / c, n + 1);
return add(co1, add(co2, co3));
} else {
ll m = (a * n + b) / c;
return add(mul(m, n), mod - f(m - 1, c, c - b - 1, a));
}
}
ll g(ll n, ll a, ll b, ll c) {
if (!a) {
return mul(n + 1, sqr(b / c));
} else if (a >= c || b >= c) {
ll co1 = g(n, a % c, b % c, c), co2 = mul(sqr(a / c), mul(n, mul(n + 1, mul(2 * n + 1, inv6))));
ll co3 = mul(n + 1, sqr(b / c)), co4 = mul(2, mul(a / c, h(n, a % c, b % c, c)));
ll co5 = mul(a / c, mul(b / c, mul(n, n + 1))), co6 = mul(2, mul(b / c, f(n, a % c, b % c, c)));
return add(co1, add(co2, add(co3, add(co4, add(co5, co6)))));
} else {
ll m = (a * n + b) / c;
ll co1 = mul(n, mul(m, m + 1)), co2 = mod - mul(2, h(m - 1, c, c - b - 1, a));
ll co3 = mod - mul(2, f(m - 1, c, c - b - 1, a)), co4 = mod - f(n, a, b, c);
return add(co1, add(co2, add(co3, co4)));
}
}
ll h(ll n, ll a, ll b, ll c) {
if (!a) {
return mul(b / c, mul(n, mul(n + 1, inv2)));
} else if (a >= c || b >= c) {
ll co1 = h(n, a % c, b % c, c), co2 = mul(a / c, mul(n, mul(n + 1, mul(2 * n + 1, inv6))));
ll co3 = mul(b / c, mul(n, mul(n + 1, inv2)));
return add(co1, add(co2, co3));
} else {
ll m = (a * n + b) / c;
ll co1 = mul(m, mul(n, mul(n + 1, inv2))), co2 = mod - mul(g(m - 1, c, c - b - 1, a), inv2);
ll co3 = mod - mul(f(m - 1, c, c - b - 1, a), inv2);
return add(co1, add(co2, co3));
}
}
*/
struct node { ll f, g, h; };
node solve(ll n, ll a, ll b, ll c) {
node ret; ll co1, co2, co3, co4, co5, co6;
if (!a) {
ret.f = mul(n + 1, b / c);
ret.g = mul(n + 1, sqr(b / c));
ret.h = mul(b / c, mul(n, mul(n + 1, inv2)));
} else if (a >= c || b >= c) {
node tmp = solve(n, a % c, b % c, c);
co1 = tmp.f, co2 = mul(mul(mul(a / c, n), n + 1), inv2), co3 = mul(b / c, n + 1);
ret.f = add(co1, add(co2, co3));
co1 = tmp.g, co2 = mul(sqr(a / c), mul(n, mul(n + 1, mul(2 * n + 1, inv6))));
co3 = mul(n + 1, sqr(b / c)), co4 = mul(2, mul(a / c, tmp.h));
co5 = mul(a / c, mul(b / c, mul(n, n + 1))), co6 = mul(2, mul(b / c, tmp.f));
ret.g = add(co1, add(co2, add(co3, add(co4, add(co5, co6)))));
co1 = tmp.h, co2 = mul(a / c, mul(n, mul(n + 1, mul(2 * n + 1, inv6))));
co3 = mul(b / c, mul(n, mul(n + 1, inv2)));
ret.h = add(co1, add(co2, co3));
} else {
ll m = (a * n + b) / c;
node tmp = solve(m - 1, c, c - b - 1, a);
co1 = mul(m, n), co2 = mod - tmp.f;
ret.f = add(co1, co2);
co1 = mul(n, mul(m, m + 1)), co2 = mod - mul(2, tmp.h);
co3 = mod - mul(2, tmp.f), co4 = mod - ret.f;
ret.g = add(co1, add(co2, add(co3, co4)));
co1 = mul(m, mul(n, mul(n + 1, inv2))), co2 = mod - mul(tmp.g, inv2);
co3 = mod - mul(tmp.f, inv2);
ret.h = add(co1, add(co2, co3));
}
return ret;
}
int main() {
int T; scanf("%d", &T);
while (T--) {
ll n, a, b, c; scanf("%lld%lld%lld%lld", &n, &a, &b, &c);
node ans = solve(n, a, b, c);
printf("%lld %lld %lld\n", ans.f, ans.g, ans.h);
}
return 0;
}
