类欧几里德算法
经典模型
\[f(a,b,c,n)=\sum_{i=0}^n \lfloor \frac{ai+b}{c} \rfloor
\]
\(a,b,c,n \le 10^9,c \not= 0\)，要求在 \(\log n\) 的时间内求出。
约定
为了方便表示，做以下规定：
\(x/y\) 表示 \(\lfloor \frac{a}{b} \rfloor\)
\(a' = a \bmod c\)
情况一
\[a,b \ge c
\]
这个时候我们可以将 \(a,b\) 中大于 \(c\) 的部分拿出来，以简化问题：
\[\begin{aligned}
f(a,b,c,n ) &= \sum_{i=0}^n \lfloor \frac{ai+b}{c} \rfloor\\
&= \frac{n \times (n + 1) \times a/c}{2} + (n+1) * (b/c)+ \sum_{i = 0}^n \lfloor \frac{a'i+b'}{c} \rfloor\\
&= \frac{n \times (n + 1) \times a/c}{2} + (n+1) * (b/c)+ f(a',b',c,n)
\end{aligned}
\]
情况二
\[a,b < c
\]
这时候我们需要一些 trick 来想方设法地交换求和号：
\[\begin{aligned}
f(a,b,c,n) &= \sum_{i=0}^n \lfloor \frac{ai+b}{c} \rfloor\\
&= \sum_{i=0}^n \sum_{j=0}^{(ai+b)/c-1}1\\
&= \sum_{j=0}^{(an+b)/c-1}\sum_{i=0}^n[j<(ai+b)/c]
\end{aligned}
\]
那个 \(j < (ai+b)/c\) 不是很好搞，但是如果关于底和顶的基本功比较好的话应该是问题不大的。
先抛俩较为基础的结论:（\(n\) 为整数，\(x\) 为实数）
\[n \le x \Leftrightarrow n \le \lfloor x \rfloor
\]
\[n > x \Leftrightarrow n > \lfloor x \rfloor
\]
然后推：
\[j < (ai+b)/c\\
j+1 \le (ai+b)/c\\
j+1 \le \frac{ai+b}{c}\\
cj+c \le ai+b\\
cj+c-b \le ai\\
cj+c-b-1 < ai\\
i > (cj+c-b-1)/a
\]
然后原式子就好看很多了。定义 \(m = (an+b)/c\)，然后：
\[\begin{aligned}
f(a,b,c,n) &= \sum_{j=0}^{(an+b)/c-1}\sum_{i=0}^n[j<(ai+b)/c]\\
&= \sum_{j=0}^{m-1}n-\lfloor \frac{cj+c-b-1}{a} \rfloor\\
&= nm-f(c,c-b-1,a,n)
\end{aligned}
\]
感觉没啥用？其实是有用的，观察到 \(a,c\) 换了位置，而情况一可以将 \(a\) 模 \(c\)，因此复杂度和辗转相除法的复杂度一样，都是 \(\log n\)。
情况三(边界)
\[a=0
\]
\[f(a,b,c,n)=(n+1)\lfloor \frac{b}{c} \rfloor
\]
拓展
现在又跑出来俩类似 \(f\) 的函数 \(g\) 和 \(h\)，\(h\) 是二次的 \(f\)，\(g\) 是 \(h\) 的辅助数组：
\[g(a,b,c,n) = \sum_{i = 0}^n i \lfloor \frac{ai+b}{c} \rfloor\\
h(a,b,c,n) = \sum_{i=0}^n \lfloor \frac{ai+b}c \rfloor ^2
\]
这俩式子的推导与 \(f\) 的思想类似，都是分为三种情况，第一种取模，第二种先办法搞出求和号然后交换以实现 \(a,c\) 互换位置的效果，第三种边界。第三种太水了，就不写了。
g函数
情况一
\[\begin{aligned}
g(a,b,c,n) &= \sum_{i=0}^n i \lfloor \frac{ai+b}c \rfloor\\
&= \sum_{i=0}^n i(i \times a/c + b/c) \lfloor \frac{a'i+b'}c \rfloor \\
&= a/c \times n(n+1)(2n+1)/6 + b/c \times n(n+1)/2 + g(a',b',c,n)
\end{aligned}
\]
情况二
\[\begin{aligned}
g(a,b,c,n) &= \sum_{i=0}^n i \sum_{j=0}^{(ai+b)/c-1}1\\
&= \sum_{j=0}^{m-1}\sum_{i=0}^ni[j < (ai+b)/c]\\
&= \sum_{j=0}^{m-1}n(n+1)/2 - (jc+c-b-1)/a \times ((jc+c-b-1)/a-1)/2\\
&= \frac{1}{2}(mn(n+1)-h(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1))
\end{aligned}
\]
h函数
情况一
公式太难打了，就直接给结论吧。反正推导全是暴力拆式子
\[h(a,b,c,n) = \\
\small (a/c)^2n(n+1)(2n+1)/6 + (a/c)(b/c)n(n+1)+(n+1)(b/c)^2+2(a/c)g(a',b',c,n)+2(b/c)f(a',b',c,n)+h(a',b',c,n)
\]
情况二
这里还需要一个小 trick：（普通幂表示成自然数幂和）
\[n^2 = (2\sum_{i=0}^{n-1}i) + n
\]
然后就可以推了：
\[\begin{aligned}
h(a,b,c,n) &= \sum_{i=0}^n 2 \sum_{j=0}^{(ai+b)/c-1}j + f(a,b,c,n)\\
&= 2 \sum_{j=0}^{m-1}j \sum_{i=0}^n[(ai+ b)/c > j] + f(a,b,c,n)\\
&= 2\sum_{j=0}^{m-1}j(n-\lfloor \frac{cj+c-b-1}a \rfloor)+f(a,b,c,n)\\
&= nm(m-1) - 2g(c,c-b-1,a,m-1)+f(a,b,c,n)
\end{aligned}
\]
由于 \(f\)， \(g\) 和 \(h\) 需要相互调用，我们可以每次打包返回三个值。
题目与模板
类欧模板（调试用）
struct node {
ll f, g, h;
node() {}
node(ll a, ll b, ll c) { f = a, g = b, h = c; }
};
node sol(int a, int b, int c, int n) {
ll nwf = 0, nwg = 0, nwh = 0;
if (!a) {
nwf = (b / c) * (n + 1) % P;
nwg = (b / c) * n % P * (n + 1) % P * inv2 % P;
nwh = (b / c) * (b / c) % P * (n + 1) % P;
return node(nwf, nwg, nwh);
}
if (a >= c || b >= c) {
node nd = sol(a % c, b % c, c, n);
ll tf = nd.f, tg = nd.g, th = nd.h;
nwf = ((a / c) * n % P * (n + 1) % P * inv2 + (b / c) * (n + 1) + tf) % P;//bug
nwg = ((a / c) * n % P * (n + 1) % P * (2ll * n + 1) % P * inv6 + (b / c) * n % P * (n + 1) % P * inv2 + tg) % P;
nwh = (a / c) * (a / c) % P * n % P * (n + 1) % P * (2ll * n + 1) % P * inv6;
nwh = (nwh + (a / c) * (b / c) % P * n % P * (n + 1) + (n + 1) * (b / c) % P * (b / c)) % P;
nwh = (nwh + 2ll * (b / c) * tf + 2ll * (a / c) * tg + th) % P;
return node(nwf, nwg, nwh);
}
ll m = (a * n + b) / c;//bug
node nd = sol(c, c - b - 1, a, m - 1);
m %= P;
ll tf = nd.f, tg = nd.g, th = nd.h;
nwf = (m * n - tf) % P;
nwg = (m * n % P * (n + 1) - th - tf) % P * inv2 % P;
nwh = (n * m % P * (m - 1) - 2ll * tg + nwf) % P;
return node(nwf, nwg, nwh);
}
注意一下取模的问题，计算答案要取模，但是向下递归的 \(m\) 不能扔个 \(m \bmod P\) 进去。
