跳转至

类欧几里得算法

推导过程
\[ \begin{align} f(a,b,c,n)&=\sum_{i=0}^n\lfloor\dfrac{a*i+b}{c}\rfloor\\ &=\sum_{i=0}^n\left\lfloor\frac{\left(\lfloor\frac{a}{c}\rfloor*c+a\ mod\ c\right)*i+\left(\lfloor\frac{b}{c}\rfloor*c+b\ mod\ c\right)}{c}\right\rfloor\\ &=\sum_{i=0}^n\left\lfloor\lfloor\frac{a}{c}\rfloor*i+\frac{\left(a\ mod\ c\right)*i}{c}+\lfloor\frac{b}{c}\rfloor+\frac{\left(b\ mod\ c\right)}{c}\right\rfloor\\ &=\frac{n*(n+1)}{2}*\lfloor\frac{a}{c}\rfloor+(n+1)*\lfloor{\frac{b}{c}}\rfloor+\sum_{i=0}^n\left\lfloor\frac{(a\ mod\ c)*i+(b\ mod\ c)}{c}\right\rfloor\\ &=\frac{n*(n+1)}{2}*\lfloor\frac{a}{c}\rfloor+(n+1)*\lfloor{\frac{b}{c}}\rfloor+f(a\ mod\ c,b\ mod\ c,c,n) \end{align} \]
\[ \begin{align} f(a,b,c,n)&=\sum_{i=0}^n\lfloor\dfrac{a*i+b}{c}\rfloor\\ &=\sum_{i=0}^n\sum_{j=0}^{\lfloor\frac{a*i+b}{c}\rfloor-1}1\\ &=\sum_{j=0}^{\lfloor\frac{a*n+b}{c}\rfloor-1}\sum_{i=0}^n\left[j<\lfloor\frac{a*i+b}{c}\rfloor\right] \end{align} \]
\[ \begin{align} j<\lfloor\frac{a*i+b}{c}\rfloor\Longleftrightarrow{j}+1\leq\frac{a*i+b}{c}\Longleftrightarrow{j}*c+c\le{a}*i+b\Longleftrightarrow\left\lfloor\frac{j*c+c-b-1}{a}\right\rfloor<i \end{align} \]
\[ \begin{align} f(a,b,c,b)&=\sum_{j=0}^{\lfloor\frac{a*n+b}{c}\rfloor-1}\sum_{i=0}^n\left[\left\lfloor\frac{j*c+c-b-1}{a}\right\rfloor<i\right]\\ &=\sum_{j=0}^{\lfloor\frac{a*n+b}{c}\rfloor-1}\left(n-\left\lfloor\frac{j*c+c-b-1}{a}\right\rfloor\right)\\ &=n*\lfloor\frac{a*n+b}{c}\rfloor-f(c,c-b-1,a,\lfloor\frac{a*n+b}{c}\rfloor-1) \end{align} \]
\[ \begin{align} g(a,b,c,n)&=\sum_{i=0}^ni*\lfloor\frac{a*i+b}{c}\rfloor\\ &=\frac{n*(n+1)*(2n+1)}{6}*\lfloor\frac{a}{c}\rfloor+\frac{n*(n+1)}{2}*\lfloor\frac{b}{c}\rfloor+g(a\ mod\ c,b\ mod\ c,c,n) \end{align} \]
\[ \begin{align} g(a,b,c,n)&=\sum_{i=0}^ni*\lfloor\dfrac{a*i+b}{c}\rfloor\\ &=\sum_{i=0}^n\sum_{j=0}^{\lfloor\frac{a*i+b}{c}\rfloor-1}i\\ &=\sum_{j=0}^{\lfloor\frac{a*n+b}{c}\rfloor-1}\sum_{i=0}^ni*\left[j<\lfloor\frac{a*i+b}{c}\rfloor\right]\\ &=\sum_{j=0}^{\lfloor\frac{a*n+b}{c}\rfloor-1}\sum_{i=0}^ni*\left[\lfloor\frac{j*c+c-b-1}{a}\rfloor<i\right]\\ &=\sum_{j=0}^{\lfloor\frac{a*n+b}{c}\rfloor-1}\frac{(\lfloor\frac{j*c+c-b-1}{a}\rfloor+1+n)*(n-\lfloor\frac{j*c+c-b-1}{a}\rfloor)}{2}\\ &=\sum_{j=0}^{\lfloor\frac{a*n+b}{c}\rfloor-1}\left(\frac{n^2+n-\lfloor\frac{j*c+c-b-1}{a}\rfloor-(\lfloor\frac{j*c+c-b-1}{a}\rfloor)^2}{2}\right)\\ &=\frac{1}{2}*n*(n+1)*\lfloor\frac{a*n+b}{c}\rfloor-\frac{1}{2}*h(c,c-b-1,a,\lfloor\frac{a*n+b}{c}\rfloor-1)\\ &-\frac{1}{2}*f(c,c-b-1,a,\lfloor\frac{a*n+b}{c}\rfloor-1) \end{align} \]
\[ \begin{align} h(a,b,c,n)&=\sum_{i=0}^n\left\lfloor\frac{a*i+b}{c}\right\rfloor^2\\ &=(n+1)*\lfloor\frac{b}{c}\rfloor^2+n*(n+1)*\lfloor\frac{a}{c}\rfloor*\lfloor\frac{b}{c}\rfloor+\frac{n*(n+1)*(2n+1)}{6}*\lfloor\frac{a}{c}\rfloor^2\\ &+2*\lfloor\frac{b}{c}\rfloor*f(a\ mod\ c,b\ mod\ c,c,n) +2*\lfloor\frac{a}{c}\rfloor*g(a\ mod\ c,b\ mod\ c,c,n) \end{align} \]
\[ \begin{align} n^2&=2*\frac{n*(n+1)}{2}-n=2*(\sum_{i=0}^ni)-n\\ h(a,b,c,n)&=\sum_{i=0}^n\left\lfloor\frac{a*i+b}{c}\right\rfloor^2\\ &=\sum_{i=0}^n\left[2*(\sum_{j=1}^{\lfloor\frac{a*i+b}{c}\rfloor}j)-\lfloor\frac{a*i+b}{c}\rfloor\right]\\ &=2*(\sum_{i=0}^n\sum_{j=1}^{\lfloor\frac{a*i+b}{c}\rfloor}j)-f(a,b,c,n)\\ \sum_{i=0}^n\sum_{j=1}^{\lfloor\frac{a*i+b}{c}\rfloor}j&=\sum_{i=0}^n\sum_{j=0}^{\lfloor\frac{a*n+b}{c}\rfloor-1}j+1\\ &=\sum_{j=0}^{\lfloor\frac{a*n+b}{c}\rfloor-1}(j+1)\sum_{i=0}^n\left[j<\lfloor\frac{a*i+b}{c}\rfloor\right]\\ &=\sum_{j=0}^{\lfloor\frac{a*n+b}{c}\rfloor-1}(j+1)\sum_{i=0}^n[\lfloor\frac{j*c+c-b-1}{a}\rfloor<i]\\ &=\sum_{j=0}^{\lfloor\frac{a*n+b}{c}\rfloor-1}(j+1)*(n-\lfloor\frac{j*c+c-b-1}{a}\rfloor)\\ &=\frac{1}{2}*n*\lfloor\frac{a*n+b}{c}\rfloor*(\lfloor\frac{a*n+b}{c}\rfloor+1)-\sum_{j=0}^{\lfloor\frac{a*n+b}{c}\rfloor-1}(j+1)*\left\lfloor\frac{j*c+c-b-1}{a}\right\rfloor\\ &=\frac{1}{2}*n*\lfloor\frac{a*n+b}{c}\rfloor*(\lfloor\frac{a*n+b}{c}\rfloor+1)-g(c,c-b-1,a,\lfloor\frac{a*n+b}{c}\rfloor-1)\\ &-f(c,c-b-1,a,\lfloor\frac{a*n+b}{c}\rfloor-1)\\ h(a,b,c,n)&=n*\lfloor\frac{a*n+b}{c}\rfloor*(\lfloor\frac{a*n+b}{c}\rfloor+1)-2*g(c,c-b-1,a,\lfloor\frac{a*n+b}{c}\rfloor-1)\\ &-2*f(c,c-b-1,a,\lfloor\frac{a*n+b}{c}\rfloor-1)-f(a,b,c,n)\\ \end{align} \]
\[ f(a,b,c,n)= \begin{cases} \frac{n*(n+1)}{2}*\lfloor\frac{a}{c}\rfloor+(n+1)*\lfloor{\frac{b}{c}}\rfloor+f(a\ mod\ c,b\ mod\ c,c,n)&a\ge{c}\text{ 或 }b\ge{c}\\ n*\lfloor\frac{a*n+b}{c}\rfloor-f(c,c-b-1,a,\lfloor\frac{a*n+b}{c}\rfloor-1)&a<c\text{ 且 }b<c\\ \end{cases} \]
\[ g(a,b,c,n)= \begin{cases} \frac{n*(n+1)*(2n+1)}{6}*\lfloor\frac{a}{c}\rfloor+\frac{n*(n+1)}{2}*\lfloor\frac{b}{c}\rfloor+g(a\ mod\ c,b\ mod\ c,c,n)&a\ge{c}\text{ 或 }b\ge{c}\\ \frac{1}{2}*n*(n+1)*\lfloor\frac{a*n+b}{c}\rfloor-\frac{1}{2}*h(c,c-b-1,a,\lfloor\frac{a*n+b}{c}\rfloor-1)\\ \ -\frac{1}{2}*f(c,c-b-1,a,\lfloor\frac{a*n+b}{c}\rfloor-1)&a<c\text{ 且 }b<c\\ \end{cases} \]
\[ h(a,b,c,n)= \begin{cases} (n+1)*\lfloor\frac{b}{c}\rfloor^2+n*(n+1)*\lfloor\frac{a}{c}\rfloor*\lfloor\frac{b}{c}\rfloor+\frac{n*(n+1)*(2n+1)}{6}*\lfloor\frac{a}{c}\rfloor^2\\ +2*\lfloor\frac{b}{c}\rfloor*f(a\ mod\ c,b\ mod\ c,c,n)+2*\lfloor\frac{a}{c}\rfloor*g(a\ mod\ c,b\ mod\ c,c,n)&a\ge{c}\text{ 或 }b\ge{c}\\ n*\lfloor\frac{a*n+b}{c}\rfloor*(\lfloor\frac{a*n+b}{c}\rfloor+1)-2*g(c,c-b-1,a,\lfloor\frac{a*n+b}{c}\rfloor-1)\\ -2*f(c,c-b-1,a,\lfloor\frac{a*n+b}{c}\rfloor-1)-f(a,b,c,n)&a<c\text{ 且 }b<c\\ \end{cases} \]
模版
代码
C++
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
using i64 = long long;
const i64 P = 998244353;
i64 qpow(i64 a, i64 b = P - 2, i64 res = 1) {
    while (b) {
        if (b & 1) {
            res = res * a % P;
        }
        a = a * a % P;
        b >>= 1;
    }
    return res;
}
const i64 i2 = qpow(2);
const i64 i6 = qpow(6);
i64 calc(i64 n, i64 a, i64 b, i64 c) {
    i64 ac = a / c, bc = b / c;
    i64 m = (a * n + b) / c;
    i64 n1 = (n + 1) % P;
    i64 n_p = n % P, m_p = m % P;
    if (a == 0) {
        return bc * n1 % P;
    }
    if (a >= c || b >= c) {
        i64 res = calc(n, a % c, b % c, c);
        return (n_p * n1 % P * i2 % P * ac % P + bc * n1 % P + res) % P;
    }
    i64 res = calc(m - 1, c, c - b - 1, a);
    return (n_p * m_p % P - res + P) % P;
}
代码
C++
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
using i64 = long long;
const i64 P = 998244353;
i64 qpow(i64 a, i64 b = P - 2, i64 res = 1) {
    while (b) {
        if (b & 1) {
            res = res * a % P;
        }
        a = a * a % P;
        b >>= 1;
    }
    return res;
}
const i64 i2 = qpow(2);
const i64 i6 = qpow(6);
struct Node {
    i64 f, g, h;
    Node() = default;
    Node(i64 f, i64 g, i64 h) : f(f), g(g), h(h) {}
};
Node calc(i64 n, i64 a, i64 b, i64 c) {
    i64 m = (a * n + b) / c;
    i64 n1 = (n + 1) % P, n2 = (n * 2 + 1) % P;
    i64 ac = a / c % P, bc = b / c % P;
    i64 n_p = n % P, m_p = m % P;
    if (a == 0) {
        return Node(
            bc * n1 % P,
            bc * n_p % P * n1 % P * i2 % P,
            bc * bc % P * n1 % P
        );
    }
    if (a >= c || b >= c) {
        Node res = calc(n, a % c, b % c, c);
        return Node(
            (n_p * n1 % P * i2 % P * ac % P + bc * n1 % P + res.f) % P,
            (ac * n_p % P * n1 % P * n2 % P * i6 % P + bc * n_p % P * n1 % P * i2 % P + res.g) % P,
            (ac * ac % P * n_p % P * n1 % P * n2 % P * i6 % P + bc * bc % P * n1 % P + ac * bc % P * n_p % P * n1 % P + res.h + 2LL * bc * res.f % P + 2LL * ac * res.g % P) % P
        );
    }
    Node res = calc(m - 1, c, c - b - 1, a);
    long long f = (n_p * m_p % P - res.f + P) % P;
    return Node(
        f,
        ((m_p * n_p % P * n1 % P - res.h - res.f) * i2 % P + P) % P,
        ((n_p * m_p % P * (m_p + 1) % P - 2LL * res.g - 2LL * res.f - f) % P + P) % P
    );
}
例题

P5170 【模板】类欧几里得算法 - 洛谷

万能欧几里得

\(\sum_{i=0}^n\frac{a*i+b}{c}\)

模版
代码
C++
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
using i64 = long long;
const i64 P = 998244353;
struct Node {
    i64 cntU, cntR;
    i64 sum, f, g, h;
    Node() : cntU(0), cntR(0) {
        sum = f = g = h = 0;
    }
    friend Node operator+(Node a, Node b) {
        Node res;
        res.cntU = (a.cntU + b.cntU) % P;
        res.cntR = (a.cntR + b.cntR) % P;
        res.sum = (a.sum + b.sum + a.cntR * b.cntR % P) % P;
        res.f = (a.f + b.f + a.cntU * b.cntR % P) % P;
        res.g = (a.g + b.g + a.cntU * a.cntU % P * b.cntR % P + 2 * a.cntU * b.f % P) % P;
        res.h = (a.h + b.h + a.cntU * a.cntR % P * b.cntR % P + a.cntU * b.sum % P + a.cntR * b.f % P) % P;
        return res;
    }
};
Node qpow(Node a, i64 b, Node res = {}) {
    while (b) {
        if (b & 1) res = res + a;
        a = a + a;
        b >>= 1;
    }
    return res;
}
i64 S(i64 P, i64 L, i64 R, i64 Q) {
    return (P * L + R) / Q;
}
Node get(i64 P, i64 Q, i64 R, i64 L, Node a, Node b) {
    if (L == 0LL) return Node();
    R %= Q;
    if (P >= Q) {
        return get(P % Q, Q, R, L, a, qpow(a, P / Q) + b);
    }
    i64 m = S(P, L, R, Q);
    if (m == 0LL) return qpow(b, L);
    i64 cnt = L - S(Q, m, -R - 1, P);
    return qpow(b, (Q - R - 1) / P) + a + get(Q, P, (Q - R - 1) % P, m - 1, b, a) + qpow(b, cnt);
}
Node solve(i64 a, i64 c, i64 b, i64 n, Node U, Node R) {
    Node res = qpow(U, b / c) + get(a, c, b, n, U, R);
    res.f = (res.f + b / c) % P;
    res.g = (res.g + (b / c) % P * (b / c) % P) % P;
    return res;
}
例题

#6440. 万能欧几里得 - LibreOJ

\[ \sum_{x=1}^LA^x*B^{\lfloor\frac{P*X+R}{Q}\rfloor} \]
代码
C++
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
using i64 = long long;
const i64 P = 998244353;
using Mat = array<array<i64, 20>, 20>;
Mat operator+(Mat x, Mat y) {
    Mat res{};
    for (int i = 0 ; i < 20 ; i++) {
        for (int j = 0 ; j < 20 ; j++) {
            res[i][j] = (x[i][j] + y[i][j]) % P;
        }
    }
    return res;
}
Mat operator*(Mat x, Mat y) {
    Mat res{};
    for (int i = 0 ; i < 20 ; i++) {
        for (int j = 0 ; j < 20 ; j++) {
            for (int k = 0 ; k < 20 ; k++) {
                res[i][j] = (res[i][j] + x[i][k] * y[k][j] % P) % P;
            }
        }
    }
    return res;
}
struct Node {
    Mat cntU, cntR, sum;
    Node() {
        cntU = cntR = sum = array<array<i64, 20>, 20>{};
        for (int i = 0 ; i < 20 ; i++) {
            cntU[i][i] = cntR[i][i] = 1;
        }
    }
    Node(Mat cntU, Mat cntR, Mat sum) : cntU(cntU), cntR(cntR), sum(sum) {}
    friend Node operator+(Node a, Node b) {
        // a 是 x 为 [1, a.cntR], (P * x + R) / Q 为 a.cntU
        // b 是 x 为 [1, b.cntR], (P * x + R) / Q 为 b.cntU
        Node res;
        res.cntU = a.cntU * b.cntU;
        res.cntR = a.cntR * b.cntR;
        // sum(A^{x+a.cntR}*B^{s+a.cntU})
        res.sum = a.sum + a.cntR * b.sum * a.cntU;
        return res;
    }
};
Node qpow(Node a, i64 b, Node res = {}) {
    while (b) {
        if (b & 1) res = res + a;
        a = a + a;
        b >>= 1;
    }
    return res;
}
i64 S(i64 P, i64 L, i64 R, i64 Q) {
    return ((__int128) P * L + R) / Q;
}
Node get(i64 P, i64 Q, i64 R, i64 L, Node a, Node b) {
    if (L == 0LL) return Node();
    if (P >= Q) {
        return get(P % Q, Q, R, L, a, qpow(a, P / Q) + b);
    }
    i64 m = S(P, L, R, Q);
    if (m == 0LL) return qpow(b, L);
    i64 cnt = L - S(Q, m, -R - 1, P);
    return qpow(b, (Q - R - 1) / P) + a + get(Q, P, (Q - R - 1) % P, m - 1, b, a) + qpow(b, cnt);
}
Node solve(i64 a, i64 c, i64 b, i64 n, Node U, Node R) {
    // 合并 [0, 0] 和 [1, n]
    Node res = qpow(U, b / c) + get(a, c, b % c, n, U, R);
    return res;
}

int main() {
    i64 P, Q, R, L;
    int n;
    cin >> P >> Q >> R >> L >> n;
    Mat A{}, B{}, E{};
    for (int i = 0 ; i < n ; i++) {
        for (int j = 0 ; j < n ; j++) {
            cin >> A[i][j];
        }
    }
    for (int i = 0 ; i < n ; i++) {
        for (int j = 0 ; j < n ; j++) {
            cin >> B[i][j];
        }
    }
    for (int i = 0 ; i < 20 ; i++) {
        E[i][i] = 1;
    }
    Node u{B, E, Mat{}}, r{E, A, A};
    auto res = solve(P, Q, R, L, u, r).sum;
    for (int i = 0 ; i < n ; i++) {
        for (int j = 0 ; j < n ; j++) {
            cout << res[i][j] << " ";
        }
        cout << "\n";
    }
    return 0;
}

Ex - Popcount Sum

\[ \begin{align} \sum_{i=1}^{n}[i\%m==r]*popcount(i)\\ \sum_{i=0}^{\lfloor\frac{n-m}{r}\rfloor}popcount(i*m+r)\\ \sum_{i=0}^{\lfloor\frac{n-m}{r}\rfloor}\sum_{j=0}^{30}[\lfloor\dfrac{i*m+r}{2^j}\rfloor\equiv1\pmod2]\\ \sum_{i=0}^{\lfloor\frac{n-m}{r}\rfloor}\sum_{j=0}^{30}\left(\lfloor\dfrac{i*m+r}{2^j}\rfloor-\lfloor\dfrac{i*m+r}{2^{j+1}}\rfloor*2\right)\\ \sum_{j=0}^{30}\left(\sum_{i=0}^{\lfloor\frac{n-m}{r}\rfloor}\lfloor\dfrac{i*m+r}{2^j}\rfloor-\sum_{i=0}^{\lfloor\frac{n-m}{r}\rfloor}\lfloor\dfrac{i*m+r}{2^{j+1}}\rfloor*2\right) \end{align} \]
代码
C++
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
using i64 = long long;
struct Node {
    i64 cntU, cntR;
    i64 f;
    Node() : cntU(0), cntR(0) {
        f = 0;
    }
    Node(i64 cntU, i64 cntR, i64 f) : cntU(cntU), cntR(cntR), f(f) {}
    friend Node operator+(Node a, Node b) {
        Node res;
        res.cntU = a.cntU + b.cntU;
        res.cntR = a.cntR + b.cntR;
        res.f = a.f + b.f + a.cntU * b.cntR;
        return res;
    }
};
Node qpow(Node a, i64 b, Node res = {}) {
    while (b) {
        if (b & 1) res = res + a;
        a = a + a;
        b >>= 1;
    }
    return res;
}
i64 S(i64 P, i64 L, i64 R, i64 Q) {
    return (P * L + R) / Q;
}
Node get(i64 P, i64 Q, i64 R, i64 L, Node a, Node b) {
    if (L == 0LL) return Node();
    R %= Q;
    if (P >= Q) {
        return get(P % Q, Q, R, L, a, qpow(a, P / Q) + b);
    }
    i64 m = S(P, L, R, Q);
    if (m == 0LL) return qpow(b, L);
    i64 cnt = L - S(Q, m, -R - 1, P);
    return qpow(b, (Q - R - 1) / P) + a + get(Q, P, (Q - R - 1) % P, m - 1, b, a) + qpow(b, cnt);
}
Node solve(i64 a, i64 c, i64 b, i64 n, Node U, Node R) {
    Node res = qpow(U, b / c) + get(a, c, b, n, U, R);
    res.f += b / c;
    return res;
}
void solve() {
    i64 n, m, r;
    cin >> n >> m >> r;
    i64 ans = 0;
    Node U = {1, 0, 0}, R = {0, 1, 0};
    for (int i = 0 ; i <= 30 ; i++) {
        i64 A = solve(m, 1LL << i, r, (n - r) / m, U, R).f;
        i64 B = solve(m, 1LL << i + 1, r, (n - r) / m, U, R).f * 2;
        ans += A - B;
    }
    cout << ans << "\n";
}
int main() {
    int T;
    cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}