[学习笔记] 类欧几里得算法
简介
类欧几里得算法其实是递归子问题巧妙运用的一个范例，主要用于计算下列柿子：
\[f(a,b,c,n)=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor
\]
\[g(a,b,c,n)=\sum_{i=0}^ni\lfloor\frac{ai+b}{c}\rfloor
\]
\[h(a,b,c,n)=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor^2
\]
这个东西看上去不是很好做，有一种特殊情况，就是可以把 \(a,b\) 模掉 \(c\) 之后递归下去。如果 \(a,b\) 比 \(c\) 小呢？那么我们就通过一些操作交换 \(a,c\) 的位置使得能像辗转相除那样递归下去，这样时间复杂度就是 \(O(\log n)\) 的了。
当然哪有说得这么简单，还是要推柿子。
我用的是传统方法，还有一个叫万能欧几里得的高科技，我太懒了不想学
f的推导
如果 \(a\geq c\) 或者 \(b\geq c\)，我们想把它取模之后递归下去：
\[\begin{aligned}
f(a,b,c,n)&=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor\\
&=\sum_{i=0}^n\lfloor\frac{(\lfloor\frac{a}{c}\rfloor c+a\bmod c)i+(\lfloor\frac{b}{c}\rfloor c+b\bmod c)}{c}\rfloor\\
&=\frac{n(n+1)}{2}\lfloor\frac{a}{c}\rfloor+(n+1)\lfloor\frac{b}{c}\rfloor+\sum_{i=0}^n\lfloor\frac{(a\bmod c)i+(b\bmod c)}{c}\rfloor\\
&=\frac{n(n+1)}{2}\lfloor\frac{a}{c}\rfloor+(n+1)\lfloor\frac{b}{c}\rfloor+f(a\bmod c,b\bmod c,c,n)\\
\end{aligned}
\]
剩下的就是 \(a<c,b<c\) 的情况了，直接推肯定是转化不动的，这里有一个神来之笔，我们引入一个新的枚举变量 \(j\)，然后交换 \(i,j\) 的求和顺序，把限制转化到 \(i\) 那里，就能找出递归子问题，看柿子理解这个思想：
\[\begin{aligned}
\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor&=\sum_{i=0}^n\sum_{j=0}^{\lfloor\frac{ai+b}{c}\rfloor-1}1\\
&=\sum_{j=0}^{\lfloor\frac{an+b}{c}\rfloor-1}\sum_{i=0}^n[j<\lfloor\frac{ai+b}{c}\rfloor]
\end{aligned}
\]
接下来要对里面的条件判断式动一下手脚了，记住我们的目的是把限制转化到 \(i\) 那里，所以要推 \(i\) 的不等式（我上面应该打等价符号，但是打不出来）：
\[j<\lfloor\frac{ai+b}{c}\rfloor\Rightarrow j+1\leq\lfloor\frac{ai+b}{c}\rfloor\Rightarrow j+1\leq\frac{ai+b}{c}
\]
\[\Rightarrow jc+c\leq ai+b\Rightarrow jc+c-b-1<ai\Rightarrow \lfloor\frac{jc+c-b-1}{a}\rfloor<i
\]
那么把这个限制条件代换到原式里面，设 \(m=\lfloor\frac{an+b}{c}\rfloor\)：
\[\begin{aligned}
f(a,b,c,n)&=\sum_{j=0}^{m-1}\sum_{i=0}^n[i>\lfloor\frac{jc+c-b-1}{a}\rfloor]\\
&=\sum_{j=0}^{m-1}n-\lfloor\frac{jc+c-b-1}{a}\rfloor\\
&=nm-f(c,c-b-1,a,m-1)
\end{aligned}
\]
你发现 \(a,c\) 交换了位置，所以复杂度就是辗转相除的 \(O(\log n)\)
g的推导
还是首先取模：
\[g(a,b,c,n)=g(a\bmod c,b\bmod c,c,n)+\lfloor\frac{a}{c}\rfloor\frac{n(n+1)(2n+1)}{6}+\lfloor\frac{b}{c}\rfloor\frac{n(n+1)}{2}
\]
设 \(m=\lfloor\frac{an+b}{c}\rfloor,t=\lfloor\frac{jc+c-b-1}{a}\rfloor\)，推导过程和上面大体是一样的：
\[\begin{aligned}
g(a,b,c,n)&=\sum_{i=0}^ni\lfloor\frac{ai+b}{c}\rfloor\\
&=\sum_{j=0}^{m-1}\sum_{i=0}^n[j<\lfloor\frac{ai+b}{c}\rfloor]\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的推导
虽然模的柿子看上去有点复杂，但其实直接暴力相乘展开就行了：
\[\begin{aligned}
h(a,b,c,n)&=h(a\bmod c,b\bmod c,c,n)\\
&+2\lfloor\frac{b}{c}\rfloor f(a\bmod c,b\bmod c,c,n)+2\lfloor\frac{a}{c}\rfloor g(a\bmod c,b\bmod c,c,n)\\
&+\lfloor\frac{a}{c}\rfloor^2\frac{n(n+1)(2n+1)}{6}+\lfloor\frac{b}{c}\rfloor^2(n+1)+\lfloor\frac{a}{c}\rfloor\lfloor\frac{b}{c}\rfloor n(n+1)
\end{aligned}
\]
这个平方要被表示成 \(j\) 的求和，要稍微操作一下：
\[n^2=(2\sum_{i=0}^ni)-n
\]
应用上面这个简单的原理来推导：
\[\begin{aligned}
h(a,b,c,n)&=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor^2=\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(a,b,c,n)
\end{aligned}
\]
化简前一个柿子就行了：
\[\begin{aligned}
\sum_{i=0}^n\sum_{j=1}^{\lfloor\frac{ai+b}{c}\rfloor}j&=\sum_{i=0}^n\sum_{j=0}^{\lfloor\frac{ai+b}{c}\rfloor-1}(j+1)\\
&=\sum_{j=0}^{m-1}(j+1)\sum_{i=0}^n[j<\lfloor\frac{ai+b}{c}\rfloor]\\
&=\sum_{j=0}^{m-1}(j+1)(n-t)\\
&=\frac{1}{2}nm(m+1)-\sum_{j=0}^{m-1}(j+1)\lfloor\frac{jc+c-b-1}{a}\rfloor\\
&=\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\n",t.f,t.h,t.g);
}
}
