類歐幾里得演算法

前幾日做了名為"It's a Mod, Mod, Mod, Mod World"(膜世界)的題

給定\(t\)組\(p,n,q\),求解

\(\sum_{i=1}^n [(pi)\,mod\,q]\)，其中\(1 \le p,n,q \le 10^6\),\(1 \le t \le 10^4\),限時\(1s\)

注意，這個\(mod\)是求和號里面的，不能對求和整體取模

這個題簡直是開幕雷擊，完全不知道怎么下手，有興趣的讀者可以思考半個小時，后來我百度了題目名，才知道要用到“類歐幾里得演算法”，
歐幾里德演算法又稱輾轉相除法，用于計算最大公因數（我們把a,b的最大公因數簡稱\(gcd(a,b)\)），計算\(gcd(a,b)\)的復雜度為\(O( log(max(a,b)) )\)，非常快，
同樣的，在“類歐幾里得演算法”當中，也用到了類似輾轉相除的操作，那么最終的復雜度也是log級別的，解決上述問題非常合適

最簡單的"類歐幾里得演算法"可以用于求解下列函式的值
\(f(a,b,c,n)=\sum_{i=1}^n \lfloor\frac{ai+b}{c} \rfloor\)
其中，a,b,c,n都是給定的整數,時間復雜度為\(O( log(max(a,c)) )\)
下面給出演算法的原理以及步驟:

首先，把這個式子的引數化小一點

\(f(a,b,c,n)=\sum_{i=1}^n \lfloor\frac{ai+b}{c}\rfloor =\sum_{i=1}^n \lfloor\frac{c\lfloor\frac{a}{c}\rfloor i+(a\,mod\,c)i+b\,mod\,c+c\lfloor \frac{b}{c} \rfloor}{c} \rfloor =\sum_{i=1}^n \lfloor\frac{(a\,mod\,c)i+b\,mod\,c}{c} \rfloor+\sum_{i=1}^n \lfloor \frac{a}{c} \rfloor i+ \lfloor \frac{b}{c} \rfloor\)
第三個等號的依據是
當\(y\)是一個整數，\(x\)是一個實數，則有
\(\lfloor x+y \rfloor =\lfloor x \rfloor + y\)
其中，第二個求和式可以用等引數列公式計算，結果為
\(\frac{1}{2} \lfloor \frac{a}{c} \rfloor n^2+\frac{1}{2}(\lfloor \frac{a}{c} \rfloor+2\lfloor\frac{b}{c}\rfloor)n\)

用\(f\)的形式簡寫，可以得到
等式1：\(f(a,b,c,n)=f(a\,mod\,c,b\,mod\,c,c,n)+\frac{1}{2} \lfloor \frac{a}{c} \rfloor n^2+\frac{1}{2}(\lfloor \frac{a}{c} \rfloor+2\lfloor\frac{b}{c}\rfloor)n\)
以上操作讓后得到的\(a\,mod\,c\),\(b\,mod\,c\)均小于c，所以在下面的討論中，不妨假設\(a,b<c\)，這是因為只要我們能求出滿足\(a,b<c\)的\(f\)函式的值，就能求出所有的\(f\)函式的值

下面分兩種情況研究\(f\)函式

1.特殊情況:

\(a=0\)

\(f(a,b,c,n)=\sum_{i=1}^n \lfloor\frac{b}{c} \rfloor=\lfloor\frac{b}{c} \rfloor n\)

這種情況對應著遞回求解的終點，在下面的一般情況的推導中，會介紹遞回演算法的遞推方程

2.一般情況：
\(a>0\)
\(\sum_{i=1}^n \lfloor \frac{ai+b}{c} \rfloor\)
\(=\sum_{i=1}^n \sum_{j=1}^{\lfloor\frac{ai+b}{c}\rfloor} 1\)
\(=\sum_{i=1}^n \sum_{j=1}^{\lfloor\frac{an+b}{c}\rfloor} (j\le \lfloor\frac{ai+b}{c}\rfloor)\)
其中\(x\le y\)為邏輯運算式，邏輯運算式為真時值為1，為假時值為0
第二層求和的上限從\(\lfloor\frac{ai+b}{c}\rfloor\)變成了\(\lfloor\frac{an+b}{c}\rfloor\)
這是為了讓兩層求和上限不相關，可以交換兩個求和號的順序
\(=\sum_{j=1}^{\lfloor\frac{an+b}{c}\rfloor} \sum_{i=1}^n (j\le \lfloor\frac{ai+b}{c}\rfloor)\)
下面，我們只需要變換\(j\le \lfloor\frac{ai+b}{c}\rfloor\)，使求和式容易計算
注意：我們對對該不等式只能進行等價變換，非等價變換會導致多算或漏算，切記
\(j\le \lfloor\frac{ai+b}{c}\rfloor\)
\(\Leftrightarrow cj\le ai+b\)
依據：當x,y為整數,z為正整數,\(x\le \lfloor \frac{y}{z} \rfloor\Leftrightarrow xz \le y\)
\(\Leftrightarrow ai \ge cj-b\)
\(\Leftrightarrow i \ge \frac{cj-b}{a}\)
\(\Leftrightarrow i \ge \lceil \frac{cj-b}{a}\rceil\)
依據：當x,y為整數,z為正整數,\(x\ge \frac{y}{z}\Leftrightarrow x\ge \lceil \frac{y}{z} \rceil\)
\(\Leftrightarrow i \ge \lceil \frac{cj-b}{a}\rceil\)
\(\Leftrightarrow i \ge \lfloor \frac{cj-b+a-1}{a}\rfloor\)
依據：當x為整數,y為正整數,\(\lceil \frac{x}{y}\rceil=\lfloor \frac{x+y-1}{y}\rfloor\)，注意：這個依據不是很顯然，讀者可以自己證明一下
所以，原求和式
\(=\sum_{j=1}^{\lfloor\frac{an+b}{c}\rfloor} \sum_{i=1}^n i \ge \lfloor \frac{cj-b+a-1}{a}\rfloor\)
對于\(i \ge \lfloor \frac{cj-b+a-1}{a}\rfloor\)
由\(a,b<c\) , \(1 \le j \le \lfloor\frac{an+b}{c}\rfloor\)可知
\(\lfloor \frac{cj-b+a-1}{a}\rfloor \ge \lfloor \frac{c-b+a-1}{a}\rfloor = 1+\lfloor \frac{c-b-1}{a}\rfloor \ge 1\)
\(\lfloor \frac{cj-b+a-1}{a}\rfloor \le \lfloor \frac{c\lfloor\frac{an+b}{c}\rfloor-b+a-1}{a}\rfloor \le \lfloor \frac{an+b-b+a-1}{a}\rfloor =n\)

即\(1 \le\lfloor \frac{cj-b+a-1}{a}\rfloor \le n\)

這一點非常重要，這確保當\(i\)可以取到\(\lfloor\frac{cj-b+a-1}{a}\rfloor \le i \le n\) 的所有取值，

讀者可以想一想如果\(\lfloor \frac{cj-b+a-1}{a}\rfloor\)可能大于\(n\)或小于\(1\)時會發生什么

所以，原求和式
\(=\sum_{j=1}^{\lfloor\frac{an+b}{c}\rfloor} n-\lfloor \frac{cj-b+a-1}{a}\rfloor+1\)

\(=\sum_{j=1}^{\lfloor\frac{an+b}{c}\rfloor} n-\lfloor \frac{cj-b-1}{a}\rfloor\)

\(=n(\lfloor\frac{an+b}{c}\rfloor)-\sum_{j=1}^{\lfloor\frac{an+b}{c}\rfloor}\lfloor \frac{cj-b-1}{a}\rfloor=n(\lfloor\frac{an+b}{c}\rfloor)-f(c,-b-1,a,\lfloor\frac{an+b}{c}\rfloor)\)

即有等式2：\(f(a,b,c,n)=(\lfloor\frac{an+b}{c}\rfloor)n-f(c,-b-1,a,\lfloor\frac{an+b}{c}\rfloor)\)

這里可以看出左右兩式\(c,a\)的位置互換了

因為\(c>a\)

所以可以用等式1處理\(f(c,-b-1,a,\lfloor\frac{an+b}{c}\rfloor)\)

得\(f(c,-b-1,a,\lfloor\frac{an+b}{c}\rfloor)=f(c\,mod\,a,(-b-1)\,mod\,a,a,\lfloor\frac{an+b}{c}\rfloor)+\frac{1}{2} \lfloor \frac{c}{a} \rfloor \lfloor\frac{an+b}{c}\rfloor^2+\frac{1}{2}(\lfloor \frac{c}{a} \rfloor+2\lfloor\frac{-b-1}{a}\rfloor)\lfloor\frac{an+b}{c}\rfloor\)

即等式3：

\(f(a,b,c,n)=(\lfloor\frac{an+b}{c}\rfloor)n\frac{1}{2} \lfloor \frac{c}{a} \rfloor \lfloor\frac{an+b}{c}\rfloor^2-\frac{1}{2}(\lfloor \frac{c}{a} \rfloor+2\lfloor\frac{-b-1}{a}\rfloor)\lfloor\frac{an+b}{c}\rfloor-f(c\,mod\,a,(-b-1)\,mod\,a,a,\lfloor\frac{an+b}{c}\rfloor)\)

這一等式中

左式中\(f\)函式第一個引數和第三個引數分別是\(a,c\),右式中\(f\)的第一個引數和第三個引數分別是\(c\,mod\,a,a\)

而歐幾里德演算法（輾轉相除法）在運算時也是利用

\(gcd(a,b)=gcd(b\,mod\,a,a)\)（假設第一個引數<第二個引數）

也就是說在遞回求解程序中，類歐幾里德演算法的第1,3引數變化和歐幾里得演算法的兩個引數變化的方式完全一致

最終的引數狀態就是第一個引數變成0，對應著“特殊情況”，標志著遞回的結束

故求解時間復雜度和歐幾里德演算法一致，是\(O(log(max(a,c)))\)的

最后總結一下遞推方程：

\(f(0,b,c,n)=\lfloor\frac{b}{c} \rfloor n\)

\(a \ge c\)時

\(f(a,b,c,n)=f(a\,mod\,c,b\,mod\,c,c,n)+\frac{1}{2} \lfloor \frac{a}{c} \rfloor n^2+\frac{1}{2}(\lfloor \frac{a}{c} \rfloor+2\lfloor\frac{b}{c}\rfloor)n\) （等式1）

\(a<c\)時

\(f(a,b,c,n)=(\lfloor\frac{an+b}{c}\rfloor)n-f(c,-b-1,a,\lfloor\frac{an+b}{c}\rfloor)\) （等式2）

由于等式3是由等式1和等式2通過代換推出的，遞回時沒有必要進行等式3的計算，

這道題要求

\(\sum_{i=1}^n [(pi)\,mod\,q]\)

即\(\sum_{i=1}^npi-q\sum_{i=1}^n \lfloor \frac{pi}{q} \rfloor\)

\(=p\frac{n(n+1)}{2}-qf(p,0,q,n)\)

遞回求解\(f(p,0,q,n)\)即可算出答案

特別提醒：在編程時，負數參與的整除和模運算會有一些問題

比如對于計算機來說，\(-5/3\)的值為\(-1\)，而\(\lfloor\frac{-5}{3}\rfloor=-2\)

模的情況也是比較復雜，讀者可以自行測驗

注：筆者在寫完這篇文章的一個月后做了一道相關題目，WA了10次，最后發現就是整除和模沒有處理好，說明這一點很容易錯，理論上整除和模都要自己手寫，除非能保證某個演算法絕對不會出現負數參與整除或模，

標籤：C++

上一篇：C++17結構化系結

下一篇：C++仿函式