简介
类欧几里得算法其实是递归子问题巧妙运用的一个范例，主要用于计算下列柿子：
[f(a,b,c,n)=sum_{i=0}^nlfloorfrac{ai+b}{c}
floor
]
[g(a,b,c,n)=sum_{i=0}^nilfloorfrac{ai+b}{c}
floor
]
[h(a,b,c,n)=sum_{i=0}^nlfloorfrac{ai+b}{c}
floor^2
]
这个东西看上去不是很好做，有一种特殊情况，就是可以把 (a,b) 模掉 (c) 之后递归下去。如果 (a,b) 比 (c) 小呢？那么我们就通过一些操作交换 (a,c) 的位置使得能像辗转相除那样递归下去，这样时间复杂度就是 (O(log n)) 的了。
当然哪有说得这么简单，还是要推柿子。
我用的是传统方法，还有一个叫万能欧几里得的高科技，我太懒了不想学
f的推导
如果 (ageq c) 或者 (bgeq c)，我们想把它取模之后递归下去：
[egin{aligned}
f(a,b,c,n)&=sum_{i=0}^nlfloorfrac{ai+b}{c}
floor\
&=sum_{i=0}^nlfloorfrac{(lfloorfrac{a}{c}
floor c+amod c)i+(lfloorfrac{b}{c}
floor c+bmod c)}{c}
floor\
&=frac{n(n+1)}{2}lfloorfrac{a}{c}
floor+(n+1)lfloorfrac{b}{c}
floor+sum_{i=0}^nlfloorfrac{(amod c)i+(bmod c)}{c}
floor\
&=frac{n(n+1)}{2}lfloorfrac{a}{c}
floor+(n+1)lfloorfrac{b}{c}
floor+f(amod c,bmod c,c,n)\
end{aligned}
]
剩下的就是 (a<c,b<c) 的情况了，直接推肯定是转化不动的，这里有一个神来之笔，我们引入一个新的枚举变量 (j)，然后交换 (i,j) 的求和顺序，把限制转化到 (i) 那里，就能找出递归子问题，看柿子理解这个思想：
[egin{aligned}
sum_{i=0}^nlfloorfrac{ai+b}{c}
floor&=sum_{i=0}^nsum_{j=0}^{lfloorfrac{ai+b}{c}
floor-1}1\
&=sum_{j=0}^{lfloorfrac{an+b}{c}
floor-1}sum_{i=0}^n[j<lfloorfrac{ai+b}{c}
floor]
end{aligned}
]
接下来要对里面的条件判断式动一下手脚了，记住我们的目的是把限制转化到 (i) 那里，所以要推 (i) 的不等式（我上面应该打等价符号，但是打不出来）：
[j<lfloorfrac{ai+b}{c}
floorRightarrow j+1leqlfloorfrac{ai+b}{c}
floorRightarrow j+1leqfrac{ai+b}{c}
]
[Rightarrow jc+cleq ai+bRightarrow jc+c-b-1<aiRightarrow lfloorfrac{jc+c-b-1}{a}
floor<i
]
那么把这个限制条件代换到原式里面，设 (m=lfloorfrac{an+b}{c} floor)：
[egin{aligned}
f(a,b,c,n)&=sum_{j=0}^{m-1}sum_{i=0}^n[i>lfloorfrac{jc+c-b-1}{a}
floor]\
&=sum_{j=0}^{m-1}n-lfloorfrac{jc+c-b-1}{a}
floor\
&=nm-f(c,c-b-1,a,m-1)
end{aligned}
]
你发现 (a,c) 交换了位置，所以复杂度就是辗转相除的 (O(log n))
g的推导
还是首先取模：
[g(a,b,c,n)=g(amod c,bmod c,c,n)+lfloorfrac{a}{c}
floorfrac{n(n+1)(2n+1)}{6}+lfloorfrac{b}{c}
floorfrac{n(n+1)}{2}
]
设 (m=lfloorfrac{an+b}{c} floor,t=lfloorfrac{jc+c-b-1}{a} floor)，推导过程和上面大体是一样的：
[egin{aligned}
g(a,b,c,n)&=sum_{i=0}^nilfloorfrac{ai+b}{c}
floor\
&=sum_{j=0}^{m-1}sum_{i=0}^n[j<lfloorfrac{ai+b}{c}
floor]cdot i\
&=sum_{j=0}^{m-1}sum_{i=0}^n[i>t]cdot i\
&=sum_{j=0}^{m-1}frac{1}{2}(t+n+1)(n-t)\
&=frac{1}{2}sum_{j=0}^{m-1}n^2+n-t-t^2\
&=frac{1}{2}[mn(n+1)-h(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1)]
end{aligned}
]
h的推导
虽然模的柿子看上去有点复杂，但其实直接暴力相乘展开就行了：
[egin{aligned}
h(a,b,c,n)&=h(amod c,bmod c,c,n)\
&+2lfloorfrac{b}{c}
floor f(amod c,bmod c,c,n)+2lfloorfrac{a}{c}
floor g(amod c,bmod c,c,n)\
&+lfloorfrac{a}{c}
floor^2frac{n(n+1)(2n+1)}{6}+lfloorfrac{b}{c}
floor^2(n+1)+lfloorfrac{a}{c}
floorlfloorfrac{b}{c}
floor n(n+1)
end{aligned}
]
这个平方要被表示成 (j) 的求和，要稍微操作一下：
[n^2=(2sum_{i=0}^ni)-n
]
应用上面这个简单的原理来推导：
[egin{aligned}
h(a,b,c,n)&=sum_{i=0}^nlfloorfrac{ai+b}{c}
floor^2=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(a,b,c,n)
end{aligned}
]
化简前一个柿子就行了：
[egin{aligned}
sum_{i=0}^nsum_{j=1}^{lfloorfrac{ai+b}{c}
floor}j&=sum_{i=0}^nsum_{j=0}^{lfloorfrac{ai+b}{c}
floor-1}(j+1)\
&=sum_{j=0}^{m-1}(j+1)sum_{i=0}^n[j<lfloorfrac{ai+b}{c}
floor]\
&=sum_{j=0}^{m-1}(j+1)(n-t)\
&=frac{1}{2}nm(m+1)-sum_{j=0}^{m-1}(j+1)lfloorfrac{jc+c-b-1}{a}
floor\
&=frac{1}{2}nm(m+1)-g(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1)
end{aligned}
]
把上面的推导综合一下，得到了我们想要的东西：
[h(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)
]
例题
先看一道模板：类欧几里得算法
为了方便写，可以直接三合一，一个递归解决问题。
#include <cstdio>
#define int long long
const int MOD = 998244353;
const int inv2 = (MOD+1)/2;
const int inv6 = (MOD+1)/6;
int read()
{
int x=0,f=1;char c;
while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
return x*f;
}
int T,n,a,b,c;
struct node
{
int f,g,h;
node() {f=g=h=0;}
};
node work(int n,int a,int b,int c)
{
int ac=a/c,bc=b/c,n1=n+1,n2=(2*n+1)%MOD,m=(a*n+b)/c;
node d;
if(!a)
{
d.f=bc*n1%MOD;
d.g=bc*n%MOD*n1%MOD*inv2%MOD;
d.h=bc*bc%MOD*n1%MOD;
return d;
}
if(a>=c || b>=c)
{
d.f=n*n1%MOD*inv2%MOD*ac%MOD+n1*bc%MOD;
d.g=ac*n%MOD*n1%MOD*n2%MOD*inv6%MOD+bc*n%MOD*n1%MOD*inv2%MOD;
d.h=ac*ac%MOD*n%MOD*n1%MOD*n2%MOD*inv6%MOD
+bc*bc%MOD*n1%MOD+ac*bc%MOD*n%MOD*n1%MOD;
d.f%=MOD;d.g%=MOD;d.h%=MOD;
node e=work(n,a%c,b%c,c);
d.f=(d.f+e.f)%MOD;
d.g=(d.g+e.g)%MOD;
d.h=(d.h+e.h+2*bc%MOD*e.f+2*ac%MOD*e.g)%MOD;
return d;
}
node e=work(m-1,c,c-b-1,a);
d.f=(n*m-e.f)%MOD;
d.g=(m*n%MOD*n1-e.h-e.f)%MOD*inv2%MOD;
d.h=(n*m%MOD*(m+1)-2*e.g-2*e.f-d.f)%MOD;
d.f=(d.f+MOD)%MOD;d.g=(d.g+MOD)%MOD;d.h=(d.h+MOD)%MOD;
return d;
}
signed main()
{
T=read();
while(T--)
{
n=read();a=read();b=read();c=read();
node t=work(n,a,b,c);
printf("%lld %lld %lld
",t.f,t.h,t.g);
}
}