Luogu 5170 【模板】类欧几里得算法

原理不难但是写起来非常复杂的东西。

我觉得讲得非常好懂的博客。   传送门

我们设

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

$$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$$

先考虑一个结论

$$\left \lfloor \frac{Ax}{y} \right \rfloor = \left \lfloor \frac{A(x \mod y)}{y} \right \rfloor + A\left \lfloor \frac{x}{y} \right \rfloor$$

那么有

$$\left \lfloor \frac{ai + b}{c} \right \rfloor = \left \lfloor \frac{(a \mod c)i + (b \mod c)}{c} \right \rfloor + i\left \lfloor \frac{a}{c} \right \rfloor + \left \lfloor \frac{b}{c} \right \rfloor$$

这个东西可以把$a \geq c$ 或者$b \geq c$的情况转化成$a,b < c$的情况。

F

先看看这个比较好做的$f$。

注意到当$a \geq c$或者$b \geq c$的时候,我们可以用以上的结论把$a$或者$b$降下来,有

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

那么当$a, b < c$的时候,有

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

为了方便把$\left \lfloor \frac{an + b}{c} \right \rfloor$记为$m$。

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

$$= \sum_{i = 0}^{n}\sum_{j = 0}^{m - 1}[ai \geq cj + c - b]$$

$$= \sum_{i = 0}^{n}\sum_{j = 0}^{m - 1}[ai > cj + c - b - 1]$$

$$= \sum_{i = 0}^{n}\sum_{j = 0}^{m - 1}[i > \left \lfloor \frac{cj + c - b - 1}{a} \right \rfloor]$$

注意到这时候$i$这一项可以直接算了。

$$f(a, b, c, n) = \sum_{j = 0}^{m - 1}(n -  \left \lfloor \frac{cj + c - b - 1}{a} \right \rfloor)$$

$$= nm - f(c, c - b - 1, a, m - 1)$$

比较方便。

G

按照套路先做$a \geq c$或者$b \geq c$的情况。

$$g(a, b, c, n) = \left \lfloor \frac{b}{c} \right \rfloor\frac{n(n + 1)}{2} + \left \lfloor \frac{a}{c} \right \rfloor \frac{n(n + 1)(2n + 1)}{6}  + g(a \mod c, b \mod c, c, n)$$

然后接着按照上述套路弄$a, b < c$的情况。

$$g(a, b, c, n) = \sum_{j = 0}^{m - 1}\frac{(n - \left \lfloor \frac{cj + c - b - 1}{a} \right \rfloor)(n + 1 + \left \lfloor \frac{cj + c - b - 1}{a} \right \rfloor)}{2}$$

$$= \frac{1}{2}(mn(n + 1) - f(c, c - b - 1, a, m - 1) - h(c, c - b - 1, a, m - 1))$$

还要解决$h$。

H

当$a \geq c$或者$b \geq c$的时候

$$h(a, b, c, n) = (n + 1)\left \lfloor \frac{b}{c} \right \rfloor^2 + \frac{n(n + 1)(2n + 1)}{6}\left \lfloor \frac{a}{c} \right \rfloor^2 + n(n + 1)\left \lfloor \frac{b}{c} \right \rfloor\left \lfloor \frac{a}{c} \right \rfloor + h(a \mod c, b \mod c, c, n) + 2\left \lfloor \frac{a}{c} \right \rfloor g(a \mod c, b \mod c, c, n) + 2\left \lfloor \frac{b}{c} \right \rfloor f(a \mod c, b \mod c, c, n)$$

好麻烦啊……

当$a, b < c$的时候需要一些操作

$$n^2 = 2 \times \frac{n(n + 1)}{2} - n = 2\sum_{i = 1}^{n}i - n$$

可以得到

$$h(a, b, c, n) = \sum_{i = 0}^{n}(2\sum_{j = 1}^{\left \lfloor \frac{ai + b}{c} \right \rfloor}j - \left \lfloor \frac{ai + b}{c} \right \rfloor)$$

$$= 2\sum_{i = 0}^{n}\sum_{j = 1}^{\left \lfloor \frac{ai + b}{c} \right \rfloor}j - f(a, b, c, n)$$

$$= 2\sum_{j = 0}^{m - 1}(j + 1)\sum_{i = 0}^{n}[\left \lfloor \frac{ai + b}{c} \right \rfloor \geq j + 1] - f(a, b, c, n)$$

$$= 2\sum_{j = 0}^{m - 1}(j + 1)(n - \left \lfloor \frac{cj + c - b - 1}{a} \right \rfloor) - f(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)$$

大概就是这样了。

 

在计算的时候如果只考虑$(a, c)$这两项的话相当于每一次把$(a, c)$变成了$(c, a \% c)$,所以时间复杂度和欧几里得算法相同,递归层数是$log$层。

时间复杂度$O(Tlogn)$。

因为递归的式子很一致,在计算的时候应当三个东西一起算比较快。

Code:

#include <cstdio>
#include <cstring>
using namespace std;
typedef long long ll;

const ll P = 998244353LL;

template <typename T>
inline void read(T &X) {
    X = 0; char ch = 0; T op = 1;
    for (; ch > '9'|| ch < '0'; ch = getchar())
        if (ch == '-') op = -1;
    for (; ch >= '0' && ch <= '9'; ch = getchar())
        X = (X << 3) + (X << 1) + ch - 48;
    X *= op;
}

template <typename T>
inline void inc(T &x, T y) {
    x += y;
    if (x >= P) x -= P;
}

template <typename T>
inline void sub(T &x, T y) {
    x -= y;
    if (x < 0) x += P;
}

inline ll fpow(ll x, ll y) {
    ll res = 1;
    for (; y > 0; y >>= 1) {
        if (y & 1) res = res * x % P;
        x = x * x % P;
    }
    return res;
}

namespace Likegcd {
    struct Node {
        ll f, g, h;    
    };
    
    #define f(now) now.f
    #define g(now) now.g
    #define h(now) now.h
    
    const ll inv2 = fpow(2, P - 2);
    const ll inv6 = fpow(6, P - 2);
    
    inline Node solve(ll a, ll b, ll c, ll n) {
        Node res;
        if (!a) {
            f(res) = (b / c) * (n + 1) % P;
            g(res) = (b / c) * (n + 1) % P * n % P * inv2 % P;
            h(res) = (b / c) * (b / c) % P * (n + 1) % P;
            return res;
        }
        
        f(res) = g(res) = h(res) = 0;
        if (a >= c || b >= c) {
            Node tmp = solve(a % c, b % c, c, n);
            inc(f(res), (a / c) * n % P * (n + 1) % P * inv2 % P);
            inc(f(res), (b / c) * (n + 1) % P);
            inc(f(res), f(tmp));
            
            inc(g(res), (a / c) * n % P * (n + 1) % P * ((2 * n + 1) % P) % P * inv6 % P);
            inc(g(res), (b / c) * n % P * (n + 1) % P * inv2 % P);
            inc(g(res), g(tmp));
            
            inc(h(res), (a / c) * (a / c) % P * n % P * (n + 1) % P * ((2 * n + 1) % P) % P * inv6 % P);
            inc(h(res), (b / c) * (b / c) % P * (n + 1) % P);
            inc(h(res), (a / c) * (b / c) % P * n % P * (n + 1) % P);
            inc(h(res), h(tmp));
            inc(h(res), 2LL * (a / c) % P * g(tmp) % P);
            inc(h(res), 2LL * (b / c) % P * f(tmp) % P);
            
            return res;
        }
        
        if (a < c && b < c) {
            ll m = (a * n + b) / c;
            Node tmp = solve(c, c - b - 1, a, m - 1);
            
            f(res) = n * m % P;
            sub(f(res), f(tmp));
            
            g(res) = n * (n + 1) % P * m % P;
            sub(g(res), f(tmp));
            sub(g(res), h(tmp));
            g(res) = g(res) * inv2 % P;
            
            h(res) = n * m % P * (m + 1) % P;
            sub(h(res), 2LL * g(tmp) % P);
            sub(h(res), 2LL * f(tmp) % P);
            sub(h(res), f(res));
            
            return res;
        }
        
        return res;
    }
        
}

int main() {
    int testCase;
    read(testCase);
    for (ll a, b, c, n; testCase--; ) {
        read(n), read(a), read(b), read(c);
        Likegcd :: Node ans = Likegcd :: solve(a, b, c, n);
        printf("%lld %lld %lld\n", ans.f, ans.h, ans.g);
    }
    return 0;
}
View Code

 

转载于:https://www.cnblogs.com/CzxingcHen/p/10365060.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值