令人心态爆炸的类欧几里得算法（大波公式预警）

前言

你以为类欧几里得算法的难度\(≈\)扩展欧几里得\(≈\)欧几里得？

实际上，除了同样是用递归实现，复杂度证明方式差不多，二者根本没有任何共同之处。

被生活所迫来学习。。。如此多的公式，真的是令人心态爆炸啊。。。

更难受的是，不知道为什么本地Typora炸了，只能在线编辑，这么多公式打得让我想死。。。

一道板子题

我们来看看这道板子题吧：【洛谷5170】【模板】类欧几里得算法。

题目大意就是让我们求\(\sum_{i=0}^n\lfloor\frac{ai+b}c\rfloor,\sum_{i=0}^n\lfloor\frac{ai+b}c\rfloor^2,\sum_{i=0}^ni\lfloor\frac{ai+b}c\rfloor\)，其中\(n,a,b,c\le10^6\)。

考虑我们设：

\[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}^ni\lfloor\frac{ai+b}c\rfloor \]

接下来就是激动人心的推式子了。

\(f(n,a,b,c)\)

让我们从最简单的一个开始吧。

对于一个函数，我们要分三种情况讨论：

当\(a=0\)时：

\[f(n,a,b,c)=\sum_{i=0}^n\lfloor\frac bc\rfloor=(n+1)\times\lfloor\frac bc\rfloor \]

当\(a\ge c\)或\(b\ge c\)时：

首先你需要知道一个很显然的性质：

\[\lfloor\frac {ai+b}c\rfloor=\lfloor\frac{(a\%c)i+(b\%c)}c\rfloor+\lfloor\frac ac\rfloor i+\lfloor\frac bc\rfloor \]

然后我们就可以对原式进行转化：

\[\begin{align}f(n,a,b,c)&=\sum_{i=0}^n\lfloor\frac{(a\%c)i+(b\%c)}c\rfloor+\lfloor\frac ac\rfloor i+\lfloor\frac bc\rfloor\\&=f(n,a\%c,b\%c,c)+\lfloor\frac ac\rfloor\times\frac{n(n+1)}2+\lfloor\frac bc\rfloor\times(n+1)\end{align} \]

容易发现，这样就可以把该情况转化为\(a<c,b<c\)的情况。

当\(a<c,b<c\)时：

首先要注意对于下取整的一个转化，这是类欧几里得算法的通用套路。

还有，下面式子中\(m=\lfloor\frac{an+b}c\rfloor\)，显然它是\(\lfloor\frac{ai+b}c\rfloor\)能取到的最大值。由于它在接下来的式子中会经常出现（基本上每个式子都有它），所以这里用一个\(m\)替代。

剩余的就不多做解释了，相信大家都是能理解的。

\[\begin{align}f(n,a,b,c) & =\sum_{i=0}^n\sum_{j=1}^m[\lfloor\frac{ai+b}c\rfloor\ge j]\\ & =\sum_{i=0}^n\sum_{j=0}^{m-1}[\lfloor\frac{ai+b}c\rfloor\ge j+1]\\ & =\sum_{i=0}^n\sum_{j=0}^{m-1}[ai+b\ge cj+c]\\ & =\sum_{i=0}^n\sum_{j=0}^{m-1}[ai+b+1>cj+c]\\ & =\sum_{i=0}^n\sum_{j=0}^{m-1}[i>\lfloor\frac{cj+c-b-1}a\rfloor]\\ & =\sum_{j=0}^{m-1}(n-\lfloor\frac{cj+c-b-1}a\rfloor)\\ & =n\times m-f(m-1,c,c-b-1,a)\end{align} \]

这样一来我们就得到了\(f(n,a,b,c)\)的递归求解式，是不是很轻松？

但接下来，最可怕的就要来了。

\(g(n,a,b,c)\)

个人认为\(g(n,a,b,c)\)是这三个当中最难、最烦的一个。

要特别注意的是，在它的推导中会用到\(h\)，而\(h\)的推导也同样会用到\(g\)，二者是相互关联的。

接下来我们同样分三种情况讨论。

当\(a=0\)时：

\[g(n,a,b,c)=\sum_{i=0}^n\lfloor\frac bc\rfloor^2=(n+1)\times\lfloor\frac bc\rfloor^2 \]

当\(a\ge c\)或\(b\ge c\)时：

\[\begin{align}g(n,a,b,c)&=\sum_{i=0}^n(\lfloor\frac{(a\%c)i+(b\%c)}c\rfloor+\lfloor\frac ac\rfloor i+\lfloor\frac bc\rfloor)^2\\&=\sum_{i=0}^n(\lfloor\frac{(a\%c)i+(b\%c)}c\rfloor^2+\lfloor\frac ac\rfloor^2i^2+\lfloor\frac bc\rfloor^2+2i\lfloor\frac ac\rfloor \lfloor\frac{(a\%c)i+(b\%c)}c\rfloor+2\lfloor\frac bc\rfloor\lfloor\frac{(a\%c)i+(b\%c)}c\rfloor+2i\lfloor\frac ac\rfloor\lfloor\frac bc\rfloor)\\&=g(n,a\%c,b\%c,c)+\frac{n(n+1)(2n+1)}6\lfloor\frac ac\rfloor^2+(n+1)\lfloor\frac bc\rfloor^2+2\lfloor\frac ac\rfloor h(n,a\%c,b\%c,c)+2\lfloor\frac bc\rfloor f(n,a\%c,b\%c,c)+n(n+1)\lfloor\frac ac\rfloor\lfloor\frac bc\rfloor\end{align} \]

这一大坨式子真是写得要心态爆炸。。。

当\(a<c,b<c\)时：

我们发现，此时如果直接做，由于是平方项，值域是\(n^2\)的，因此我们需要一个转化：

\[x^2=2\sum_{i=1}^xi-x \]

证明？直接高斯求和就好了吧。。。

那么原式可以转化为这样：

\[\begin{align}g(n,a,b,c)&=\sum_{i=0}^n(2\sum_{j=1}^{\lfloor\frac{ai+b}c\rfloor}j-\lfloor\frac{ai+b}c\rfloor)\\&=2\sum_{i=0}^n\sum_{j=1}^{\lfloor\frac{ai+b}c\rfloor}j-f(n,a,b,c)\end{align} \]

单独提出前半部分，继续推式子：（你会发现和上面很像）

\[\begin{align}2\sum_{i=0}^n\sum_{j=1}^{\lfloor\frac{ai+b}c\rfloor}j&=2\sum_{i=0}^n\sum_{j=1}^mj[j\le\lfloor\frac{ai+b}c\rfloor]\\&=2\sum_{i=0}^n\sum_{j=0}^{m-1}(j+1)[i>\lfloor\frac{cj+c-b-1}a\rfloor]\\&=2\sum_{j=0}^{m-1}(j+1)(n-\lfloor\frac{cj+c-b-1}a\rfloor)\\&=2\sum_{j=0}^{m-1}(n(j+1)-j\lfloor\frac{cj+c-b-1}a\rfloor-\lfloor\frac{cj+c-b-1}a\rfloor)\\&=nm(m+1)-2h(m-1,c,c-b-1,a)-2f(m-1,c,c-b-1,a)\end{align} \]

代入回原式：

\[g(n,a,b,c)=nm(m+1)-2h(m-1,c,c-b-1,a)-2f(m-1,c,c-b-1,a)-f(n,a,b,c) \]

\(h(n,a,b,c)\)

当\(a=0\)时：

\[h(n,a,b,c)=\sum_{i=0}^ni\lfloor\frac bc\rfloor=\frac{n(n+1)}2\times\lfloor\frac bc\rfloor \]

当\(a\ge c\)或\(b\ge c\)时：

\[\begin{align}h(n,a,b,c)&=\sum_{i=0}^ni(\lfloor\frac{(a\%c)i+(b\%c)}c\rfloor+\lfloor\frac ac\rfloor i+\lfloor\frac bc\rfloor)\\&=\sum_{i=0}^n(i\lfloor\frac{(a\%c)i+(b\%c)}c\rfloor+\lfloor\frac ac\rfloor i^2+\lfloor\frac bc\rfloor i)\\&=h(n,a\%c,b\%c,c)+\frac{n(n+1)(2n+1)}6\lfloor\frac ac\rfloor+\frac{n(n+1)}2\lfloor\frac bc\rfloor\end{align} \]

当\(a<c,b<c\)时：

\[\begin{align}h(n,a,b,c)&=\sum_{i=0}^ni\sum_{j=1}^m[\lfloor\frac{ai+b}c\rfloor\ge j]\\&=\sum_{i=0}^ni\sum_{j=0}^{m-1}[i>\lfloor\frac{cj+c-b-1}a\rfloor]\\&=\sum_{j=0}^{m-1}\sum_{i=0}^ni[i>\lfloor\frac{cj+c-b-1}a\rfloor]\\&=\sum_{j=0}^{m-1}\frac{(\lfloor\frac{cj+c-b-1}a\rfloor+1+n)(n-\lfloor\frac{cj+c-b-1}a\rfloor)}2\\&=\frac12\sum_{j=0}^{m-1}(n(n+1)-\lfloor\frac{cj+c-b-1}a\rfloor-\lfloor\frac{cj+c-b-1}a\rfloor^2)\\&=\frac12(nm(n+1)-f(m-1,c,c-b-1,a)-g(m-1,c,c-b-1,a)\end{align} \]

具体实现

和求\(gcd\)过程差不多，都是递归实现的，相信大家都会。

尽管此题代码非常不可读（几乎全是1LL和取模），但还是贴一份吧。

代码

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define X 998244353
#define I2 499122177
#define I6 166374059
using namespace std;
int n,a,b,c;struct S {int f,g,h;};
class FastIO
{
	private:
		#define FS 100000
		#define tc() (A==B&&(B=(A=FI)+fread(FI,1,FS,stdin),A==B)?EOF:*A++)
		#define pc(c) (C==E&&(clear(),0),*C++=c)
		#define D isdigit(c=tc())
		int T;char c,*A,*B,*C,*E,FI[FS],FO[FS],S[FS];
	public:
		I FastIO() {A=B=FI,C=FO,E=FO+FS;}
		Tp I void read(Ty& x) {x=0;W(!D);W(x=(x<<3)+(x<<1)+(c&15),D);}
		Ts I void read(Ty& x,Ar&... y) {read(x),read(y...);}
		Tp I void write(Ty x) {W(S[++T]=x%10+48,x/=10);W(T) pc(S[T--]);}
		Tp I void write(Ty x,Ty y,Ty z) {write(x),pc(' '),write(y),pc(' '),write(z),pc('\n');}
		I void clear() {fwrite(FO,1,C-FO,stdout),C=FO;}
}F;
I S Get(CI n,CI a,CI b,CI c)//递归求解
{
	if(!a) return (S){1LL*(n+1)*(b/c)%X,1LL*(n+1)*(b/c)%X*(b/c)%X,1LL*n*(n+1)%X*I2%X*(b/c)%X};//a=0
	S t,res;if(a>=c||b>=c)//a>=c或b>=c
	{
		t=Get(n,a%c,b%c,c),res.f=(1LL*n*(n+1)%X*I2%X*(a/c)%X+1LL*(n+1)*(b/c)%X+t.f)%X,
		res.g=(1LL*n*(n+1)%X*(2*n+1)%X*I6%X*(a/c)%X*(a/c)%X+1LL*(n+1)*(b/c)%X
			*(b/c)%X+2LL*(a/c)*t.h%X+2LL*(b/c)*t.f%X+1LL*n*(n+1)%X*(a/c)%X*(b/c)%X+t.g)%X;
		res.h=(1LL*n*(n+1)%X*(2*n+1)%X*I6%X*(a/c)%X+(1LL*n*(n+1)/2)%X*(b/c)%X+t.h)%X;return res;
	}
	RI m=(1LL*a*n+b)/c;t=Get(m-1,c,c-b-1,a),res.f=(1LL*n*m%X-t.f+X)%X,//a<c且b<c
	res.g=((1LL*n*m%X*(m+1)%X-2LL*t.h%X-2LL*t.f%X-res.f)%X+X)%X,
	res.h=1LL*I2*((1LL*n*m%X*(n+1)%X-t.f-t.g)%X+X)%X;return res;
}
int main()
{
	RI Qt;S k;F.read(Qt);W(Qt--) F.read(n,a,b,c),k=Get(n,a,b,c),F.write(k.f,k.g,k.h);
	return F.clear(),0;
}