原理不难但是写起来非常复杂的东西。
我觉得讲得非常好懂的博客。 传送门
我们设
$$f(a, b, c, n) = sum_{i = 0}^{n}left lfloor frac{ai + b}{c} ight floor$$
$$g(a, b, c, n) = sum_{i = 0}^{n}ileft lfloor frac{ai + b}{c} ight floor$$
$$h(a, b, c, n) = sum_{i = 0}^{n}left lfloor frac{ai + b}{c} ight floor^2$$
先考虑一个结论
$$left lfloor frac{Ax}{y} ight floor = left lfloor frac{A(x mod y)}{y} ight floor + Aleft lfloor frac{x}{y} ight floor$$
那么有
$$left lfloor frac{ai + b}{c} ight floor = left lfloor frac{(a mod c)i + (b mod c)}{c} ight floor + ileft lfloor frac{a}{c} ight floor + left lfloor frac{b}{c} ight floor$$
这个东西可以把$a geq c$ 或者$b geq c$的情况转化成$a,b < c$的情况。
F
先看看这个比较好做的$f$。
注意到当$a geq c$或者$b geq c$的时候，我们可以用以上的结论把$a$或者$b$降下来，有
$$f(a, b, c, n) = left lfloor frac{b}{c} ight floor(n + 1) + left lfloor frac{a}{c} ight floor frac{n(n + 1)}{2} + f(a mod c, b mod c, c, n)$$
那么当$a, b < c$的时候，有
$$f(a, b, c, n) = sum_{i = 0}^{n}sum_{j = 1}^{m}[left lfloor frac{ai + b}{c} ight floor geq j]$$
为了方便把$left lfloor frac{an + b}{c} ight floor$记为$m$。
$$f(a, b, c, n) = sum_{i = 0}^{n}sum_{j = 0}^{m - 1}[left lfloor frac{ai + b}{c} ight floor geq j + 1]$$
$$= sum_{i = 0}^{n}sum_{j = 0}^{m - 1}[ai geq cj + c - b]$$
$$= sum_{i = 0}^{n}sum_{j = 0}^{m - 1}[ai > cj + c - b - 1]$$
$$= sum_{i = 0}^{n}sum_{j = 0}^{m - 1}[i > left lfloor frac{cj + c - b - 1}{a} ight floor]$$
注意到这时候$i$这一项可以直接算了。
$$f(a, b, c, n) = sum_{j = 0}^{m - 1}(n - left lfloor frac{cj + c - b - 1}{a} ight floor)$$
$$= nm - f(c, c - b - 1, a, m - 1)$$
比较方便。
G
按照套路先做$a geq c$或者$b geq c$的情况。
$$g(a, b, c, n) = left lfloor frac{b}{c} ight floorfrac{n(n + 1)}{2} + left lfloor frac{a}{c} ight floor frac{n(n + 1)(2n + 1)}{6} + g(a mod c, b mod c, c, n)$$
然后接着按照上述套路弄$a, b < c$的情况。
$$g(a, b, c, n) = sum_{j = 0}^{m - 1}frac{(n - left lfloor frac{cj + c - b - 1}{a} ight floor)(n + 1 + left lfloor frac{cj + c - b - 1}{a} ight floor)}{2}$$
$$= frac{1}{2}(mn(n + 1) - f(c, c - b - 1, a, m - 1) - h(c, c - b - 1, a, m - 1))$$
还要解决$h$。
H
当$a geq c$或者$b geq c$的时候
$$h(a, b, c, n) = (n + 1)left lfloor frac{b}{c} ight floor^2 + frac{n(n + 1)(2n + 1)}{6}left lfloor frac{a}{c} ight floor^2 + n(n + 1)left lfloor frac{b}{c} ight floorleft lfloor frac{a}{c} ight floor + h(a mod c, b mod c, c, n) + 2left lfloor frac{a}{c} ight floor g(a mod c, b mod c, c, n) + 2left lfloor frac{b}{c} ight floor f(a mod c, b mod c, c, n)$$
好麻烦啊……
当$a, b < c$的时候需要一些操作
$$n^2 = 2 imes frac{n(n + 1)}{2} - n = 2sum_{i = 1}^{n}i - n$$
可以得到
$$h(a, b, c, n) = sum_{i = 0}^{n}(2sum_{j = 1}^{left lfloor frac{ai + b}{c} ight floor}j - left lfloor frac{ai + b}{c} ight floor)$$
$$= 2sum_{i = 0}^{n}sum_{j = 1}^{left lfloor frac{ai + b}{c} ight floor}j - f(a, b, c, n)$$
$$= 2sum_{j = 0}^{m - 1}(j + 1)sum_{i = 0}^{n}[left lfloor frac{ai + b}{c} ight floor geq j + 1] - f(a, b, c, n)$$
$$= 2sum_{j = 0}^{m - 1}(j + 1)(n - left lfloor frac{cj + c - b - 1}{a} ight floor) - f(a, b, c, n)$$
$$= nm(m + 1) - 2g(c, c - b - 1, a, m - 1) - 2f(c, c - b - 1, a, m - 1) - f(a, b, c, n)$$
大概就是这样了。
在计算的时候如果只考虑$(a, c)$这两项的话相当于每一次把$(a, c)$变成了$(c, a \% c)$，所以时间复杂度和欧几里得算法相同，递归层数是$log$层。
时间复杂度$O(Tlogn)$。
因为递归的式子很一致，在计算的时候应当三个东西一起算比较快。
Code：
#include <cstdio> #include <cstring> using namespace std; typedef long long ll; const ll P = 998244353LL; template <typename T> inline void read(T &X) { X = 0; char ch = 0; T op = 1; for (; ch > '9'|| ch < '0'; ch = getchar()) if (ch == '-') op = -1; for (; ch >= '0' && ch <= '9'; ch = getchar()) X = (X << 3) + (X << 1) + ch - 48; X *= op; } template <typename T> inline void inc(T &x, T y) { x += y; if (x >= P) x -= P; } template <typename T> inline void sub(T &x, T y) { x -= y; if (x < 0) x += P; } inline ll fpow(ll x, ll y) { ll res = 1; for (; y > 0; y >>= 1) { if (y & 1) res = res * x % P; x = x * x % P; } return res; } namespace Likegcd { struct Node { ll f, g, h; }; #define f(now) now.f #define g(now) now.g #define h(now) now.h const ll inv2 = fpow(2, P - 2); const ll inv6 = fpow(6, P - 2); inline Node solve(ll a, ll b, ll c, ll n) { Node res; if (!a) { f(res) = (b / c) * (n + 1) % P; g(res) = (b / c) * (n + 1) % P * n % P * inv2 % P; h(res) = (b / c) * (b / c) % P * (n + 1) % P; return res; } f(res) = g(res) = h(res) = 0; if (a >= c || b >= c) { Node tmp = solve(a % c, b % c, c, n); inc(f(res), (a / c) * n % P * (n + 1) % P * inv2 % P); inc(f(res), (b / c) * (n + 1) % P); inc(f(res), f(tmp)); inc(g(res), (a / c) * n % P * (n + 1) % P * ((2 * n + 1) % P) % P * inv6 % P); inc(g(res), (b / c) * n % P * (n + 1) % P * inv2 % P); inc(g(res), g(tmp)); inc(h(res), (a / c) * (a / c) % P * n % P * (n + 1) % P * ((2 * n + 1) % P) % P * inv6 % P); inc(h(res), (b / c) * (b / c) % P * (n + 1) % P); inc(h(res), (a / c) * (b / c) % P * n % P * (n + 1) % P); inc(h(res), h(tmp)); inc(h(res), 2LL * (a / c) % P * g(tmp) % P); inc(h(res), 2LL * (b / c) % P * f(tmp) % P); return res; } if (a < c && b < c) { ll m = (a * n + b) / c; Node tmp = solve(c, c - b - 1, a, m - 1); f(res) = n * m % P; sub(f(res), f(tmp)); g(res) = n * (n + 1) % P * m % P; sub(g(res), f(tmp)); sub(g(res), h(tmp)); g(res) = g(res) * inv2 % P; h(res) = n * m % P * (m + 1) % P; sub(h(res), 2LL * g(tmp) % P); sub(h(res), 2LL * f(tmp) % P); sub(h(res), f(res)); return res; } return res; } } int main() { int testCase; read(testCase); for (ll a, b, c, n; testCase--; ) { read(n), read(a), read(b), read(c); Likegcd :: Node ans = Likegcd :: solve(a, b, c, n); printf("%lld %lld %lld ", ans.f, ans.h, ans.g); } return 0; }