直接用例子开始吧：
\[给定 a,b,c,n，求解\sum_{i = 0}^n\left\lfloor\frac{ai+b}{c}\right\rfloor
\]
我们记 \(f(a,b,c,n) = \sum_{i = 0}^n\left\lfloor\frac{ai+b}{c}\right\rfloor\)
首先考虑 \(a \ge c\) 或 \(b \ge c\) 的情况：
先记 \(q_a = \left\lfloor\frac{a}{c}\right\rfloor,r_a = a \bmod c,q_b = \left\lfloor\frac{b}{c}\right\rfloor,r_b = b \bmod c\)
\[\begin{aligned}
\sum_{i = 0}^n\left\lfloor\frac{ai+b}{c}\right\rfloor &= \sum_{i = 0}^n\left\lfloor\frac{(q_a\cdot c+r_a)i+(q_b\cdot c+r_b)}{c}\right\rfloor\\
&= \sum_{i = 0}^n\left\lfloor\frac{r_ai+r_b}{c}+q_ai+q_b\right\rfloor\\
&= \sum_{i = 0}^n\left\lfloor\frac{r_ai+r_b}{c}\right\rfloor+q_ai+q_b
\end{aligned}
\]
于是有：
\[f(a,b,c,n) = \frac{n(n+1)}{2}q_a+(n+1)q_b+f(r_a,r_b,c,n)
\]
接下来考虑 \(a < c, b < c\) 时
记
\[m = \left\lfloor\frac{an+b}{c}\right\rfloor
\]
\[\left\lfloor\frac{ai+b}{c}\right\rfloor = \sum_{j = 0}^{m-1}\left[j < \left\lfloor\frac{ai+b}{c}\right\rfloor\right]
\]
放回原式中：
\[\sum_{i = 0}^n\sum_{j = 0}^{m-1}\left[j < \left\lfloor\frac{ai+b}{c}\right\rfloor\right] = \sum_{j = 0}^{m-1}\sum_{i = 0}^n\left[j < \left\lfloor\frac{ai+b}{c}\right\rfloor\right]
\]
考虑对 \(j\) 的条件变为对 \(i\) 的条件：
\[\begin{aligned}
\left[j < \left\lfloor\frac{ai+b}{c}\right\rfloor\right] &= \left[j+1 \le \left\lfloor\frac{ai+b}{c}\right\rfloor\right]\\
&= \left[j+1 \le \frac{ai+b}{c}\right]\\
&= \left[\frac{cj+c-b}{a} \le i\right]\\
&= \left[\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor < i\right]
\end{aligned}
\]
带回前式：
\[\begin{aligned}
f(a,b,c,n) &= \sum_{j = 0}^{m-1}\sum_{i = 0}^n \left[\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor < i\right]\\
&= \sum_{j = 0}^{m-1} n-\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor\\
&= nm - f(c,c-b-1,a,m-1)
\end{aligned}
\]
于是我们得到了一个递归的方法，由于 \(a\) 和 \(c\) 互换了位置，之后我们还会去用 \(a\) 模 \(c\)，复杂度和欧几里得方法是一样的。
然后是类欧的洛谷板子题
这道题里面我们还需要求：
\[\sum_{i = 0}^n \left\lfloor\frac{ai+b}{c}\right\rfloor^2,\sum_{i = 0}^n i\left\lfloor\frac{ai+b}{c}\right\rfloor
\]
凑出可行的形式就好了，由于这两个式子在递归的时候还需要求对方以及之前的式子，所以我们把三个要求的东西一块递归
代码见下：
#include <bits/stdc++.h>
using namespace std;
using llt = long long;
using ull = unsigned long long;
using uit = unsigned int;
using cit = const int;
using cllt = const long long;
int quick_pow(int _x,int _y,cit& _mod = 998244353) {int _res = 1;while(_y) {if(_y&1) _res = (ull)_res*_x%_mod;_x = (ull)_x*_x%_mod;_y >>= 1;}return _res;
}
int inv(cit& _x,cit& _mod = 998244353) {assert(_x);return quick_pow(_x,_mod-2,_mod);
}
template<cit _mod>
class mod_int {
private :constexpr static int mod = _mod;int _val;
public ://构造函数mod_int () : _val(0) {}mod_int (cit& _x) : _val(norm(_x)) {} mod_int (cllt& _x) : _val(norm(_x)) {}inline int &val() { return _val;}inline int norm(const llt &_x) { return (_x%mod+mod)%mod;}mod_int inv() const {assert(_val);mod_int res;res._val = ::inv(_val,mod);return res;}void set(cit &_x) { _val = _x;}#define cmi const mod_intmod_int (cmi& _x) : _val(_x._val) {}mod_int operator - () const { mod_int res(*this); res._val = mod-res._val; return res;}mod_int& operator = (cmi& _y) { _val = _y._val; return *this;}mod_int& operator += (cmi& _y) { _val += _y._val; if(_val >= mod) _val -= mod; return *this;}mod_int& operator -= (cmi& _y) { _val -= _y._val; if(_val < 0) _val += mod; return *this;}mod_int& operator *= (cmi& _y) { _val = (ull)_val*_y._val%mod; return *this;}mod_int& operator /= (cmi& _y) { _val = (ull)_val*_y.inv()%mod; return *this;}friend mod_int operator + (cmi& _a,cmi& _b) { mod_int res(_a); if((res._val += _b._val) >= mod) res._val -= mod; return res;}friend mod_int operator - (cmi& _a,cmi& _b) { mod_int res(_a); if((res._val -= _b._val) < 0) res._val += mod; return res;}friend mod_int operator * (cmi& _a,cmi& _b) { mod_int res; res._val = (ull)_a._val*_b._val%mod; return res;}friend mod_int operator / (cmi& _a,cmi& _b) { mod_int res; res._val = (ull)_a._val*::inv(_b._val,mod)%mod; return res;}friend istream &operator >> (istream &is,mod_int& a) { is >> a._val; a._val = a.norm(a._val); return is;}friend ostream &operator << (ostream &os,cmi& a) { return os << a._val;}friend bool operator > (cmi &_a,cmi &_b) { return _a._val > _b._val;}friend bool operator < (cmi &_a,cmi &_b) { return _a._val < _b._val;}friend bool operator >= (cmi &_a,cmi &_b) { return _a._val >= _b._val;}friend bool operator <= (cmi &_a,cmi &_b) { return _a._val <= _b._val;}operator bool () { return _val != 0;}friend mod_int operator + (cmi &_a,cllt &_b) { return mod_int(_a._val+_b);}friend mod_int operator + (cllt &_b,cmi &_a) { return mod_int(_a._val+_b);}friend mod_int operator - (cmi &_a,cllt &_b) { return mod_int(_a._val-_b);}friend mod_int operator - (cllt &_b,cmi &_a) { return mod_int(_b-_a._val);}friend mod_int operator * (cmi &_a,cllt &_b) { return mod_int(1ll*_a._val*(_b%mod));}friend mod_int operator * (cllt &_b,cmi &_a) { return mod_int(1ll*_a._val*(_b%mod));}friend mod_int operator / (cmi &_a,cllt &_b) { return mod_int(_a._val*::inv(_b,mod));}friend mod_int operator / (cllt &_b,cmi &_a) { return mod_int((_b%mod)*::inv(_a._val));}mod_int &operator = (cllt &_x) { _val = norm(_x); return *this;}#undef cmi
};
const int mod = 998244353;
const int mod_g = 3;
using mint = mod_int<mod>;
mint _1(1), _2(2), _6(6);
mint i2((mod+1)/2), i6((mod+1)/6);
mint sum_1(mint n) {return n*(n+_1)*i2;
}
mint sum_2(mint n) {return n*(n+_1)*(_2*n+_1)*i6;
}
struct NODE {mint f, g, h;NODE (const mint &_f = 0,const mint &_g = 0,const mint &_h = 0) : f(_f), g(_g), h(_h) {}
};
NODE solve(llt a,llt b,llt c,llt n) {mint n2(sum_1(n));mint n3(sum_2(n));NODE ans(0,0,0);if(a >= c||b >= c) {NODE res = solve(a%c,b%c,c,n);llt ta = a/c, tb = b/c;ans.f = res.f+n2*ta+(n+_1)*tb;ans.g = res.g+ta*n3+tb*n2;ans.h = res.h+tb*res.f*_2+ta*res.g*_2+n3*ta*ta+(n+_1)*tb*tb+n2*ta*tb*_2;return ans;}llt m = (1ll*a*n+b)/c;if(!m) return NODE();NODE res = solve(c,c-b-1,a,m-1);m %= mod;ans.f = n*m-res.f;ans.g = m*sum_1(n)-res.f*i2-res.h*i2;ans.h = n*m%mod*m-res.f-res.g*_2;return ans;
}
int T;
int a, b, c, n;
int main() {cin >> T;while(T--) {cin >> n >> a >> b >> c;NODE ans = solve(a,b,c,n);cout << ans.f << ' ' << ans.h << ' ' << ans.g << endl;}return 0;
}
感觉我的 mod_int 融入得很丑陋，哎