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

前言

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

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

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

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

一道板子题

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

题目大意就是让我们求(sum_{i=0}^nlfloorfrac{ai+b}c floor,sum_{i=0}^nlfloorfrac{ai+b}c floor^2,sum_{i=0}^nilfloorfrac{ai+b}c floor)，其中(n,a,b,cle10^6)。

考虑我们设：

[f(n,a,b,c)=sum_{i=0}^nlfloorfrac{ai+b}c floor ]

[g(n,a,b,c)=sum_{i=0}^nlfloorfrac{ai+b}c floor^2 ]

[h(n,a,b,c)=sum_{i=0}^nilfloorfrac{ai+b}c floor ]

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

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

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

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

当(a=0)时：

[f(n,a,b,c)=sum_{i=0}^nlfloorfrac bc floor=(n+1) imeslfloorfrac bc floor ]

当(age c)或(bge c)时：

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

[lfloorfrac {ai+b}c floor=lfloorfrac{(a\%c)i+(b\%c)}c floor+lfloorfrac ac floor i+lfloorfrac bc floor ]

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

[egin{align}f(n,a,b,c)&=sum_{i=0}^nlfloorfrac{(a\%c)i+(b\%c)}c floor+lfloorfrac ac floor i+lfloorfrac bc floor\&=f(n,a\%c,b\%c,c)+lfloorfrac ac floor imesfrac{n(n+1)}2+lfloorfrac bc floor imes(n+1)end{align} ]

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

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

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

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

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

[egin{align}f(n,a,b,c) & =sum_{i=0}^nsum_{j=1}^m[lfloorfrac{ai+b}c floorge j]\ & =sum_{i=0}^nsum_{j=0}^{m-1}[lfloorfrac{ai+b}c floorge j+1]\ & =sum_{i=0}^nsum_{j=0}^{m-1}[ai+bge cj+c]\ & =sum_{i=0}^nsum_{j=0}^{m-1}[ai+b+1>cj+c]\ & =sum_{i=0}^nsum_{j=0}^{m-1}[i>lfloorfrac{cj+c-b-1}a floor]\ & =sum_{j=0}^{m-1}(n-lfloorfrac{cj+c-b-1}a floor)\ & =n imes 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}^nlfloorfrac bc floor^2=(n+1) imeslfloorfrac bc floor^2 ]

当(age c)或(bge c)时：

[egin{align}g(n,a,b,c)&=sum_{i=0}^n(lfloorfrac{(a\%c)i+(b\%c)}c floor+lfloorfrac ac floor i+lfloorfrac bc floor)^2\&=sum_{i=0}^n(lfloorfrac{(a\%c)i+(b\%c)}c floor^2+lfloorfrac ac floor^2i^2+lfloorfrac bc floor^2+2ilfloorfrac ac floor lfloorfrac{(a\%c)i+(b\%c)}c floor+2lfloorfrac bc floorlfloorfrac{(a\%c)i+(b\%c)}c floor+2ilfloorfrac ac floorlfloorfrac bc floor)\&=g(n,a\%c,b\%c,c)+frac{n(n+1)(2n+1)}6lfloorfrac ac floor^2+(n+1)lfloorfrac bc floor^2+2lfloorfrac ac floor h(n,a\%c,b\%c,c)+2lfloorfrac bc floor f(n,a\%c,b\%c,c)+n(n+1)lfloorfrac ac floorlfloorfrac bc floorend{align} ]

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

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

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

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

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

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

[egin{align}g(n,a,b,c)&=sum_{i=0}^n(2sum_{j=1}^{lfloorfrac{ai+b}c floor}j-lfloorfrac{ai+b}c floor)\&=2sum_{i=0}^nsum_{j=1}^{lfloorfrac{ai+b}c floor}j-f(n,a,b,c)end{align} ]

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

[egin{align}2sum_{i=0}^nsum_{j=1}^{lfloorfrac{ai+b}c floor}j&=2sum_{i=0}^nsum_{j=1}^mj[jlelfloorfrac{ai+b}c floor]\&=2sum_{i=0}^nsum_{j=0}^{m-1}(j+1)[i>lfloorfrac{cj+c-b-1}a floor]\&=2sum_{j=0}^{m-1}(j+1)(n-lfloorfrac{cj+c-b-1}a floor)\&=2sum_{j=0}^{m-1}(n(j+1)-jlfloorfrac{cj+c-b-1}a floor-lfloorfrac{cj+c-b-1}a floor)\&=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}^nilfloorfrac bc floor=frac{n(n+1)}2 imeslfloorfrac bc floor ]

当(age c)或(bge c)时：

[egin{align}h(n,a,b,c)&=sum_{i=0}^ni(lfloorfrac{(a\%c)i+(b\%c)}c floor+lfloorfrac ac floor i+lfloorfrac bc floor)\&=sum_{i=0}^n(ilfloorfrac{(a\%c)i+(b\%c)}c floor+lfloorfrac ac floor i^2+lfloorfrac bc floor i)\&=h(n,a\%c,b\%c,c)+frac{n(n+1)(2n+1)}6lfloorfrac ac floor+frac{n(n+1)}2lfloorfrac bc floorend{align} ]

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

[egin{align}h(n,a,b,c)&=sum_{i=0}^nisum_{j=1}^m[lfloorfrac{ai+b}c floorge j]\&=sum_{i=0}^nisum_{j=0}^{m-1}[i>lfloorfrac{cj+c-b-1}a floor]\&=sum_{j=0}^{m-1}sum_{i=0}^ni[i>lfloorfrac{cj+c-b-1}a floor]\&=sum_{j=0}^{m-1}frac{(lfloorfrac{cj+c-b-1}a floor+1+n)(n-lfloorfrac{cj+c-b-1}a floor)}2\&=frac12sum_{j=0}^{m-1}(n(n+1)-lfloorfrac{cj+c-b-1}a floor-lfloorfrac{cj+c-b-1}a floor^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('
');}
		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;
}