类欧几里得算法

问题引入

\[f(a, b, c, n) = \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor \]

其中 \(a, b, c, n\) 是常数，需要 \(\mathcal O(\log n)\) 的做法。

若 \(a \geq c\) 或 \(b \geq c\)，我们可以将 \(a, b\) 对 \(c\) 取模以简化问题。

考虑到

\[x = \left \lfloor \frac{x}{c}\right \rfloor c + x \bmod 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 \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 \\ &= \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\)。

\[S(i)={\left\lfloor\frac{ai + b}{c}\right\rfloor} - 1 \]

再进行转化

\[\begin{split} \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor &= \sum_{i=0}^n \sum_{j=0}^{S(i)} 1 \\ &= \sum_{j=0}^{S(n)} \sum_{i=0}^n \left[ j \leq S(i) \right] \end{split} \]

考虑到

\[\begin{split} j \leq S(i) &\iff j + 1 \leq {\left\lfloor\frac{ai + b}{c}\right\rfloor} \\ 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 \leq ai + b \\ jc + c \leq ai + b &\iff jc + c - b \leq ai \\ jc + c - b \leq ai &\iff jc + c - b - 1 < ai \\ jc + c - b - 1 < ai &\iff \left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor < i \end{split} \]

\[\begin {split} \sum_{j=0}^{S(n)} \sum_{i=0}^n \left[ j \leq S(i) \right] &= \sum_{j=0}^{S(n)} \sum_{i=0}^n \left[ i > \left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor\right] \\ &= \sum_{j=0}^{S(n)} \left(n - \left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor \right) \\ &= (S(n) + 1)n - \sum_{j=0}^{S(n)} \left \lfloor \frac{c j + (c - b - 1)}{a} \right \rfloor \\ &= (S(n) + 1)n - f(c, c - b - 1, a, S(n)) \end {split} \]

\[f(a, b, c, n) = (S(n) + 1)n - f(c, c - b - 1, a, S(n)) \]

可以发现，上述式子是一个递归式，我们不断重复上述过程，先取模，后递归，其实就是辗转相除的过程，时间复杂度 \(\mathcal O(\log n)\)。

拓展

我们再来推导两个变种求和式

\[g(a, b, c, n) = \sum_{i=0}^n i\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\)

引理 1.

\[\sum_{i=0}^n i^2 = \frac{n (n+1)(2n+1)}{6} \]

证明如下：

考虑到

\[(n+1)^3=n^3+3n^2+3n+1 \]

\[\begin {split} (n+1)^3-n^3 &= 3n^2+3n+1 \\ n^3-(n-1)^3 &=3(n-1)^2+3(n-1)+1 \\ &\cdots \\ 2^3 - 1^3&=3 \times (2-1)^2+3 \times(2-1) + 1 \\ \end {split} \]

将这 \(n\) 个等式左右两边相加，得到

\[(n+1)^3 - 1 = 3(1^2 + 2^2 + \cdots +n^2) + 3(1 + 2 + \cdots +n) + n \]

\[n^3+3n^2+3n = 3(1^2 + 2^2 + \cdots +n^2) + 3\frac{n(1 +n)}{2} + n \]

整理后得

\[\sum_{i=0}^n i^2 = \frac{n (n+1)(2n+1)}{6} \]

首先和 \(f\) 一样，对其取模（根据引理 1. 可将其展开）。

\[\begin{split} g(a, b, c, n) &= \sum_{i=0}^n i\left\lfloor\frac{ai + b}{c}\right\rfloor \\ &= \sum_{i=0}^n i\left\lfloor\frac{(\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 \\ &= \frac{n(n + 1)(2n+1)}{6} \left \lfloor \frac{a}{c}\right \rfloor + \frac{n(n + 1)}{2} \left \lfloor \frac{b}{c}\right \rfloor + g(a \bmod c , b\bmod c, c, n) \end{split} \]

其他部分推导与 \(f\) 类似。

\[S(i)={\left\lfloor\frac{ai + b}{c}\right\rfloor} - 1 \]

\[\begin {split} g(a, b, c, n) &= \sum_{j=0}^{S(n)} \sum_{i=0}^n i\left[ j \leq S(i) \right] \\ &= \sum_{j=0}^{S(n)} \sum_{i=0}^n i\left[ i > \left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor\right] \\ \end {split} \]

\[t(j) =\left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor \]

则有

\[\begin{split} g(a, b, c, n) &= \sum_{j=0}^{S(n)} \left((t(j) + 1) + (t(j) + 2) + \cdots + n\right) \\ &= \sum_{j=0}^{S(n)} \left(\frac{(t(j)+1+n)\times(n-t(j))}{2}\right) \\ &= \frac{1}{2}\left[ (S(n) + 1) n(n+1) - \sum_{j=0}^{S(n)}(t(j))^2 - \sum_{j=0}^{S(n)}t(j)\right] \\ &= \frac{1}{2} \left[(S(n) + 1) n(n+1) - h(c, c - b - 1, a, S(n)) - f(c, c - b - 1, a, S(n))\right] \end {split} \]

推导 \(h\)

同样套路，先取模

\[\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\lfloor\frac{(\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^2 \\ &= \sum_{i=0}^n \left(i \left \lfloor \frac{a}{c}\right \rfloor + \left \lfloor \frac{b}{c}\right \rfloor + \frac{(a \bmod c)i + (b\bmod c)}{c} \right)^2 \\ \end{split} \]

\[Q(i)=\frac{(a \bmod c)i + (b\bmod c)}{c} \]

拆开则有

\[\begin {split} h(a, b, c, n) &= \sum_{i=0}^n \left(i^2\left \lfloor \frac{a}{c}\right \rfloor^2 + \left \lfloor \frac{b}{c}\right \rfloor^2 + (Q(i))^2 + 2 i\left \lfloor \frac{a}{c}\right \rfloor \left \lfloor \frac{b}{c}\right \rfloor + 2i\left \lfloor \frac{a}{c}\right \rfloor Q(i) + 2\left \lfloor \frac{b}{c}\right \rfloor Q(i)\right) \\ &= \frac{n(n + 1)(2n + 1)}{6} \left \lfloor \frac{a}{c}\right \rfloor^2 + (n + 1)\left \lfloor \frac{b}{c}\right \rfloor^2 + h(a \bmod c, b \bmod c, c, n) + \\ & n(n+1)\left \lfloor \frac{a}{c}\right \rfloor \left \lfloor \frac{b}{c}\right \rfloor + 2\left \lfloor \frac{a}{c}\right \rfloor g(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) \end {split} \]

照样令

\[S(i)={\left\lfloor\frac{ai + b}{c}\right\rfloor} - 1 \\ t(j) =\left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor \]

发现平方不好处理，将平方转化为加法

\[\begin{split} n^2 &= 2\frac{n(n + 1)}{2} - n \\ &= \left(2\sum_{i=0}^n i \right) - n \end{split} \]

\[\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=0}^{S(i) + 1} j \right) - \left\lfloor\frac{ai + b}{c}\right\rfloor\right] \\ &= \left(2 \sum_{i=0}^n \sum_{j=0} ^{S(i) + 1} j\right) - f(a, b, c, n) \end{split} \]

现在只要解决前面的式子即可

\[\begin{split} \sum_{i=0}^n \sum_{j=0} ^{S(i) + 1} j &= \sum_{i=0}^n \sum_{j=0} ^{S(i)} (j+1) \\ &= \sum_{j=0}^{S(n)} (j+1)\sum_{i=0}^n \left[j < \left\lfloor\frac{ai + b}{c}\right\rfloor \right] \\ &= \sum_{j=0}^{S(n)} (j+1) (n - t(j)) \\ &= (S(n) + 1)n + \frac{S(n)(S(n) + 1)}{2} n - g(c, c - b - 1, a, S(n)) - f(c, c - b - 1, a, S(n)) \end{split} \]

综上所述

\[h(a, b, c, n) =(S(n)+1)(S(n) + 2) n - 2g(c, c - b - 1, a, S(n)) - 2f(c, c - b - 1, a, S(n)) - f(a, b, c, n) \]

当 \(a=0\) 时，上述三个式子都可以 \(\mathcal O(1)\) 计算。

代码实现时，因为三个函数交错递归，考虑三个一起整体递归，同步计算，时间复杂度为 \(\mathcal O(\log n)\)。

Luogu5170

#include 
const int mod = 998244353;
int InvTwo, InvSix;
inline int read()
{
    int cnt = 0; char ch = getchar(); bool op = 1;
    for (; ! isdigit(ch); ch = getchar())
        if (ch == '-') op = 0;
    for (; isdigit(ch); ch = getchar())
        cnt = cnt * 10 + ch - 48;
    return op ? cnt : - cnt;
}

inline int quick_pow(int a, int b)
{
    int ret = 1;
    for (; b; b >>= 1)
    {
        if (b & 1) ret = 1ll * ret * a % mod;
        a = 1ll * a * a % mod;
    }
    return ret % mod;
}

struct eu
{
    int f, g, h;
};

inline int S(int a, int b, int c, int i)
{
    return ((1ll * a * i + b) / c) - 1;
}

inline eu solve(int a, int b, int c, int n)
{
    if (a == 0)
    {
        eu now = {0, 0, 0};
        now.f = 1ll * (b / c) * (n + 1) % mod;
        now.h = 1ll * (b / c) * (b / c) % mod * (n + 1) % mod;
        now.g = 1ll * n * (n + 1) % mod * InvTwo % mod * (b / c) % mod;
        return now;
    }
    if (a >= c || b >= c)
    {
        eu now = {0, 0, 0};
        now.f = (now.f + 1ll * n * (n + 1) % mod * InvTwo % mod * (a / c) % mod) % mod;
        now.f = (now.f + 1ll * (n + 1) * (b / c) % mod) % mod;
        now.g = (now.g + 1ll * n * (n + 1) % mod * (2 * n + 1) % mod * InvSix % mod * (a / c) % mod) % mod;
        now.g = (now.g + 1ll * n * (n + 1) % mod * InvTwo % mod * (b / c) % mod) % mod;
        now.h = (now.h + 1ll * n * (n + 1) % mod * (2 * n + 1) % mod * InvSix % mod * (a / c) % mod * (a / c) % mod) % mod;
        now.h = (now.h + 1ll * (n + 1) * (b / c) % mod * (b / c) % mod) % mod;
        now.h = (now.h + 1ll * n * (n + 1) % mod * (a / c) % mod * (b / c) % mod) % mod;
        eu nxt = solve(a % c, b % c, c, n);
        now.f = (now.f + nxt.f) % mod;
        now.g = (now.g + nxt.g) % mod;
        now.h = (now.h + nxt.h) % mod;
        now.h = (now.h + 1ll * 2 * (a / c) % mod * nxt.g % mod) % mod;
        now.h = (now.h + 1ll * 2 * (b / c) % mod * nxt.f % mod) % mod;
        return now;
    }
    else
    {
        eu now = {0, 0, 0};
        eu nxt = solve(c, c - b - 1, a, S(a, b, c, n));
        now.f = 1ll * (S(a, b, c, n) + 1) * n % mod;
        now.f = (now.f - nxt.f) % mod;
        now.f = (now.f + mod) % mod;
        now.g = 1ll * InvTwo * ((((1ll * (S(a, b, c, n) + 1) * n % mod * (n + 1) % mod - nxt.h) % mod) - nxt.f) % mod) % mod;
        now.g = (now.g + mod) % mod;
        now.h = (1ll * (S(a, b, c, n) + 1) * (S(a, b, c, n) + 2) % mod * n % mod) % mod;
        now.h = (now.h - 1ll * 2 * nxt.g) % mod;
        now.h = (now.h - 1ll * 2 * nxt.f) % mod;
        now.h = (now.h - now.f) % mod;
        now.h = (now.h + mod) % mod;
        return now;
    }
}

int main()
{
    int t = read();
    InvTwo = quick_pow(2, mod - 2);
    InvSix = quick_pow(6, mod - 2);
    while (t --)
    {
        int n, a, b, c;
        n = read(), a = read(), b = read(), c = read();
        eu ans = solve(a, b, c, n);
        printf("%d %d %d\n", ans.f, ans.h, ans.g);
    }
    return 0;
}

Reference

oi-wiki