题解 P5170 【【模板】类欧几里得算法】

· · 题解

Updated: 2020/03/02 修复了一些错误

提示

前置姿势：如何求gcd；和式的性质；基本数论姿势

阅读本文大约需要20 ~ 60min，由于我比较菜，因此部分东西可能讲的有些繁琐，还请各位神仙轻点吊打。

题目大意

求

\sum^n_{i=0}\left\lfloor\frac{ai+b}{c}\right\rfloor \bmod {998244353} \sum^n_{i=0}\left\lfloor\frac{ai+b}{c}\right\rfloor^2 \bmod {998244353} \sum^n_{i=0}i\left\lfloor\frac{ai+b}{c}\right\rfloor \bmod {998244353}

分析

我们设

\varphi(a,b,c,n)=\sum^n_{i=0}\left\lfloor\frac{ai+b}{c}\right\rfloor \psi (a,b,c,n)=\sum^n_{i=0}\left\lfloor\frac{ai+b}{c}\right\rfloor^2 \phi(a,b,c,n)=\sum^n_{i=0}i\left\lfloor\frac{ai+b}{c}\right\rfloor

这几个看起来最可做的应该就是\varphi了，所以我们不妨先看一下\varphi，发现它的主要部分就是后面这个东西，所以我们令\xi(i)=\left\lfloor\frac{ai+b}{c}\right\rfloor。我们思考，当i增大时，\xi会如何变化。

注意到\xi(i)=\left\lfloor\frac{ai}{c}+\frac{b}{c}\right\rfloor，而其又涉及到了取整，因此这个式子的变化率与a,b和c的大小关系有关。

思考当a\geqslant c时，\frac{a(i+1)}{c}=\frac{ai}{c}+\frac{a}{c}\geqslant \frac{ai}{c}+1，同理b\geqslant c时，显然有\frac bc\geqslant 1。

由取整函数显然的性质：

\left\lfloor r+N \right\rfloor=N+\left\lfloor r\right\rfloor, \forall r\in\mathbb R, N\in\mathbb N

暗示了我们要将取整函数内的式子拆开。

首先注意以下初等结论：

结论1.1：\left\lfloor\frac{Ax}{y}\right\rfloor=\left\lfloor\frac{A(x \bmod y)}{y}\right\rfloor+A\left\lfloor\frac xy\right\rfloor

立即可以看出

\xi(i)=\left\lfloor\frac{ai+b}{c}\right\rfloor=\left\lfloor\frac{(a\bmod c)i+(b\bmod c)}{c}\right\rfloor+i\left\lfloor\frac ac\right\rfloor + \left\lfloor\frac bc\right\rfloor

替换可得以下引理

引理1.1： \sum^n_{i=1}\left\lfloor\frac{ai+b}{c}\right\rfloor=\sum^n_{i=1}\left\lfloor\frac{(a\bmod c)i+(b\bmod c)}{c}\right\rfloor+i\left\lfloor\frac ac\right\rfloor + \left\lfloor\frac bc\right\rfloor

通过引理，立即就可以得到下面的显然结论。

定理1.1：\varphi(a,b,c,n)=\varphi(a\bmod c, b\bmod c, c, n)+\frac{n(n+1)}{2}\left\lfloor\frac{a}{c}\right\rfloor+(n+1)\frac{b}{c}

定理1.1是至关重要的，因为他可以把a,b\geqslant c的情形转化为a,b < c的情形。

通过这一步转化，我们发现\xi(i)的取值不可能超过n。这十分重要，因为接下来的变换会使得上限变为\xi(i)，这保证了我们可以将原问题转化为规模更小的子问题。

下面，重新考虑\varphi函数，不难看出：

\sum^n_{i=0}\left\lfloor\frac{ai+b}{c}\right\rfloor=\sum^n_{i=0}\sum^{\left\lfloor\frac{an+b}{c}\right\rfloor}_{d=1}\left[\left\lfloor\frac{ai+b}{c}\right\rfloor\geqslant d\right] =\sum^n_{i=0}\sum^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}_{d=0}\left[a^{-1}c\left\lfloor\frac{ai+b}{c}\right\rfloor\geqslant a^{-1}c(d+1)>a^{-1}(cd+c-1)\right] =\sum^n_{i=0}\sum^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}_{d=0}\left[i> \frac{cd+c-b-1}a\right]

考虑到右侧的谓词函数（艾弗森括号）实质上统计的是大于\frac{cd+c-b-1}a的i的个数，其恰好可以看做n-\left\lfloor\frac{cd+c-b-1}a\right\rfloor，因此，不难看出我们可以把原式化为与i无关的和式。

\varphi(a,b,c,n)=\sum^n_{i=0}\sum^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}_{d=0}\left[i> \frac{cd+c-b-1}a\right]=\sum^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}_{d=0}{n-\left\lfloor \frac{cd+c-b-1}a\right\rfloor} =n\left\lfloor\frac{an+b}{c}\right\rfloor-\sum^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}_{d=0}{\left\lfloor \frac{cd+c-b-1}a\right\rfloor}=n\left\lfloor\frac{an+b}{c}\right\rfloor-\varphi(c, c-b-1, a, \left\lfloor\frac{an+b}{c}\right\rfloor-1)

因此，可得到如下定理：

定理1.2：\varphi(a,b,c,n)=\varphi(c,c-b-1, a, \xi(n)-1)

这个定理是我们这道题的基础，是极其重要的。现在，通过定理1.2，我们可以在O(\log a)的时间复杂度内求出\varphi函数的值。（为什么？考虑你计算gcd的时间复杂度，会发现他们本质是一样的。）

当然，你要注意，我们在推导过程中除以了a，因此你需要特判a=0的情况。

40分-参考代码

ll f(ll a, ll b, ll c, ll n) {
    if (a == 0)
        return (n + 1) * (b / c) % mod;
    if (a >= c || b >= c)
        return (f(a % c, b % c, c, n) + (a / c) * n % mod * (n + 1) % mod * inv2 % mod + (b / c) * (n + 1) % mod) % mod;
    ll m = (a * n + b) / c;
    return (n * m % mod - f(c, c - b - 1, a, m - 1)) % mod;
}

下面，我们思考函数\psi，有了刚才的经验，显然我们还是需要先将a,b\geqslant c的情形转换。

由于推导过程非常简单与显然，各位神仙也不用听我在这里啰嗦，这里不再赘述，只给出一个大致的推导过程。

\psi(a,b,c,n)=\sum^n_{i=0}{\left(\left\lfloor\frac{(a\bmod c)i+b\bmod c}{c}\right\rfloor+i\left\lfloor\frac{a}{c}\right\rfloor+\left\lfloor\frac{b}{c}\right\rfloor\right)^2} =\sum^n_{i=0}\left(i^2\left\lfloor\frac{a}{c}\right\rfloor^2+\left\lfloor\frac{b}{c}\right\rfloor^2+\left\lfloor\frac{(a\bmod c)i+b\bmod c}{c}\right\rfloor^2+2i\left\lfloor\frac{a}{c}\right\rfloor\left\lfloor\frac{(a\bmod c)i+b\bmod c}{c}\right\rfloor+2\left\lfloor\frac{b}{c}\right\rfloor\left\lfloor\frac{(a\bmod c)i+b\bmod c}{c}\right\rfloor+2i\left\lfloor\frac{a}{c}\right\rfloor\left\lfloor\frac{b}{c}\right\rfloor\right)

观察这个乱七八糟的式子，有关i的项当中，i^2通过初等知识O(1)得到它的前缀和，而\left\lfloor\frac{(a\bmod c)i+b\bmod c}{c}\right\rfloor的前缀和显然就是\varphi(a\bmod c, b\bmod c, c,n)，\left\lfloor\frac{(a\bmod c)i+b\bmod c}{c}\right\rfloor^2的前缀和显然是\psi(a\bmod c, b\bmod c, c, n)，至于i\left\lfloor\frac{(a\bmod c)i+b\bmod c}{c}\right\rfloor，……我们发现这其实就是我们还没讨论到的函数\phi！也就是说，我们想求\psi，就必须要求出函数\phi。不妨先假设上帝给了你求这个函数的方法，我们先把注意力集中在\psi上。

原式可重新写作

\psi(a,b,c,n)=\psi(a\bmod c, b\bmod c, c, n)+\left\lfloor\frac{a}{c}\right\rfloor^2S_2(n)+(n+1)\left\lfloor\frac{b}{c}\right\rfloor^2+2\left\lfloor\frac{a}{c}\right\rfloor \phi(a\bmod c, b\bmod c, c, n)+2\left\lfloor\frac{b}{c}\right\rfloor \varphi(a \bmod c, b\bmod c, c, n) + 2\left\lfloor\frac{a}{c}\right\rfloor\left\lfloor\frac{b}{c}\right\rfloor S_1(n)

其中S_k(n)=\sum^n_{i=1}i^k，显然可以O(1)得出答案。

同样，运用这种手段，我们可以把a,b\geqslant c的情形转化到a,b<c。

随后，我们还是按照惯例试着去推一下式子。

\psi(a,b,c,n)=\sum^n_{i=0}\sum^{\left\lfloor \frac{an+b}{c} \right\rfloor^2}_{d=1}\left[\left\lfloor \frac{ai+b}{c} \right\rfloor^2\geqslant d \right]

注意到这个式子我们是不能使用的，原因是我们只保证了\xi(n)=\left\lfloor \frac{an+b}{c} \right\rfloor\leqslant n，在最坏情况下，\xi^2(n)会达到n^2级别，这将使得我们的状态永远无法转移。因此，我们需要使用一个小技巧。

结论1.2：X^2=-X+2\sum^X_{i=0} i

因此，应用这个结论，原式等价于

\psi(a,b,c,n)=\sum^n_{i=0} 2\sum^{\left\lfloor\frac{ai+b}{c}\right\rfloor}_{d=1}d-\left\lfloor\frac{ai+b}{c}\right\rfloor =2\sum^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}_{d=0}(d+1)\sum^n_{i=0}\left[\frac{(ai+b)}{c}\geqslant d+1\right]-\left\lfloor\frac{ai+b}{c}\right\rfloor =2\sum^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}_{d=0}(d+1)\sum^n_{i=0}\left[a^{-1}ai\geqslant a^{-1}(c(d+1)-b) > a^{-1}(cd+c-b-1)\right]-\varphi(a,b,c,n) =2\sum^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}_{d=0}(d+1)\sum^n_{i=0}\left[i > \frac{cd+c-b-1}{a}\right]-\varphi(a,b,c,n) =2\sum^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}_{d=0}(d+1)\left(n-\left\lfloor\frac{cd+c-b-1}{a}\right\rfloor\right)-\varphi(a,b,c,n) =2n\sum^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}_{d=0}(d+1)-2\sum^{\left\lfloor\frac{an+b}{c}\right\rfloor-1}_{d=0}(d+1)\left\lfloor\frac{cd+c-b-1}{a}\right\rfloor-\varphi(a,b,c,n) =n\xi(n)(\xi(n)+1)-2\phi(c,c-b-1,a,\xi(n)-1)-2\varphi(c,c-b-1, a, \xi(n)-1) - \varphi(a,b,c,n)

因此，我们可以通过这个式子，把计算规模为N的\psi，转化为了计算规模为n的\phi，且保证了n<N

定理1.3：\varphi(a,b,c,n) = n\xi(n)(\xi(n)+1)-2\phi(c,c-b-1,a,\xi(n)-1)-2\varphi(c,c-b-1, a, \xi(n)-1) - \varphi(a,b,c,n)

最后，我们只需要求出\phi即可，这非常的简单。按照惯例，我们依然将a,b\geqslant c的情形进行转换。由于前面已经进行了详细的解释，这里的思想是一样的，因此我们直接给出结论。

结论1.3：\phi(a,b,c,n)=\left\lfloor\frac{a}{c}\right\rfloor S_2(n)+\left\lfloor\frac bc\right\rfloor S_1(n) + \phi(a\bmod c, b\bmod c, c, n)

然后，根据a,b<c的情形，用\xi替代原式，这里同样十分简单，过程与上面相同，直接给出结论。

定理1.4：\phi(a,b,c,n)=\frac{2\xi(n)S_1(n)-\varphi(c, c - b - 1, a, \xi(n)-1) - \psi(c, c - b - 1, a, \xi(n) - 1)}{2}

同理，运用定理1.4，可将计算规模为N'的\phi，转化为计算规模为n'的\psi

由定理1.3，显然可得n'<N'=n<N'，实际上，有以下结论

定理 1.5：\lim \frac{N'}{n'}=0

证明也并不困难，只需要考虑进行一次变换后，子问题的规模至多为原文题的O(\log *)级别。这里不再给出。实际上，这说明了这个算法的时间复杂度为O(\log a)

实现

实际上，不难写出这个暴力程序：

ll f(ll a, ll b, ll c, ll n) {
    if (a == 0)
        return (n + 1) * (b / c) % mod;
    if (n == 0) return (b / c);
    if (a >= c || b >= c)
        return (f(a % c, b % c, c, n) + (a / c) * n % mod * (n + 1) % mod * inv2 % mod + (b / c) * (n + 1) % mod) % mod;
    ll m = (a * n + b) / c;
    return (n * m % mod - f(c, c - b - 1, a, m - 1)) % mod;
}

ll g(ll a, ll b, ll c, ll n) {
    if (a == 0) return (b / c) * n % mod * (n + 1) % mod * inv2 % mod;
    if (n == 0) return 0;
    if (a >= c || b >= c) return (g(a % c, b % c, c, n) + (a / c) * n % mod * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod + (b / c) * n % mod * (n + 1) % mod * inv2 % mod) % mod;
    ll m = (a * n + b) / c;
    return (n * (n + 1) % mod * m % mod - f(c, c - b - 1, a, m - 1) - h(c, c - b - 1, a, m - 1)) % mod * inv2 % mod;
}

ll h(ll a, ll b, ll c, ll n) {
    if (a == 0) return (n + 1) * (b / c) % mod * (b / c) % mod;
    if (n == 0) return (b / c) * (b / c) % mod;
    if (a >= c || b >= c)
        return ((a / c) * (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 + 
                h(a % c, b % c, c, n) + 2 * (a / c) % mod * g(a % c, b % c, c, n) % mod + 
                2 * (b / c) % mod * f(a % c, b % c, c, n) % mod) % mod;
    ll m = (a * n + b) / c;
    return (n * m % mod * (m + 1) % mod - 2 * g(c, c - b - 1, a, m - 1) - 2 * f(c, c - b - 1, a, m - 1) - f(a, b, c, n)) % mod;
}

但这是无法通过的。原因在于，每次计算时，许多东西被重复计算了。这不是很好处理，因为你无法记录哪些信息被计算过了。一个更好的解决方案是，将\varphi, \psi, \phi一起计算，并统一使用一个结构体来维护。

有了上面的讨论，您应该可以很快A掉这题

下面给出一个参考程序

#pragma GCC optimize(3)
#include <bits/stdc++.h>

typedef long long ll;
constexpr int mod = 998244353;

constexpr ll inv2 = 499122177;
constexpr ll inv6 = 166374059;

ll f(ll a, ll b, ll c, ll n);
ll g(ll a, ll b, ll c, ll n);
ll h(ll a, ll b, ll c, ll n);

struct Query {
    ll f, g, h;
};

Query solve(ll a, ll b, ll c, ll n) {
    Query ans, tmp;
    if (a == 0) {
        ans.f = (n + 1) * (b / c) % mod;
        ans.g = (b / c) * n % mod * (n + 1) % mod * inv2 % mod;
        ans.h = (n + 1) * (b / c) % mod * (b / c) % mod;
        return ans;
    }
    if (a >= c || b >= c) {
        tmp = solve(a % c, b % c, c, n);
        ans.f = (tmp.f + (a / c) * n % mod * (n + 1) % mod * inv2 % mod + (b / c) * (n + 1) % mod) % mod;
        ans.g = (tmp.g + (a / c) * n % mod * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod + (b / c) * n % mod * (n + 1) % mod * inv2 % mod) % mod;
        ans.h = ((a / c) * (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 + 
                tmp.h + 2 * (a / c) % mod * tmp.g % mod + 2 * (b / c) % mod * tmp.f % mod) % mod;
        return ans;
    }
    ll m = (a * n + b) / c;
    tmp = solve(c, c - b - 1, a, m - 1);
    ans.f = (n * (m % mod) % mod - tmp.f) % mod;
    ans.g = (n * (n + 1) % mod * (m % mod) % mod - tmp.f - tmp.h) % mod * inv2 % mod;
    ans.h = (n * (m % mod) % mod * ((m + 1) % mod) % mod - 2 * tmp.g - 2 * tmp.f - ans.f) % mod;
    return ans;
}

inline char nc() {
    static char buf[1000000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 1000000, stdin), p1 == p2) ? EOF : *p1++;
}

inline ll read() {
    ll res = 0;
    char ch;
    do ch = nc(); while (ch < 48 || ch > 57);
    do res = res * 10 + ch - 48, ch = nc(); while (ch >= 48 && ch <= 57);
    return res;
}

char pbuf[1 << 20], *pp = pbuf;
inline void push(const char &c) {
    if (pp - pbuf == 1 << 20) fwrite(pbuf, 1, 1 << 20, stdout), pp = pbuf;
    *pp++ = c;
}
inline void write(int x) {
    static ll sta[35];
    ll top = 0;
    do {
        sta[top++] = x % 10, x /= 10;
    } while (x);
    while (top) push(sta[--top] + '0');
}

int main() {
    ll t = read();
    ll n, a, b, c;
    while (t--) {
        n = read(), a = read(), b = read(), c = read();
        Query ans = solve(a, b, c, n);
        printf("%lld %lld %lld\n", (ans.f + mod) % mod, (ans.h + mod) % mod, (ans.g + mod) % mod);
    }
    return 0;
}