跳转至

数论分块

数论分块可以快速计算一些形如

i=1nf(i)g(ni)

的和式.如果可以在 O(1) 时间内计算出 i=lrf(i) 或已经预处理出 f 的前缀和时,数论分块就可以在 O(n) 时间内计算出上述和式的值.

数论分块常与 莫比乌斯反演 等技巧结合使用.

思路

首先,本文通过一个简单的例子来说明数论分块的思路.假设要计算如图所示的双曲线下的整点个数:

双曲线下整点

这相当于在计算和式:

i=11111i.

这是前文所示和式在 f(k)=1, g(k)=k 时的特殊情况.

最简单的做法当然是逐列计算然后求和,但是这样需要计算第 i=1,2,,11 列中每列整点的个数.观察图示可以发现,这些整点列的可以分成 5 块,每块内点列的高度是一致的,形成一个矩形点阵.所以,只要能够知道这些块的宽度,就能够通过计算这些矩形块的大小快速完成统计.

这就是整除分块的基本思路.

性质

本节讨论关于双曲线 y=nx 下整点分块的若干结论.具体地,需要将 1n 之间的整数按照 ni 的取值分为若干块.设

Dn={ni:1in, iN+}.

这就是 ni 所有可能的取值集合.

首先,这样的不同取值只有 O(n) 个.所以,数论分块得到的块的数目只有 O(n) 个.

性质 1

|Dn|2n

证明

分两种情况讨论:

  • in 时,i 的取值至多 n 个,所以 ni 的取值也至多只有 n 个.
  • i>n 时,nini<n 也至多只有 n 个取值.

因此,所有可能取值的总数 |Dn|2n

然后,单个块的左右端点都很容易确定.

性质 2

对于 dDn,所有满足 ni=d 的整数 i 的取值范围为

nd+1+1ind.
证明

因为 ni=d 相当于不等式

dni<d+1.

这进一步等价于

nd+1<ind.

利用 iN+ 这一点,可以对该不等式取整,它就等价于

nd+1+1ind.

这一性质还体现了图像的对称性:每个块的右端点(前文图中的绿点)的集合恰为 Dn.这很容易理解,因为整个图像关于直线 y=x 是对称的.

除了这两条性质外,集合 Dn 还具有良好的递归性质.这基于如下引理:

引理

对于 a,b,cN+,有

a/bc=abc.
证明

a 对于 b 做带余除法有 a=qb+r,其中 q,rNr<b.需要证明的是

qc=qc+rbc.

这等价于

qcqc+rbc<qc+1.

左侧不等式是显然的.关键是要证明右侧不等式.为此,令 qc 做带余除法,就有

q=qcc+r,

其中,0rc1.所以,有

qqcc+c1q+1cqc+1.

进而,利用 r<b,就有

qc+rbc<q+1cqc+1.

这就证明右侧不等式,进而完成本命题的证明.

这一引理经常出现在数论分块的各种应用中.利用该引理可以得到如下性质:

性质 3

对于 mD(n),有 D(m)D(n)

证明

m=nk,那么,因为对于所有 iN+,都有

mi=n/ki=nkiD(n),

所以,D(m)D(n)

前文已经提到,D(n) 既是每块中 ni 的取值集合,也是每块的右端点集合.这意味着,如果要递归地应用数论分块(即函数在 n 处的值依赖于它在 mD(n){n} 处的值),那么在整个计算过程中所涉及的取值集合和右端点集合,其实都是 D(n).一个典型的例子是 杜教筛

过程

利用上一节叙述的结论,就得到数论分块的具体过程.

为计算和式

i=1nf(i)g(ni)

的取值,可以依据 ni 的取值将标号 i=1,2,,n 分块.因为 ni 取值相同的标号是一段连续整数 [l,r],所以该块中和式的取值为

(i=lrf(i))g(nl).

为了快速计算该和式,通常需要能够快速计算左侧关于 f 的和式.有些时候,该和式的表达式是已知的,可以在 O(1) 时间内完成单次计算;有些时候,可以预处理出它的前缀和,仍然可以在 O(1) 时间内完成单次查询.

在顺次计算每块左右端点时,当前块的左端点 l 就等于前一块的右端点再加 1,而当前块的右端点就等于 nn/l.由此,可以得到如下伪代码:

Algorithm Sum(f,g,n):Input. n, s(k)=i=1kf(k), g(k).Output. S(n)=i=1nf(i)g(n/i).Method.1l12result03while ln do4rnn/l5resultresult+(s(r)s(l1))g(nl)6lr+17end while8return result

假设单次计算 s() 的时间复杂度为 O(1),则整个过程的时间复杂度为 O(n)

扩展

前文讨论的是数论分块的最常见也最基本的形式.本节进一步讨论数论分块的扩展形式.

向上取整的数论分块

数论分块可以用于计算含有向上取整的和式:

i=1nf(i)g(ni).

因为 ni=n1i+1,该和式可以转化为向下取整的情形:

f(n)g(1)+i=1n1f(i)g(n1i+1).

注意到求和的上限发生了变化,以及 i=n 时单独的一项.

多维数论分块

数论分块还可以用于处理包含不只有一个取整式的和式:

i=1nf(i)g(n1i,n2i,,nmi).

为了应用数论分块的思想,需要保证每块中所有取整式 n1i,n2i,,nmi 的取值都不发生变化,也就是说,多维的块应当是所有一维的块的交集.为此,对于已知的左端点 l,相应的右端点为

min{n1n1/l,n2n2/l,,nmnm/l}.

也就是说,对于所有一维分块的右端点取最小值,作为多维分块的右端点.可以借助下图理解:

多维数论分块图解

较为常见的是二维形式,此时可将前述伪代码中 rnn/l 替换成

rmin{n1n1/l,n2n2/l}.

任意指数数论分块

数论分块可以用于含有任意指数的取整式的和式计算:

i=1nα/βf(i)g(nαiβ).

其中,α,β 为正实数.本文讨论的基本形式中,α=β=1

性质

对于正整数 n 和正实数 α,β,设

Dn,α,β={nαiβ:i=1,2,,nα/β}.

那么,有

  1. |Dn,α,β|2nα/(1+β)
  2. 对于 dDn,α,β,使得 nαiβ=d 成立的 i 的取值范围为

    nα/β(d+1)1/β+1inα/βd1/β.
证明

对于第一点,分两种情况:

  • inαiβ 时,有 inα/(1+β),所以 nαiβ 至多有 nα/(1+β) 种取值.
  • i>nαiβ 时,有 i>nα/(1+β),进而有 nαiβ<nα/(1+β),所以 nαiβ 也至多只有 nα/(1+β) 种取值.

综合两种情形,就有 |Dn,α,β|2nα/(1+β)

对于第二点,nαiβ=d 就等价于

dnαiβ<d+1nα/β(d+1)1/β<inα/βd1/β.

对该不等式取整,就得到第二个命题.

利用这些性质,就可以在 O(nα/(1+β)) 的时间复杂度下实现对任意指数的数论分块.

例子

例如,对于 α=β=1/2 时的如下和式

i=1nf(i)g(ni),

可以通过数论分块在 O(n1/3) 时间内解决.已知块的左端点为 l 时,可以计算右端点为 r=nn/l2

例题

UVa11526 H(n)

T 组数据,每组一个整数 n.对于每组数据,输出 i=1nni

解答

根据前文分析,可以对于每一块相同的 ni 一起计算.时间复杂度为 O(Tn)

实现
 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
#include <iostream>

long long H(int n) {
  long long res = 0;  // 储存结果
  int l = 1, r;       // 块左端点与右端点
  while (l <= n) {
    r = n / (n / l);  // 计算当前块的右端点
    // 累加这一块的贡献到结果中。乘上 1LL 防止溢出
    res += 1LL * (r - l + 1) * (n / l);
    l = r + 1;  // 左端点移到下一块
  }
  return res;
}

int main() {
  std::ios::sync_with_stdio(false);
  std::cin.tie(nullptr);
  int t, n;
  std::cin >> t;
  while (t--) {
    std::cin >> n;
    std::cout << H(n) << '\n';
  }
  return 0;
}
Codeforces 1954E Chain Reaction

有一排 n 只怪兽,每只怪兽初始血量为 ai.一次攻击会使一段连续的存活的怪兽血量减 k,血量不大于 0 视作死亡.对于所有 k 求出击杀所有怪兽所需攻击次数.其中,n,ai105

解答

a0=0.假设击杀所有前 (i1) 只怪兽需要 T(k,i1) 次攻击,第 i 只怪兽的血量为 ai.由于击杀第 (i1) 只怪兽时,需要攻击它 ai1/k 次,这些攻击都可以延伸到第 i 只怪兽.因此,要击杀第 i 只怪兽,只需要再攻击 max{0,ai/kai1/k} 次.由此,总的攻击次数为

T(k,n)=i=1nmax(0,aikai1k).

由于题目涉及的 n,k 都比较大,对每个 k 分别计算该和式并不可行.可以考虑对每个 i=1,2,,n,都维护数列 {T(k,i)}k.初始时,设 T(k,0)0.假设数列 {T(k,i1)}k 已知,考虑如何对它进行修改才能得到数列 {T(k,i)}k.根据前文分析,只需要对数列的第 k 项增加 max(0,aikai1k) 即可.利用二维数论分块,可以将这一修改操作拆分成 O(ai1+ai) 段区间修改操作,且每段区间上增加的值是固定的.最后得到的数列 {T(k,n)}k 就是答案.

由于题目涉及一系列区间加操作,且查询只发生所有修改完成后.所以,可以通过维护差分序列进行区间加修改,最后通过求前缀和得到所求数列.总的时间复杂度为 O(ai).本题也存在其他解法.

实现
 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
#include <algorithm>
#include <iostream>

constexpr int N = 100009;
int n, a[N], maxn;
long long ans[N];

int main() {
  std::ios::sync_with_stdio(false);
  std::cin.tie(nullptr);
  std::cin >> n;
  for (int i = 1; i <= n; ++i) {
    std::cin >> a[i];
    maxn = std::max(maxn, a[i]);
  }
  for (int i = 0; i < n; ++i)
    for (int l = 1, r;; l = r + 1) {
      r = std::min(l < a[i] ? (a[i] - 1) / ((a[i] - 1) / l) : N,
                   l < a[i + 1] ? (a[i + 1] - 1) / ((a[i + 1] - 1) / l)
                                : N);  // 二维数论分块
      if (r == N) break;
      int x = (a[i + 1] - 1) / l - std::max(a[i] - 1, 0) / l;
      if (x > 0) ans[l] += x, ans[r + 1] -= x;  // 累加贡献
    }
  ++ans[0];  // ⌈a/l⌉=(a-1)/l+1的式子当a=0时不成立,需要修正
  for (int i = 1; i <= maxn; ++i)
    std::cout << (ans[i] += ans[i - 1]) << " \n"[i == maxn];
  return 0;
}

习题

参考资料与注释