類歐幾里德算法
定義
類歐幾里德算法是洪華敦在 2016 年冬令營營員交流中提出的內容。
其本質可以理解為，使用一個類似輾轉相除法的方法來進行函數求和。
引入
設
\[
f(a,b,c,n)=\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor
\]
其中 \(a,b,c,n\) 是常數。需要一個 \(O(\log n)\) 的算法。
這個式子和我們以前見過的式子都長得不太一樣。帶向下取整的式子容易讓人想到數論分塊，然而數論分塊似乎不適用於這個求和。但是我們是可以做一些預處理的。
如果説 \(a\ge c\) 或者 \(b\ge c\) ，意味着可以將 \(a,b\) 對 \(c\) 取模以簡化問題：
\[
\begin{aligned}
f(a,b,c,n)&=\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor\\
&=\sum_{i=0}^n\left\lfloor
\frac{\left(\left\lfloor\frac{a}{c}\right\rfloor c+a\bmod c\right)i+\left(\left\lfloor\frac{b}{c}\right\rfloor c+b\bmod c\right)}{c}\right\rfloor\\
&=\frac{n(n+1)}{2}\left\lfloor\frac{a}{c}\right\rfloor+(n+1)\left\lfloor\frac{b}{c}\right\rfloor+
\sum_{i=0}^n\left\lfloor\frac{\left(a\bmod c\right)i+\left(b\bmod c\right)}{c}
\right\rfloor\\
&=\frac{n(n+1)}{2}\left\lfloor\frac{a}{c}\right\rfloor
+(n+1)\left\lfloor\frac{b}{c}\right\rfloor+f(a\bmod c,b\bmod c,c,n)
\end{aligned}
\]
那麼問題轉化為了 \(a<c,b<c\) 的情況。觀察式子，你發現只有 \(i\) 這一個變量。因此要推就只能從 \(i\) 下手。在推求和式子中有一個常見的技巧，就是條件與貢獻的放縮與轉化。具體地説，在原式 \(\displaystyle f(a,b,c,n)=\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor\) 中，\(0\le i\le n\) 是條件，而 \(\left\lfloor \dfrac{ai+b}{c} \right\rfloor\) 是對總和的貢獻。
要加快一個和式的計算過程，所有的方法都可以歸約為 貢獻合併計算 。但你發現這個式子的貢獻難以合併，怎麼辦？將貢獻與條件做轉化 得到另一個形式的和式。具體地，我們直接把原式的貢獻變成條件：
\[
\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor
=\sum_{i=0}^n\sum_{j=0}^{\left\lfloor \frac{ai+b}{c} \right\rfloor-1}1
\]
現在多了一個變量 \(j\) ，既然算 \(i\) 的貢獻不方便，我們就想辦法算 \(j\) 的貢獻。因此想辦法搞一個和 \(j\) 有關的貢獻式。這裏有另一個家喻户曉的變換方法，筆者概括為限制轉移。具體來説，在上面的和式中 \(n\) 限制 \(i\) 的上界，而 \(i\) 限制 \(j\) 的上界。為了搞 \(j\) ，就先把 j 放到貢獻的式子裏，於是我們交換一下 \(i,j\) 的求和算子，強制用 \(n\) 限制 \(j\) 的上界。
\[
=\sum_{j=0}^{\left\lfloor \frac{an+b}{c} \right\rfloor-1}\sum_{i=0}^n\left[j<\left\lfloor \frac{ai+b}{c} \right\rfloor\right]
\]
這樣做的目的是讓 \(j\) 擺脱 \(i\) 的限制，現在 \(i,j\) 都被 \(n\) 限制，而貢獻式看上去是一個條件，但是我們仍把它叫作貢獻式，再對貢獻式做變換後就可以改變 \(i,j\) 的限制關係。於是我們做一些放縮的處理。首先把向下取整的符號拿掉
\[
j<\left\lfloor \frac{ai+b}{c} \right\rfloor
\iff j+1\leq \left\lfloor \frac{ai+b}{c} \right\rfloor
\iff j+1\leq \frac{ai+b}{c}
\]
然後可以做一些變換
\[
j+1\leq \frac{ai+b}{c} \iff jc+c\le ai+b \iff jc+c-b-1< ai
\]
最後一步，向下取整得到：
\[
jc+c-b-1< ai\iff \left\lfloor\frac{jc+c-b-1}{a}\right\rfloor< i
\]
這一步的重要意義在於，我們可以把變量 \(i\) 消掉了！具體地，令 \(m=\left\lfloor \frac{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\frac{jc+c-b-1}{a}\right\rfloor \right]\\
&=\sum_{j=0}^{m-1}
(n-\left\lfloor\frac{jc+c-b-1}{a}\right\rfloor)\\
&=nm-f\left(c,c-b-1,a,m-1\right)
\end{aligned}
\]
這是一個遞歸的式子。並且你發現 \(a,c\) 分子分母換了位置，又可以重複上述過程。先取模，再遞歸。這就是一個輾轉相除的過程，這也是類歐幾里德算法的得名。
容易發現時間複雜度為 \(O(\log n)\) 。
擴展
理解了最基礎的類歐幾里德算法，我們再來思考以下兩個變種求和式：
\[
g(a,b,c,n)=\sum_{i=0}^ni\left\lfloor \frac{ai+b}{c} \right\rfloor
\]
\[
h(a,b,c,n)=\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor^2
\]
推導 g
我們先考慮 \(g\) ，類似地，首先取模：
\[
g(a,b,c,n)
=g(a\bmod c,b\bmod c,c,n)+\left\lfloor\frac{a}{c}\right\rfloor\frac{n(n+1)(2n+1)}{6}+\left\lfloor\frac{b}{c}\right\rfloor\frac{n(n+1)}{2}
\]
接下來考慮 \(a<c,b<c\) 的情況，令 \(m=\left\lfloor\frac{an+b}{c}\right\rfloor\) 。之後的過程比較簡略，因為方法和上文略同：
\[
\begin{aligned}
g(a,b,c,n)&=\sum_{i=0}^ni\left\lfloor \frac{ai+b}{c} \right\rfloor\\
&=\sum_{j=0}^{m-1}
\sum_{i=0}^n\left[j<\left\lfloor\frac{ai+b}{c}\right\rfloor\right]\cdot i
\end{aligned}
\]
這時我們設 \(t=\left\lfloor\frac{jc+c-b-1}{a}\right\rfloor\) ，可以得到
\[
\begin{aligned}
g(a,b,c,n)&=\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}\left[mn(n+1)-\sum_{j=0}^{m-1}t^2-\sum_{j=0}^{m-1}t\right]\\
&=\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\left\lfloor\frac{b}{c}\right\rfloor f(a\bmod c,b\bmod c,c,n)
+2\left\lfloor\frac{a}{c}\right\rfloor g(a\bmod c,b\bmod c,c,n)\\
&+\left\lfloor\frac{a}{c}\right\rfloor^2\frac{n(n+1)(2n+1)}{6}+\left\lfloor\frac{b}{c}\right\rfloor^2(n+1)
+\left\lfloor\frac{a}{c}\right\rfloor\left\lfloor\frac{b}{c}\right\rfloor n(n+1)
\end{aligned}
\]
考慮 \(a<c,b<c\) 的情況，\(m=\left\lfloor\dfrac{an+b}{c}\right\rfloor, t=\left\lfloor\dfrac{jc+c-b-1}{a}\right\rfloor\) .
我們發現這個平方不太好處理，於是可以這樣把它拆成兩部分：
\[
n^2=2\dfrac{n(n+1)}{2}-n=\left(2\sum_{i=0}^ni\right)-n
\]
這樣做的意義在於，添加變量 \(j\) 的時侯就只會變成一個求和算子，不會出現 \(\sum\times \sum\) 的形式：
\[
\begin{aligned}
h(a,b,c,n)&=\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor^2\\
&=\sum_{i=0}^n\left[\left(2\sum_{j=1}^{\left\lfloor \frac{ai+b}{c} \right\rfloor}j \right)-\left\lfloor\frac{ai+b}{c}\right\rfloor\right]\\
&=\left(2\sum_{i=0}^n\sum_{j=1}^{\left\lfloor \frac{ai+b}{c} \right\rfloor}j\right) -f(a,b,c,n)\\
\end{aligned}
\]
接下來考慮化簡前一部分：
\[
\begin{aligned}
\sum_{i=0}^n\sum_{j=1}^{\left\lfloor \frac{ai+b}{c} \right\rfloor}j&=\sum_{i=0}^n\sum_{j=0}^{\left\lfloor \frac{ai+b}{c} \right\rfloor-1}(j+1)\\
&=\sum_{j=0}^{m-1}(j+1)
\sum_{i=0}^n\left[j<\left\lfloor \frac{ai+b}{c} \right\rfloor\right]\\
&=\sum_{j=0}^{m-1}(j+1)\sum_{i=0}^n[i>t]\\
&=\sum_{j=0}^{m-1}(j+1)(n-t)\\
&=\frac{1}{2}nm(m+1)-\sum_{j=0}^{m-1}(j+1)\left\lfloor \frac{jc+c-b-1}{a} \right\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)
\]
在代碼實現的時侯，因為 \(3\) 個函數各有交錯遞歸，因此可以考慮三個一起整體遞歸，同步計算，否則有很多項會被多次計算。這樣實現的複雜度是 \(O(\log n)\) 的。
實現
模板題代碼實現 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
43
44
45
46
47
48
49
50
51
52
53
54 #include <bits/stdc++.h>
#define int long long
using namespace std ;
const int P = 998244353 ;
int i2 = 499122177 , i6 = 166374059 ;
struct data {
data () { f = g = h = 0 ; }
int f , g , h ;
}; // 三個函數打包
data calc ( int n , int a , int b , int c ) {
int ac = a / c , bc = b / c , m = ( a * n + b ) / c , n1 = n + 1 , n21 = n * 2 + 1 ;
data d ;
if ( a == 0 ) { // 迭代到最底層
d . f = bc * n1 % P ;
d . g = bc * n % P * n1 % P * i2 % P ;
d . h = bc * bc % P * n1 % P ;
return d ;
}
if ( a >= c || b >= c ) { // 取模
d . f = n * n1 % P * i2 % P * ac % P + bc * n1 % P ;
d . g = ac * n % P * n1 % P * n21 % P * i6 % P + bc * n % P * n1 % P * i2 % P ;
d . h = ac * ac % P * n % P * n1 % P * n21 % P * i6 % P +
bc * bc % P * n1 % P + ac * bc % P * n % P * n1 % P ;
d . f %= P , d . g %= P , d . h %= P ;
data e = calc ( n , a % c , b % c , c ); // 迭代
d . h += e . h + 2 * bc % P * e . f % P + 2 * ac % P * e . g % P ;
d . g += e . g , d . f += e . f ;
d . f %= P , d . g %= P , d . h %= P ;
return d ;
}
data e = calc ( m - 1 , c , c - b - 1 , a );
d . f = n * m % P - e . f , d . f = ( d . f % P + P ) % P ;
d . g = m * n % P * n1 % P - e . h - e . f , d . g = ( d . g * i2 % P + P ) % P ;
d . h = n * m % P * ( m + 1 ) % P - 2 * e . g - 2 * e . f - d . f ;
d . h = ( d . h % P + P ) % P ;
return d ;
}
int T , n , a , b , c ;
signed main () {
scanf ( "%lld" , & T );
while ( T -- ) {
scanf ( "%lld%lld%lld%lld" , & n , & a , & b , & c );
data ans = calc ( n , a , b , c );
printf ( "%lld %lld %lld \n " , ans . f , ans . h , ans . g );
}
return 0 ;
}
更新历史 在 GitHub 上编辑此页！ sshwy, FFjet CC BY-SA 4.0 和 SATA