【学习笔记】类欧几里得算法
概述
主要是求以下三个式子：
\[f(a,b,c,n)=\sum_{i=0}^n \left\lfloor\dfrac{ai+b}{c}\right\rfloor
\]
\[g(a,b,c,n)=\sum_{i=0}^n i\left\lfloor\dfrac{ai+b}{c}\right\rfloor
\]
\[h(a,b,c,n)=\sum_{i=0}^n \left\lfloor\dfrac{ai+b}{c}\right\rfloor^2
\]
求 \(f(a,b,c,n)\)
首先规范成 \(a<c,b<c\) 的情况，可以推得：
\[\begin{aligned}
f(a,b,c,n)&=\sum_{i=0}^n \left\lfloor\dfrac{(\left\lfloor\frac{a}{c}\right\rfloor c+a\bmod c)i+(\left\lfloor\frac{b}{c}\right\rfloor c+b\bmod c)}{c}\right\rfloor\\
&=\left\lfloor\dfrac{a}{c}\right\rfloor\dfrac{n(n+1)}{2}+\left\lfloor\dfrac{b}{c}\right\rfloor(n+1)+\sum_{i=0}^n \left\lfloor\dfrac{(a\bmod c)i+(b\bmod c)}{c}\right\rfloor\\
&=\left\lfloor\dfrac{a}{c}\right\rfloor\dfrac{n(n+1)}{2}+\left\lfloor\dfrac{b}{c}\right\rfloor(n+1)+f(a\bmod c,b\bmod c,c,n)
\end{aligned}
\]
之后再推导：
\[\begin{aligned}
f(a,b,c,n)&=\sum_{i=0}^n \left\lfloor\dfrac{ai+b}{c}\right\rfloor\\
&=\sum_{i=0}^n\sum_{j=0}^{\left\lfloor\frac{ai+b}{c}\right\rfloor-1} 1\\
&=\sum_{j=0}^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}\sum_{i=0}^n \left[j<\left\lfloor\dfrac{ai+b}{c}\right\rfloor\right]\\
\end{aligned}
\]
可以经过一系列变换：
\[j<\left\lfloor\dfrac{ai+b}{c}\right\rfloor\Leftrightarrow j+1\le \left\lfloor\dfrac{ai+b}{c}\right\rfloor\Leftrightarrow j+1\le \dfrac{ai+b}{c}
\]
整理可以得到：
\[jc+c-b\le ai\Leftrightarrow jc+c-b-1<ai\Leftrightarrow \left\lfloor\dfrac{jc+c-b-1}{a}\right\rfloor<i
\]
设 \(m=\left\lfloor\dfrac{an+b}{c}\right\rfloor\)，可以得到：
\[\begin{aligned}
f(a,b,c,n)&=\sum_{j=0}^{m-1} \sum_{i=0}^n \left[i>\left\lfloor\dfrac{jc+c-b-1}{a}\right\rfloor\right]\\
&=\sum_{j=0}^{m-1} \left(n-\left\lfloor\dfrac{jc+c-b-1}{a}\right\rfloor\right)\\
&=mn-f(c,c-b-1,a,m-1)
\end{aligned}
\]
注意到递归过程 \(a,c\) 交换，类似欧几里得算法（因此叫类欧几里得算法），复杂度是 \(O(\log n)\)。
上面的推导过程重点在：
-
设立第二维并交换求和号。
-
把艾弗森括号内的式子放缩并变形，方便化简。
求 \(g(a,b,c,n)\)
设 \(m=\left\lfloor\dfrac{an+b}{c}\right\rfloor\)，\(k=\left\lfloor\dfrac{jc+c-b-1}{a}\right\rfloor\)，类比过程。
先是取模：
\[g(a,b,c,n)=\left\lfloor\dfrac{a}{c}\right\rfloor\dfrac{n(n+1)(2n+1)}{6}+\left\lfloor\dfrac{b}{c}\right\rfloor \dfrac{n(n+1)}{2}+g(a\bmod c,b\bmod c,c,n)
\]
之后是：
\[\begin{aligned}
g(a,b,c,n)&=\sum_{i=0}^n i\left\lfloor\dfrac{ai+b}{c}\right\rfloor\\
&=\sum_{i=0}^n i\sum_{j=0}^{\left\lfloor\frac{ai+b}{c}\right\rfloor-1} 1\\
&=\sum_{j=0}^{m-1} \sum_{i=0}^n i[i>k]\\
&=\sum_{j=0}^{m-1} \dfrac{(k+1+n)(n-k)}{2}\\
&=\dfrac{1}{2}\sum_{j=0}^{m-1} n(n+1)-k-k^2\\
&=\dfrac{1}{2}(mn(n+1)-f(c,c-b-1,a,m-1)-h(c,c-b-1,a,m-1))
\end{aligned}
\]
求 \(h(a,b,c,n)\)
直接上过程了：
\[\begin{aligned}
h(a,b,c,n)=&\left\lfloor\dfrac{a}{c}\right\rfloor^2 \dfrac{n(n+1)(2n+1)}{6}+\left\lfloor\dfrac{b}{c}\right\rfloor^2 (n+1)+h(a\bmod c,b\bmod c,c,n)\\&+2\left\lfloor\dfrac{a}{c}\right\rfloor g(a\bmod c,b\bmod c,c,n)+2\left\lfloor\dfrac{b}{c}\right\rfloor f(a\bmod c,b\bmod c,c,n)+\left\lfloor\dfrac{a}{c}\right\rfloor\left\lfloor\dfrac{b}{c}\right\rfloor(n+1)\end{aligned}
\]
之后略微有变化，由于式子中带着平方，为了避免出现 \(\sum\times \sum\) 之类的情况，使用 \(n^2=2\sum_{i=0}^n i-n\) 来改写。
于是：
\[\begin{aligned}
h(a,b,c,n)&=\sum_{i=0}^n \left\lfloor\dfrac{ai+b}{c}\right\rfloor^2\\
&=\sum_{i=0}^n 2\sum_{j=0}^{\left\lfloor\frac{ai+b}{c}\right\rfloor} j-\left\lfloor\dfrac{ai+b}{c}\right\rfloor\\
&=2\sum_{i=0}^n \sum_{j=0}^{\left\lfloor\frac{ai+b}{c}\right\rfloor-1}(j+1)-f(a,b,c,n)\\
&=2\sum_{j=0}^{m-1} (j+1)\sum_{i=0}^n [i>k]-f(a,b,c,n)\\
&=2\sum_{j=0}^{m-1} (j+1)(n-k)-f(a,b,c,n)\\
&=2\sum_{j=0}^{m-1} n(j+1)-k-jk-f(a,b,c,n)\\
&=2\left(\dfrac{m(m+1)n}{2}-f(c,c-b-1,a,m-1)-g(c,c-b-1,a,m-1)\right)-f(a,b,c,n)\\
&=m(m+1)n-2f(c,c-b-1,a,m-1)-2g(c,c-b-1,a,m-1)-f(a,b,c,n)
\end{aligned}
\]
注意到 \((a,b,c,n)\) 的答案只和 \((c,c-b-1,a,m-1)\) 有关，所以只递归一次即可。
点击查看代码
struct Data{
ll f,g,h;
Data()=default;
Data(ll f_,ll g_,ll h_):f(f_),g(g_),h(h_){}
};
Data solve(int a,int b,int c,int n){
Data now=Data(0,0,0),res=Data(0,0,0),last=Data(0,0,0);
ll S0=(n+1)%mod,S1=1ll*n*(n+1)%mod*inv2%mod,S2=1ll*n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
int d1=a/c,d2=b/c;
a%=c,b%=c;
int m=(1ll*a*n+b)/c;
if(a) last=solve(c,c-b-1,a,m-1);
now.f=(1ll*n*m%mod-last.f+mod)%mod;
res.f=(S1*d1%mod+S0*d2%mod+now.f)%mod;
now.g=(1ll*m*n%mod*(n+1)%mod-last.f+mod-last.h+mod)%mod*inv2%mod;
res.g=(S2*d1%mod+S1*d2%mod+now.g)%mod;
now.h=(1ll*n*m%mod*(m+1)%mod-2*last.f%mod+mod-2*last.g%mod+mod-now.f+mod)%mod;
res.h=(S2*d1%mod*d1%mod+S0*d2%mod*d2%mod+now.h+2ll*d1*now.g%mod+2ll*d2*now.f%mod+2*S1*d1%mod*d2%mod)%mod;
return res;
}
参考资料
- OI Wiki
