类欧几里德算法

定义

类欧几里德算法是洪华敦在 2016 年冬令营营员交流中提出的内容。

其本质可以理解为,使用一个类似辗转相除法的方法来进行函数求和。

引入

其中 是常数。需要一个 的算法。

这个式子和我们以前见过的式子都长得不太一样。带向下取整的式子容易让人想到数论分块,然而数论分块似乎不适用于这个求和。但是我们是可以做一些预处理的。

如果说 或者 ,意味着可以将 取模以简化问题:

那么问题转化为了 的情况。观察式子,你发现只有 这一个变量。因此要推就只能从 下手。在推求和式子中有一个常见的技巧,就是条件与贡献的放缩与转化。具体地说,在原式 中, 是条件,而 是对总和的贡献。

要加快一个和式的计算过程,所有的方法都可以归约为 贡献合并计算。但你发现这个式子的贡献难以合并,怎么办?将贡献与条件做转化 得到另一个形式的和式。具体地,我们直接把原式的贡献变成条件:

现在多了一个变量 ,既然算 的贡献不方便,我们就想办法算 的贡献。因此想办法搞一个和 有关的贡献式。这里有另一个家喻户晓的变换方法,笔者概括为限制转移。具体来说,在上面的和式中 限制 的上界,而 限制 的上界。为了搞 ,就先把 j 放到贡献的式子里,于是我们交换一下 的求和算子,强制用 限制 的上界。

这样做的目的是让 摆脱 的限制,现在 都被 限制,而贡献式看上去是一个条件,但是我们仍把它叫作贡献式,再对贡献式做变换后就可以改变 的限制关系。于是我们做一些放缩的处理。首先把向下取整的符号拿掉

然后可以做一些变换

最后一步,向下取整得到:

这一步的重要意义在于,我们可以把变量 消掉了!具体地,令 ,那么原式化为

这是一个递归的式子。并且你发现 分子分母换了位置,又可以重复上述过程。先取模,再递归。这就是一个辗转相除的过程,这也是类欧几里德算法的得名。

容易发现时间复杂度为

扩展

理解了最基础的类欧几里德算法,我们再来思考以下两个变种求和式:

推导 g

我们先考虑 ,类似地,首先取模:

接下来考虑 的情况,令 。之后的过程我会写得很简略,因为方法和上文略同:

这时我们设 ,可以得到

推导 h

同样的,首先取模:

考虑 的情况,.

我们发现这个平方不太好处理,于是可以这样把它拆成两部分:

这样做的意义在于,添加变量 的时侯就只会变成一个求和算子,不会出现 的形式:

接下来考虑化简前一部分:

因此

在代码实现的时侯,因为 个函数各有交错递归,因此可以考虑三个一起整体递归,同步计算,否则有很多项会被多次计算。这样实现的复杂度是 的。

实现

模板题代码实现
 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
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int P = 998244353;
int i2 = 499122177, i6 = 166374059;

struct data {
  data() { f = g = h = 0; }

  int f, g, h;
};  // 三个函数打包

data calc(int n, int a, int b, int c) {
  int ac = a / c, bc = b / c, m = (a * n + b) / c, n1 = n + 1, n21 = n * 2 + 1;
  data d;
  if (a == 0) {  // 迭代到最底层
    d.f = bc * n1 % P;
    d.g = bc * n % P * n1 % P * i2 % P;
    d.h = bc * bc % P * n1 % P;
    return d;
  }
  if (a >= c || b >= c) {  // 取模
    d.f = n * n1 % P * i2 % P * ac % P + bc * n1 % P;
    d.g = ac * n % P * n1 % P * n21 % P * i6 % P + bc * n % P * n1 % P * i2 % P;
    d.h = ac * ac % P * n % P * n1 % P * n21 % P * i6 % P +
          bc * bc % P * n1 % P + ac * bc % P * n % P * n1 % P;
    d.f %= P, d.g %= P, d.h %= P;

    data e = calc(n, a % c, b % c, c);  // 迭代

    d.h += e.h + 2 * bc % P * e.f % P + 2 * ac % P * e.g % P;
    d.g += e.g, d.f += e.f;
    d.f %= P, d.g %= P, d.h %= P;
    return d;
  }
  data e = calc(m - 1, c, c - b - 1, a);
  d.f = n * m % P - e.f, d.f = (d.f % P + P) % P;
  d.g = m * n % P * n1 % P - e.h - e.f, d.g = (d.g * i2 % P + P) % P;
  d.h = n * m % P * (m + 1) % P - 2 * e.g - 2 * e.f - d.f;
  d.h = (d.h % P + P) % P;
  return d;
}

int T, n, a, b, c;

signed main() {
  scanf("%lld", &T);
  while (T--) {
    scanf("%lld%lld%lld%lld", &n, &a, &b, &c);
    data ans = calc(n, a, b, c);
    printf("%lld %lld %lld\n", ans.f, ans.h, ans.g);
  }
  return 0;
}