类欧几里得算法

本题以 P5170 【模板】类欧几里得算法 为例，深入讲解类欧几里得算法的本质。

普通类欧几里得算法

首先，设 $f(n,a,b,c)=\sum_{i=0}^n \lfloor \frac {ai+b}{c} \rfloor$，我们想在 $\log$ 时间范围内求出这个函数的值。

首先考虑把 $a \ge c$ 或者 $b \ge c$ 的情况规约到 $a<c$ 并且 $b<c$ 的情况，容易发现，如果 $a>c$，式子可以化成这样：
$$
\sum_{i=0}^n \lfloor \frac {(a \bmod c)i+b}{c} \rfloor + \sum_{i=0}^n \lfloor \frac ac \rfloor i
$$
后面的一堆可以用等差数列的求和公式计算，如果 $b>c$，式子可以化成这样：
$$
\sum_{i=0}^n \lfloor \frac {ai+b \bmod c}{c} \rfloor+\sum_{i=0}^n \lfloor \frac bc \rfloor
$$
后面一堆可以直接计算，于是我们只需要考虑 $a<c$ 并且 $b<c$ 的情况。

我们来推一下式子：
$$
\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}
$$

再设 $m=\frac {an+b}{c}$，则原式变为 $\sum_{j=0}^{m-1} \sum_{i=0}^n [j<\lfloor \frac {ai+b}{c} \rfloor]$。

考虑对艾佛森括号内的式子进行变换，因为两边都是整数，所以可以写成 $[j+1 \le \lfloor \frac {ai+b}{c} \rfloor]$，又因为 $c$ 一定是正整数，根据带余除法，设 $ai+b=kc+q(0 \le q<c)$，则 $\lfloor \frac {ai+b}{c} \rfloor=k$，所以两边乘 $c$ 是没有问题的，即 $[cj+c \le ai+b]$，就相当于若 $j+1 \le k$，则此式子一定成立，否则此式子一定不成立。

这时我们继续变换之后下取整得到 $cj+c \le ai+b \rightarrow cj+c-b \le ai \rightarrow cj+c-b-1<ai \rightarrow \lfloor\frac {cj+c-b-1}{a}\rfloor<i$。

代回原式，得 $\sum_{j=0}^{m-1} \sum_{i=0}^n [\lfloor \frac {cj+c-b-1}{a}\rfloor<i]$，这个函数后面一个 $\sum$ 的值我们可以简单地计算得出，所以原式变为：
$$
\sum_{j=0}^{m-1} n-\lfloor \frac {cj+c-b-1}{a} \rfloor=nm-\sum_{j=0}^{m-1} \lfloor \frac {cj+c-b-1}{a} \rfloor
$$
于是我们可以把函数写成一个常数加上另外一个函数的和，如下：
$$
f(n,a,b,c)=nm-f(m-1,c,c-b-1,a)
$$
因为每次我们都会让 $a \gets a \bmod c$，所以这个函数的计算是 $\log$ 级别的，时间复杂度分析很像欧几里得算法，故称为类欧几里得算法。

边界条件：$a=0$ 的时候函数等于 $\sum_{i=0}^n \lfloor \frac bc \rfloor$，直接计算即可，这个部分代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

#include<bits/stdc++.h>
#define mod 998244353
#define ll long long
using namespace std;
ll T,n,a,b,c,inv;
ll solve(ll n,ll a,ll b,ll c){
	if(a==0) return (b/c)*(n+1)%mod;
	if(a>=c||b>=c) return (solve(n,a%c,b%c,c)+(a/c)*n%mod*(n+1)%mod*inv%mod+(b/c)*(n+1)%mod)%mod;
	ll m = (a*n+b)/c;
	return (n*m%mod-solve(m-1,c,c-b-1,a)+mod)%mod;
}
int main(){
	inv = (mod+1)/2;
	ios::sync_with_stdio(false);
	cin>>T;
	while(T--){
		cin>>n>>a>>b>>c;
		cout<<solve(n,a,b,c)<<endl;
	}
	return 0;
} 

平方类欧几里得算法

设 $h(n,a,b,c)=\sum_{i=0}^n {\lfloor \frac {ai+b}{c} \rfloor}^2$，我们仍然需要在 $\log$ 时间范围内求出这个式子的值。

首先像普通类欧那样让 $a,b \gets a \bmod c,b \bmod c$。

考虑化一下，设 $a=pc+p’(0\le p’<c),b=qc+q’(0 \le q’ < c)$，则：
$$
\begin{aligned}
h(n,a,b,c)&=\sum_{i=0}^n {\lfloor \frac {ai+b}{c} \rfloor}^2 \\
&=\sum_{i=0}^n {(\lfloor \frac {p’i+q’}{c} \rfloor+pi+q)}^2 \\
&=\sum_{i=0}^n {\lfloor \frac {p’i+q’}{c} \rfloor}^2+\sum_{i=0}^n p^2i^2+\sum_{i=0}^n q^2+2\sum_{i=0}^n pqi+2p\sum_{i=0}^n i\lfloor \frac{p’i+q’}{c} \rfloor+2q\sum_{i=0}^n \lfloor \frac{p’i+q’}{c} \rfloor
\end{aligned}
$$
其中第五个项是等会要讲的带权类欧几里得算法，暂且设 $g(n,a,b,c)=\sum_{i=0}^n i\lfloor \frac {ai+b}{c} \rfloor$，则原式变为：
$$
h(n,a,b,c)=h(n,p’,q’,c)+2qf(n,p’,q’,c)+2pg(n,p’,q’,c)+p^2\frac {n(n+1)(2n+1)}{6}+q^2(n+1)+pqn(n+1)
$$
接下来我们只需要计算 $h(n,p’,q’,c)$ 的值，因为 $p’,q’<c$，我们可以继续推导。
$$
h(n,a,b,c)=\sum_{i=0}^n {\lfloor \frac {ai+b}{c} \rfloor}^2
$$
首先 $t^2=2 \times \frac {t(t+1)}{2}-t$，所以：
$$
h(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)
$$
后面可以直接计算，看前面的：
$$
\begin{aligned}
\sum_{i=0}^n \sum_{j=1}^{\lfloor \frac {ai+b}{c} \rfloor} j &= \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) \sum_{i=0}^n [\lfloor \frac {cj+c-b-1}{a} \rfloor <i] \\
&= \sum_{j=0}^{m-1} (j+1)(n-\lfloor \frac {cj+c-b-1}{a} \rfloor) \\
&= nm+n \frac {m(m-1)}{2}-f(m-1,c,c-b-1,a)-g(m-1,c,c-b-1,a)
\end{aligned}
$$
带回原式，得：
$$
h(n,a,b,c)=nm(m+1)-2f(m-1,c,c-b-1,a)-2g(m-1,c,c-b-1,a)-f(n,a,b,c)
$$

带权类欧几里得算法

在 $\log$ 时间复杂度内求解 $g(n,a,b,c)=\sum_{i=0}^n i \lfloor \frac {ai+b}{c} \rfloor$。

先化掉 $a \ge c,b \ge c$ 的情况，设 $a=pc+p’,b=qc+q’(0 \le p’,q’<c)$，则：
$$
\sum_{i=0}^n i\lfloor \frac {ai+b}{c} \rfloor= \sum_{i=0}^n i(\lfloor \frac {p’i+q’}{c} \rfloor+pi+q)=\sum_{i=0}^n i \lfloor \frac {p’i+q’}{c} \rfloor+\sum_{i=0}^n pi^2+\sum_{i=0}^n qi
$$
后面两个和式可以 $O(1)$ 计算。

设 $m=\lfloor \frac {an+b}{c} \rfloor$，还是化一下式子：
$$
\begin{aligned}
\sum_{i=0}^n i \lfloor \frac {ai+b}{c} \rfloor &= \sum_{i=0}^n i\sum_{j=0}^{\lfloor \frac {ai+b}{c} \rfloor-1} 1 \\
&= \sum_{j=0}^{m-1} \sum_{i=0}^n i[j<\lfloor \frac {ai+b}{c} \rfloor]
\end{aligned}
$$
采用之前的方法，化掉艾佛森括号：
$$
\begin{aligned}
\sum_{j=0}^{m-1} \sum_{i=0}^n i[j<\lfloor \frac {ai+b}{c} \rfloor] &= \sum_{j=0}^{m-1} \sum_{i=0}^n i[\lfloor \frac {cj+c-b-1}{a} \rfloor<i] \\
t \gets \lfloor \frac {cj+c-b-1}{a} \rfloor \\
&= \sum_{j=0}^{m-1} \frac {(t+1+n)(n-t)}{2} \\
&= \frac 12(\sum_{j=0}^{m-1}(n^2+n)-\sum_{j=0}^{m-1} (t^2+t)) \\
&= \frac 12 ((n^2+n)m-\sum_{j=0}^{m-1} t^2-\sum_{j=0}^{m-1}t)
\end{aligned}
$$
因此，和式都可以递归进行计算，于是：
$$
g(n,a,b,c)=\frac 12(n(n+1)m-h(m-1,c,c-b-1,a)-f(m-1,c,c-b-1,a))
$$

其他事项

单个函数计算显然会超时，于是我们考虑用一个结构体存下 $n,a,b,c$ 对应的 $f,g,h$ 值，这样就可以在 $\log$ 时间内快速计算三个函数的值了，代码中的 $h$ 和 $g$ 函数交换了含义，读的时候要注意一下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

#include<bits/stdc++.h>
#define mod 998244353
#define ll long long
using namespace std;
struct node{ll f,g,h;};
ll T,n,a,b,c,inv2,inv6;
inline ll qmi(ll a,ll b,ll p){
	ll res = 1%p,t=a;
	while(b){
		if(b&1) res=res*t%p;
		t=t*t%p;
		b>>=1;
	}
	return res;
}
node solve(ll n,ll a,ll b,ll c){
	if(a==0) return (node){(b/c)*(n+1)%mod,(b/c)%mod*(b/c)%mod*(n+1)%mod,(b/c)*n%mod*(n+1)%mod*inv2%mod};
	if(a>=c||b>=c){
		node temp = solve(n,a%c,b%c,c);
		temp.g = (temp.g+2*(b/c)%mod*temp.f%mod+2*(a/c)%mod*temp.h%mod+(a/c)%mod*(a/c)%mod*n%mod*(n+1)%mod*(2*n+1)%mod*inv6%mod+(b/c)*(b/c)%mod*(n+1)%mod+(a/c)*(b/c)%mod*n%mod*(n+1)%mod)%mod;
		temp.h = (temp.h+(a/c)*n%mod*(n+1)%mod*(2*n+1)%mod*inv6%mod+(b/c)*n%mod*(n+1)%mod*inv2%mod)%mod;
		temp.f = (temp.f+(a/c)*n%mod*(n+1)%mod*inv2%mod+(b/c)*(n+1)%mod)%mod;
		return temp;
	}
	ll m = (a*n+b)/c;
	node temp = solve(m-1,c,c-b-1,a),ans;
	ans.h = ((n*m%mod*(n+1)%mod-temp.g-temp.f)*inv2%mod+mod)%mod;
	ans.f = (n*m%mod-temp.f+mod)%mod;
	ans.g = ((n*m%mod*(m+1)%mod-2*temp.h%mod-2*temp.f%mod-ans.f)%mod+mod)%mod;
	return ans;
}
int main(){
	inv2 = qmi(2,mod-2,mod),inv6 = qmi(6,mod-2,mod);
	ios::sync_with_stdio(false);
	cin>>T;
	while(T--){
		cin>>n>>a>>b>>c;
		node ans = solve(n,a,b,c);
		cout<<ans.f<<" "<<ans.g<<" "<<ans.h<<endl;
	}
	return 0;
} 

万能欧几里得算法

待补。