C++类欧几里得算法

（本文用 $绿色$标注一切函数，用 $紫色$标注部分可以 $O(1)$计算的常量式，便于观察）（不知道为什么加了颜色过后括号不能自动变大了，将就看吧= =）

概述

类欧几里得算法用于求形如$f(a,b,c,n)=\sum _{i=0}^n\lfloor \frac{ai+b}{c}\rfloor$的式子。这个式子有个几何意义，表示求下图中阴影部分中的整点个数（直线为$y=ax+b$）：

算法

令$a^{\prime }=a\%c,b^{\prime }=b\%c$，由定义式得$$\begin{array}{rl}f(a,b,c,n)& =\sum _{i=0}^n\left(\lfloor \frac{a^{\prime }i+b^{\prime }}{c}\rfloor +\lfloor \frac{a}{c}\rfloor i+\lfloor \frac{b}{c}\rfloor \right)\\ & =f(a^{\prime },b^{\prime },c,n)+\lfloor \frac{a}{c}\rfloor \sum _{i=0}^{n}i+(n+1)\lfloor \frac{b}{c}\rfloor \end{array}$$于是只需要处理$a,b<c$的情况。

设$m=\lfloor \frac{an+b}{c}\rfloor$，则$$\begin{array}{rl}f(a,b,c,n)& =\sum _{i=0}^n\sum _{j=0}^m\left[j\le \lfloor \frac{ai+b}{c}\rfloor \right]\\ & =\sum _{i=0}^n\sum _{j=0}^{m-1}\left[j<\lfloor \frac{ai+b}{c}\rfloor \right]\\ & =\sum _{j=0}^{m-1}\sum _{i=0}^n\left[j<\lfloor \frac{ai+b}{c}\rfloor \right]\end{array}$$

现在要想办法把里面变成$i$为主元。
如下图所示，可以看出$0\le j<5$，结合图示可以得到$$j<\lfloor \frac{ai+b}{c}\rfloor \iff (j+1)c\le ai+b$$
即$(j+1)c<ai+b+1$，即$$i>\lfloor \frac{cj-b+c-1}{a}\rfloor$$于是可以继续化简了$$\begin{array}{rl}f(a,b,c,n)& =\sum _{j=0}^{m-1}\sum _{i=0}^n\left[i>\lfloor \frac{cj-b+c-1}{a}\rfloor \right]\\ & =\sum _{j=0}^{m-1}\left(n-\lfloor \frac{cj-b+c-1}{a}\rfloor \right)\\ & =mn-\sum _{j=0}^{m-1}\lfloor \frac{cj-b+c-1}{a}\rfloor \\ & =mn-f(c,c-b-1,a,m-1)\end{array}$$
仔细康康这个递推式，发现是由$f(a,b,c,n)$变成了$f(c,c-b\%c-1,a\%c,m-1)$，$a,c$两个参数的变化跟欧几里得算法一模一样，所以它的复杂度跟欧几里得求最大公约数的算法一样，也因此称它为类欧。

版题

洛谷P5170 【模板】类欧几里得算法

它还让你求$g(a,b,c,n)=\sum _{i=0}^n{\lfloor \frac{ai+b}{c}\rfloor }^2$，以及$h(a,b,c,n)=\sum _{i=0}^n\left(i\lfloor \frac{ai+b}{c}\rfloor \right)$，直接刚，令$a^{\prime }=a\%c,b^{\prime }=b\%c$：$$\begin{array}{rl}g(a,b,c,n)=& \sum _{i=0}^n{\left(\lfloor \frac{a^{\prime }i+b^{\prime }}{c}\rfloor +\lfloor \frac{a}{c}\rfloor i+\lfloor \frac{b}{c}\rfloor \right)}^2\\ =& g(a^{\prime },b^{\prime },c,n)+{\lfloor \frac{a}{c}\rfloor }^2\sum _{i=0}^n{i}^2+\sum _{i=0}^n{\lfloor \frac{b}{c}\rfloor }^2+2\sum _{i=0}^n\left(\lfloor \frac{a^{\prime }i+b^{\prime }}{c}\rfloor \lfloor \frac{a}{c}\rfloor i\right)\\ & +2\sum _{i=0}^n\left(\lfloor \frac{a^{\prime }i+b^{\prime }}{c}\rfloor \lfloor \frac{b}{c}\rfloor \right)+2\sum _{i=0}^n\left(i\lfloor \frac{a}{c}\rfloor \lfloor \frac{b}{c}\rfloor \right)\\ =& g(a^{\prime },b^{\prime },c,n)+{\lfloor \frac{a}{c}\rfloor }^2\sum _{i=0}^n{i}^2+\sum _{i=0}^n{\lfloor \frac{b}{c}\rfloor }^2+2\lfloor \frac{a}{c}\rfloor h(a^{\prime },b^{\prime },c,n)\\ & +2\lfloor \frac{b}{c}\rfloor f(a^{\prime },b^{\prime },c,n)+2\sum _{i=0}^n\left(i\lfloor \frac{a}{c}\rfloor \lfloor \frac{b}{c}\rfloor \right)\\ =& g(a^{\prime },b^{\prime },c,n)+2\lfloor \frac{a}{c}\rfloor h(a^{\prime },b^{\prime },c,n)+2\lfloor \frac{b}{c}\rfloor f(a^{\prime },b^{\prime },c,n)\\ & +\sum _{i=0}^n{\lfloor \frac{b}{c}\rfloor }^2+{\lfloor \frac{a}{c}\rfloor }^2\sum _{i=0}^n{i}^2+2\sum _{i=0}^n\left(i\lfloor \frac{a}{c}\rfloor \lfloor \frac{b}{c}\rfloor \right)\\ h(a,b,c,n)=& \sum _{i=0}^n\left(i\lfloor \frac{a^{\prime }i+b^{\prime }}{c}\rfloor +\lfloor \frac{a}{c}\rfloor i^2+\lfloor \frac{b}{c}\rfloor i\right)\\ =& h(a^{\prime },b^{\prime },c,n)+\lfloor \frac{a}{c}\rfloor \sum _{i=0}^n{i}^2+\lfloor \frac{b}{c}\rfloor \sum _{i=0}^{n}i\end{array}$$接下来考虑$a,b<c$的情况，这里有一个显然而神奇的变形：$$x^2=2\cdot \frac{x(x+1)}{2}-x=2\sum _{i=0}^{n}i-x$$于是，令$m=\lfloor \frac{an+b}{c}\rfloor$：$$\begin{array}{rl}g(a,b,c,n)=& \sum _{i=0}^n{\lfloor \frac{ai+b}{c}\rfloor }^2\\ =& \sum _{i=0}^n\left(2\sum _{j=0}^{\lfloor \frac{ai+b}{c}\rfloor }j-\lfloor \frac{ai+b}{c}\rfloor \right)\\ =& 2\sum _{i=0}^n\sum _{j=0}^{\lfloor \frac{ai+b}{c}\rfloor }j-f(a,b,c,n)\\ =& 2\sum _{j=0}^{m-1}\left((j+1)\sum _{i=0}^n\left[j<\lfloor \frac{ai+b}{c}\rfloor \right]\right)-f(a,b,c,n)\\ =& 2\sum _{j=0}^{m-1}\left((j+1)\sum _{i=0}^n\left[i>\lfloor \frac{cj-b+c-1}{a}\rfloor \right]\right)-f(a,b,c,n)\\ =& 2\sum _{j=0}^{m-1}\left((j+1)(n-\lfloor \frac{cj-b+c-1}{a}\rfloor )\right)-f(a,b,c,n)\\ =& 2n\sum _{j=0}^{m-1}(j+1)-2\sum _{j=0}^{m-1}\left((j+1)\lfloor \frac{cj-b+c-1}{a}\rfloor \right)-f(a,b,c,n)\\ =& 2n\sum _{j=0}^{m-1}(j+1)-2(h(c,c-b-1,a,m-1)+f(c,c-b-1,a,m-1))-f(a,b,c,n)\\ h(a,b,c,n)=& \sum _{i=0}^n\left(i\lfloor \frac{ai+b}{c}\rfloor \right)\\ =& \sum _{j=0}^{m-1}\sum _{i=0}^n\left(i\left[j<\lfloor \frac{ai+b}{c}\rfloor \right]\right)\\ =& \sum _{j=0}^{m-1}\sum _{i=0}^n\left(i\left[i<\lfloor \frac{cj-b+c-1}{a}\rfloor \right]\right)\\ =& \sum _{j=0}^{m-1}\sum _{i=\lfloor \frac{cj-b+c-1}{a}\rfloor +1}^{n}i\\ =& \sum _{j=0}^{m-1}\left(\frac{1}{2}\left(n+\lfloor \frac{cj-b+c-1}{a}\rfloor +1\right)\left(n-\lfloor \frac{cj-b+c-1}{a}\rfloor \right)\right)\\ =& \frac{1}{2}\sum _{j=0}^{m-1}(n(n+1))+\frac{1}{2}n\sum _{j=0}^{m-1}\lfloor \frac{cj-b+c-1}{a}\rfloor -\frac{1}{2}\sum _{j=0}^{m-1}{\lfloor \frac{cj-b+c-1}{a}\rfloor }^2\\ & -\frac{1}{2}(n+1)\sum _{j=0}^{m-1}\lfloor \frac{cj-b+c-1}{a}\rfloor \\ =& \frac{1}{2}\sum _{j=0}^{m-1}(n(n+1))+\frac{1}{2}g(c,c-b-1,a,m-1)-\frac{1}{2}f(c,c-b-1,a,m-1)\end{array}$$

代码
用结构体存，同步计算。
每个字符都要看清楚！

#include <cstdio>

typedef long long LL;

const LL MOD = 998244353;
const LL Inv2 = 499122177, Inv6 = 166374059;

struct Data {
	LL f, g, h;
	Data() { f = g = h = 0; }
	Data(LL _f, LL _g, LL _h) {
		f = _f, g = _g, h = _h;
	}
};

Data Euclid(LL a, LL b, LL c, LL n) {
	LL a1 = a / c, b1 = b / c, m = (a * n + b) / c;
	if (a == 0)
		return Data(b1 * (n + 1) % MOD,
					b1 * b1 % MOD * (n + 1) % MOD,
					b1 * n % MOD * (n + 1) % MOD * Inv2 % MOD);
	if (a >= c || b >= c) {
		Data T = Euclid(a % c, b % c, c, n);
		return Data((T.f + a1 * n % MOD * (n + 1) % MOD * Inv2 % MOD + (n + 1) * b1 % MOD) % MOD,
					(T.g + 2 * a1 % MOD * T.h % MOD + 2 * b1 % MOD * T.f % MOD + (n + 1) * b1 % MOD * b1 % MOD + a1 * b1 % MOD * n % MOD * (n + 1) % MOD + a1 * a1 % MOD * n % MOD * (n + 1) % MOD * (2 * n + 1) % MOD * Inv6 % MOD) % MOD,
					(T.h + a1 * n % MOD * (n + 1) % MOD * (2 * n + 1) % MOD * Inv6 % MOD + b1 * n % MOD * (n + 1) % MOD * Inv2 % MOD) % MOD);
	}
	Data T = Euclid(c, c - b - 1, a, m - 1);
	Data R(((m * n % MOD - T.f) % MOD + MOD) % MOD, 0, 0);
	return Data(R.f,
				((n * m % MOD * (m + 1) % MOD - 2 * T.h % MOD - 2 * T.f % MOD - R.f) % MOD + MOD) % MOD,
				((m * n % MOD * (n + 1) % MOD - T.g - T.f) % MOD + MOD) % MOD * Inv2 % MOD);
}

int main() {
	int T; scanf("%d", &T);
	while (T--) {
		int N, A, B, C;
		scanf("%d%d%d%d", &N, &A, &B, &C);
		Data Ans = Euclid(A, B, C, N);
		printf("%lld %lld %lld\n", Ans.f, Ans.g, Ans.h);
	}
}

---

强推OI WiKi的类欧几里得算法！