DDOSvoid's Blog

类欧几里得

发表于 更新于 分类于
本文字数： 5.8k 阅读时长 ≈ 5 分钟

简介

希望退役之前能用到一次

另外推导过程应该没什么用吧 = =

类欧几里得算法

大概就是用来求这三个东西的算法

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

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

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

其中 $a,b,c,n\ge 0$ 且为整数

首先是求 $f$

如果 $a\ge c$ 或者 $b\ge c$

那么我们能够得到

$f(a,b,c,n)=\sum_{i=0}^n\lfloor\frac{(a\%c)i+b\%c}{c}\rfloor+i\lfloor\frac{a}{c}\rfloor+\lfloor\frac bc\rfloor=f(a\%c,b\%c,c,n)+\frac{n*(n+1)}{2}\lfloor\frac ac\rfloor+(n+1)\lfloor\frac bc\rfloor$

直接递归计算即可

那么剩下的情况就是 $a<c$ 并且 $b<c$ 的情况了

我们来愉快的化式子，首先令 $m=\lfloor\frac{an+b}{c}\rfloor$

由于下次 $f(c,c-b-1,a,m-1)$ 一定变成第一种方式计算，所以复杂度还是 $O(\log n)$

边界条件 $a=0$

那么现在来求 $g$

对于 $a\ge c$ 或者是 $b\ge c$ 的情况

类似讨论能够得到

$g(a,b,c,n)=\frac{n(n+1)(2n+1)}{6}\lfloor \frac ac\rfloor+\frac{n(n+1)}{2}\lfloor \frac bc\rfloor+g(a\%c,b\%c,c,n)$

对 $a<c$ 并且 $b<c$ 的情况

后面是一个等差数列，所以 $g(a,b,c,n)=\frac12\sum_{j=0}^{m-1}(n-\lfloor\frac{cj+c-b-1}{a}\rfloor)(n+\lfloor\frac{cj+c-b-1)}{a}\rfloor+1)$

化简得到 $g(a,b,c,n)=\frac{nm(n+1)}{2}-\frac 1 2f(c,c-b-1,a,m-1)-\frac 1 2h(c,c-b-1,a,m-1)$

好像还要求 $h$

那么我们现在求 $h$

对于 $a\ge c$ 或者 $b \ge c$ 的情况

类似讨论能够得到

$h(a,b,c,n)= h(a\%c,b\%c,c,n) + \frac{n(n+1)(2n+1)}{6}\lfloor \frac ac\rfloor^2+(n+1)\lfloor \frac bc\rfloor^2+2\lfloor \frac ac\rfloor g(a\%c,b\%c,c,n)+2\lfloor \frac bc\rfloor f(a\%c,b\%c,c,n)+n(n+1)\lfloor \frac ac\rfloor\lfloor \frac bc\rfloor$

注意到 $\lfloor \frac {ai+b}{c}\rfloor^2$ 不太好整

我们考虑如何将平方去掉，能够发现 $2\sum_{i=1}^ni-n=n(n+1)-n=n^2$

贴一下代码，时间复杂度均为 $O(\log n)$，默认 $a,b,c,n$ 同阶

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
struct Likegcd {
    ll inv2, inv6; 
    void init() {
        inv2 = pow_mod(2, p - 2);
        inv6 = pow_mod(6, p - 2); 
    }

    inline ll S(ll n) {
        return n * (n + 1) % p * inv2 % p; 
    }

    inline ll S2(ll n) {
        return n * (2 * n + 1) % p * (n + 1) % p * inv6 % p; 
    }

    inline ll F(ll n) { return n * n % p; } 

    tuple<ll, ll, ll> calc(ll a, ll b, ll c, ll n) {
        ll f, g, h, tf, tg, th, m = (a * n + b) / c;
        if (!a) {
            f = (n + 1) * (b / c) % p;
            g = S(n) * (b / c) % p;
            h = F(b / c) * (n + 1) % p;
            return make_tuple(f, g, h);
        }
        if (a >= c || b >= c) {
            tie(tf, tg, th) = calc(a % c, b % c, c, n);
            f = (tf + S(n) * (a / c) + (n + 1) * (b / c)) % p; 
            g = (tg + S2(n) * (a / c) + S(n) * (b / c)) % p;
            h = (th + S2(n) * F(a / c) + (n + 1) * F(b / c)) % p;
            h = (h + 2 * (a / c) * tg + 2 * (b / c) * tf + n * (n + 1) % p * (a / c) % p * (b / c)) % p;
            return make_tuple(f, g, h); 
        }
        tie(tf, tg, th) = calc(c, c - b - 1, a, m - 1); 
        f = (n * m - tf) % p;
        g = (n * m % p * (n + 1) - tf - th) % p * inv2 % p;
        h = (n * m % p * (m + 1) - 2 * tf - 2 * tg - f) % p;
        return make_tuple(f, g, h); 
    }
} likegcd;

万能欧几里得算法

我们从类欧最经典的题目入手，考虑如何求这个东西 $\sum_{i=0}^n{\lfloor\frac{ai+b}{c}\rfloor}$

我们从 $ax+b$​ 这条直线出发，如果经过一个纵坐标为整数的点，我们记录一个 $U$，如果经过一个横坐标为整数的点，我们记录一个 $R$，如果经过一个整点，即横纵坐标都是整数的点，我们先记录一个 $U$，然后在记录一个 $R$

容易发现，每次我们记录一个 $R$ 的时候，就是我们走到了一个新的整数 $i$（感觉这么理解有一点怪），然后这个时候 ${\lfloor\frac{ai+b}{c}\rfloor}$ 就是之前记录的 $L$ 的个数，也就是说这整个操作序列的答案就是每个 $R$ 前面的 $U$ 的个数的和

那么这个东西有什么用呢，我们考虑将之前的东西转换为这样三元组 $(U,R,ans)$ 分别表示 $U$ 的个数和$R$​ 的个数以及答案

我们考虑三元组之间合并，容易得到 $(U,R,ans)\times (U’,R’,ans’)=(U+U’,R+R’,ans+ans’+U\times R’)$，并且这个东西是​​满足结合律的

现在我们考虑利用这个东西来求解 $\sum_{i=0}^n{\lfloor\frac{ai+b}{c}\rfloor}$，由于我们的算法是递归的，我们定义 $solve(a,b,c,n,A,B)$，其中 $A$ 表示遇到一个 $U$ 以后的变换，$B$ 表示遇到 $R$ 的变换

如果 $a\ge c$，那么容易得到两个 $R$ 之间一定有 $\lfloor\frac{a}{c}\rfloor$ 个连续的 $U$，所以我们将这些 $U$ 和这一个 $R$ 合并，得到 $solve(a\bmod c,b,c,A,A^{\lfloor\frac{a}{c}\rfloor}B)$

若 $a <c$​，通过 一系列操作，我们能够得到

具体证明参考：https://www.luogu.com.cn/blog/ix-35/solution-p5170

可以证明时间复杂度为 $O(F\log n)$，其中 $n,a,b,c$ 同阶，$F$ 为做一次乘法的时间复杂度

1
2
3
4
5
6
node calc(ll a, ll b, ll c, ll n, const node &A, const node &B) {
    ll m = (a * n + b) / c; 
    if (!m) return pow(B, n);
    if (a >= c) return calc(a % c, b, c, n, A, pow(A, a / c) * B);
    else return pow(B, (c - b - 1) / a) * A * calc(c, (c - b - 1) % a, a, m - 1, B, A) * pow(B, n - (m * c - b - 1) / a);
}

知道如何求 $\sum_{i=0}^n{\lfloor\frac{ai+b}{c}\rfloor}$ 后，我们来尝试一些更加复杂的形式，例如 $\sum_{i=0}^ni{\lfloor\frac{ai+b}{c}\rfloor}$

我们尝试记录 $U$​ 的数量 $x$，$R$ 的数量 $y$，${\lfloor\frac{ai+b}{c}\rfloor}$ 的和，$i$ 的和，$i{\lfloor\frac{ai+b}{c}\rfloor}$​ 的和

我们发现合并的时候大概就是右边的贡献应该为 $(i+y)({\lfloor\frac{ai+b}{c}\rfloor}+x)$​，我们将括号拆开即可

现在我们不妨来考虑更加一般化的形式，$\sum_{i=0}^n({\lfloor\frac{ai+b}{c}\rfloor})^{k_1}i^{k_2}$​​​，注意到右边的贡献变为 $({\lfloor\frac{ai+b}{c}\rfloor}+x)^a(i+y)^b$​​​，我们直接二项式展开，可以得到 $res.ans[a][b]=\sum_{i=0}^a\sum_{j=0}^{b}\binom{a}{i}\binom{b}{j}(u.x)^i(u.y)^jv.ans[a-i][b-j]$​​，由于我们还需要枚举 $a$​​ 和 $b$​​，所以总的时间复杂度为 $O(k^4\log n)$​​​

例题

  1. 简要题意：求 $\sum_{i=0}^{n}A^iB^{\lfloor\frac{ai+b}{c}\rfloor}$，其中 $A$ 和 $B$ 均为 $N\times N$ 的矩阵

    $a,b,c,n,\lfloor\frac{an+b}{c}\rfloor\le 10^{18},N\le 20$

    简要题解：我们考虑万能欧几里得算法，我们定义三元组 $(B^x,A^y,ans)$，容易得到 $(B^x,A^y,ans)\times (B^{x’},A^{y’},ans)=(B^{x+x’},A^{y+y’},A^yansB^x)$，时间复杂度 $O(N^3\log n)$​

    Loj 6440. 万能欧几里得

  2. 简要题意：有 $T$ 组询问，每组询问给出 $n,a,b,c,k_1,k_2$，求 $\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor^{k_1}i^{k_2}$​

    $T=1000,n,a,b,c\le 10^9,k_1+k_2\le 10$

    简要题解：考虑万欧，定义三元组 $(x,y,ans[a][b])$​

    容易得到转移 $res.ans[a][b]=\sum_{i=0}^a\sum_{j=0}^{b}\binom{a}{i}\binom{b}{j}(u.x)^i(u.y)^jv.ans[a-i][b-j]$，就是直接二项式展开，时间复杂度 $O(k^4\log n)$

    Loj 138. 类欧几里得算法