【算法學習筆記】類歐幾里得算法

一個基礎的數論問題。

試求 \(\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor\) 的值，其中：\(a,b \ge 0\)，\(n,c >0\)

在Atcoder的AC庫中有這樣一個函數能夠在 \(\mathcal{O}(lg(n + c + a + b))\) 的時間內解決問題。

函數代碼 ↓

using ll = long long;
ll floor_sum(ll n, ll m, ll a, ll b) {
    ll ans = 0;
    if (a >= m) {
        ans += (n - 1) * n * (a / m) / 2;
        a %= m;
    }
    if (b >= m) {
        ans += n * (b / m);
        b %= m;
    }

    ll y_max = (a * n + b) / m, x_max = (y_max * m - b);
    if (y_max == 0) return ans;
    ans += (n - (x_max + a - 1) / a) * y_max;
    ans += floor_sum(y_max, a, m, (a - x_max % a) % a);
    return ans;
}

好奇它的證實過程，而後在 OI wiki 上找到相應文檔，這個算法名爲：類歐幾里德算法

我的證實

另外補上我的證實：

當 \(a \ge c\) 或 \(b\ge c\) 時，

\[f(a,b,c,n) = \frac{n(n + 1)}{2}*\left\lfloor \frac{a}{c} \right\rfloor + (n + 1) * \left\lfloor \frac{b}{c} \right\rfloor + f(a\ mod\ c,b\ mod\ c,c,n). \]

當 \(a < c\) 而且 \(b < c\) 時，

設 \(m = \left\lfloor \frac{an+b}{c} \right\rfloor\)，

\[f(a,b,c,n) = mn - f(c,c - b - 1,a,m-1). \]

而後遞歸至 \(a = 0\) 便可.

例題

luoguP5171 Earthquake

數形結合, 把式子稍微簡單轉換一下, 套用類歐幾里得算法便可.

引入

設

\[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{split} 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{split} \]

那麼問題轉化爲了 \(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\) 的上界。

\[\begin{split} &=\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]\\ \end{split} \]

這樣作的目的是讓 \(j\) 擺脫 \(i\) 的限制，如今 \(i,j\) 都被 \(n\) 限制，而貢獻式看上去是一個條件，可是咱們仍把它叫做貢獻式，再對貢獻式作變換後就能夠改變 \(i,j\) 的限制關係。因而咱們作一些放縮的處理。首先把向下取整的符號拿掉

\[j<\left\lfloor \frac{ai+b}{c} \right\rfloor \Leftrightarrow j+1\leq \left\lfloor \frac{ai+b}{c} \right\rfloor \Leftrightarrow j+1\leq \frac{ai+b}{c}\\ \]

而後能夠作一些變換

\[j+1\leq \frac{ai+b}{c} \Leftrightarrow jc+c\le ai+b \Leftrightarrow jc+c-b-1< ai \]

最後一步，向下取整獲得：

\[jc+c-b-1< ai\Leftrightarrow \left\lfloor\frac{jc+c-b-1}{a}\right\rfloor< i \]

這一步的重要意義在於，咱們能夠把變量 \(i\) 消掉了！具體地，令 \(m=\left\lfloor \frac{an+b}{c} \right\rfloor\)，那麼原式化爲

\[\begin{split} 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{split} \]

這是一個遞歸的式子。而且你發現 \(a,c\) 分子分母換了位置，又能夠重複上述過程。先取模，再遞歸。這就是一個展轉相除的過程，這也是類歐幾里德算法的得名。

容易發現時間複雜度爲 \(\mathcal{O}(lg(n + m + a + b))\)。

同時關於 類歐幾里德算法 有兩個函數的拓展

擴展

理解了最基礎的類歐幾里德算法，咱們再來思考如下兩個變種求和式：

\[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{split} &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{split} \]

這時咱們設 \(t=\left\lfloor\frac{jc+c-b-1}{a}\right\rfloor\)，能夠獲得

\[\begin{split} &=\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{split} \]

推導 h

一樣的，首先取模：

\[\begin{split} 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{split} \]

考慮 \(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{split} &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{split} \]

接下來考慮化簡前一部分：

\[\begin{split} &\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{split} \]

所以

\[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)\) 的。

模板代碼實現

#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;
}