8.12.10 ACM-ICPC数学 数论 类欧几里德算法

8.12.10 ACM-ICPC数学 数论 类欧几里德算法

定义

类欧几里德算法由洪华敦在2016年冬令营中提出，本质上是一个使用类似于辗转相除法的方法进行函数求和的算法。

引入

考虑求和函数：

其中 𝑎,𝑏,𝑐,𝑛a,b,c,n 是常数，我们需要一个 𝑂(log⁡𝑛)O(logn) 的算法来计算这个求和。

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

如果 𝑎≥𝑐a≥c 或者 𝑏≥𝑐b≥c，意味着可以将 𝑎,𝑏a,b 对 𝑐c 取模以简化问题：

这等于：

问题转化为 𝑎<𝑐,𝑏<𝑐a<c,b<c 的情况。我们发现只有变量 𝑖i，因此要推就只能从 𝑖i 下手。在推求和式子中有一个常见的技巧，就是条件与贡献的放缩与转化。具体地说，在原式 𝑓(𝑎,𝑏,𝑐,𝑛)f(a,b,c,n) 中，0≤𝑖≤𝑛0≤i≤n 是条件，而 ⌊𝑎𝑖+𝑏𝑐⌋⌊cai+b​⌋ 是对总和的贡献。

贡献与条件的转化

我们把原式的贡献变成条件：

我们交换求和的顺序，让 𝑗j 摆脱 𝑖i 的限制：

这时我们对条件与贡献做进一步的变换：

向下取整得到：

这一步的重要意义在于，我们可以把变量 𝑖i 消掉了！

递归式

令 𝑚=⌊𝑎𝑛+𝑏𝑐⌋m=⌊can+b​⌋，那么原式化为：

这是一个递归的式子。并且你发现 𝑎,𝑐a,c 分子分母换了位置，又可以重复上述过程。先取模，再递归。这就是一个辗转相除的过程，这也是类欧几里德算法的得名。容易发现时间复杂度为 𝑂(log⁡𝑛)O(logn)。

扩展

理解了最基础的类欧几里德算法后，我们可以探索它在求解其他类型函数（如 𝑔(𝑎,𝑏,𝑐,𝑛)g(a,b,c,n) 和 ℎ(𝑎,𝑏,𝑐,𝑛)h(a,b,c,n)）中的应用，进一步扩展其应用场景和解决问题的能力。

#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;
}