NeoMind

Liuguang_Ji的算法笔记

数论函数专题

数论函数专题

1. 数论函数

  • 以正整数集的一个子集 DD 作为定义域的函数称为数论函数,值域 YCY \in C

1.1 积性函数和加性函数

若一个数论函数 f(x)f(x) 满足定义域 DD 对乘法封闭,且 m,nD\forall m, n \in D(m,n)=1(m, n) = 1 时,总有 f(mn)=f(n)f(m)f(mn) = f(n)f(m) 那么称这个数论函数为积性函数 更进一步的若满足 m,nD\forall m, n \in D,总有 f(mn)=f(n)f(m)f(mn) = f(n)f(m) 则称该数论函数为完全积性函数

若一个数论函数 f(x)f(x) 满足定义域 DD 对乘法封闭,且 m,nD\forall m, n \in D(m,n)=1(m, n) = 1 时,总有 f(mn)=f(n)+f(m)f(mn) = f(n) + f(m) 那么称这个数论函数为加性函数 更进一步的若满足 m,nD\forall m, n \in D,总有 f(mn)=f(n)+f(m)f(mn) = f(n) + f(m) 则称该数论函数为完全加性函数

积性函数性质

  1. 积性函数的积和商,实数幂依然是积性函数

  2. f(n)f(n) 是积性函数的一个充要条件是 f(1)=1f(1) = 1,且对于 n>1n > 1 时标准素因数分解 n=p1α1p2α2...psαsn = p_1^{\alpha_1} p_2^{\alpha_2} … p_s^{\alpha_s}f(n)=i=1sf(piαi)f(n) = \prod_{i = 1}^s f(p_i^{\alpha_i})

  3. 类似地,f(n)f(n) 是完全积性函数的一个充要条件是 f(1)=1f(1) = 1,且对于 n>1n > 1 时标准素因数分解 n=p1α1p2α2...psαsn = p_1^{\alpha_1} p_2^{\alpha_2} … p_s^{\alpha_s}f(n)=i=1sf(pi)αif(n) = \prod_{i = 1}^s f(p_i)^{\alpha_i}

1.2 常见数论函数基本性质介绍

1.2.1 元函数 ε\varepsilon

元函数 ε\varepsilon 定义为 ε(n)=[1n]=[n=1]\varepsilon(n) = [\frac{1}{n}] = [n = 1],这表明函数取 11 当且仅当自变量为 11,其他情况函数值为 00

  • 元函数是完全积性函数

1.2.2 幂函数 IkI_k

幂函数 IkI_k 定义为 Ik(n)=nkI_k(n) = n^k 特别地, 当 k=0k = 0 时为常数函数 1(n)=11(n) = 1k=1k = 1 时为恒等函数 I(n)=nI(n) = n

  • 幂函数是完全积性函数

1.2.3 因数个数函数 d(n)/τ(n)d(n) / \tau(n)

因数个数函数d(n)/τ(n)d(n) / \tau(n) 定义为 dn1(d)\sum_{d \mid n} 1(d),即 nn 的全体因数的个数

数论性质
  • 因数个数函数是积性函数
  • n=p1α1p2α2...psαsn = p_1^{\alpha_1} p_2^{\alpha_2} … p_s^{\alpha_s},则 d(n)=i=1s(αi+1)d(n) = \prod_{i = 1}^s (\alpha_i + 1)
分析学性质
  • d(n)d(n) 函数的值域相当离散,对于 f(n)=max1ind[i]f(n) = max_{1 \leq i \leq n}{d[i]} 这个函数来说,它的阶小于任意幂函数但大于对数函数,但在算法竞赛的小范围里 f(104)=64,f(105)=128,f(106)=240f(10^4) = 64, f(10^5) = 128, f(10^6) = 240 f(107)=448,f(109)=1344,f(1012)=6720f(10^7) = 448, f(10^9) = 1344, f(10^{12}) = 6720
  • d(n)d(n) 的前 nn 项和却比较稳定,阶约为 O(nlogn)O(nlogn)

1.2.4 因数函数 σk(n)\sigma_k(n)

因数幂和函数 σk(n)\sigma_k(n),定义为 dknI(d)\sum_{d_k \mid n} I(d),即 nn 的全体因数之和 特别的, 当 k=0k = 0 时,σ0(n)\sigma_0(n) 就是刚刚介绍的因数个数函数 d(n)d(n)k=1k = 1 时,σ1(n)\sigma_1(n) 记作 σ(n)\sigma(n),称为因数和函数

  • 因数幂和函数是积性函数
  • n=p1α1p2α2...psαsn = p_1^{\alpha_1} p_2^{\alpha_2} … p_s^{\alpha_s},则 σ(n)=i=1spiαi1pi1\sigma(n) = \prod_{i = 1}^s \frac{p_i^{\alpha_i} - 1}{p_i - 1}

1.2.5 欧拉函数 ϕ(n)\phi(n)

欧拉函数 ϕ(n)\phi(n) 定义为 ϕ(n)=i=1n[(i,n)=1]\phi(n) = \sum_{i = 1}^n [(i, n) = 1],即不超过 nnnn 互质的正整数的个数

数论性质
  • 欧拉函数是积性函数
  • 欧拉函数的一个可以基于容斥得到的公式 ϕ(n)=npn(1p1)\phi(n) = n\prod_{p \mid n}(1 - p^{-1})
  • 依赖于欧拉函数的欧拉定理:在 (a,n)=1(a, n) = 1 时 总有 aϕ(n)1(modn)a^{\phi(n)} \equiv 1 \pmod n
  • nn 互质的数字之和为 nϕ(n)2\frac{n\cdot \phi(n)}{2}
分析学性质
  • 欧拉函数的上界显然是 nn,且存在无穷项接近它(ϕ(p)=p1\phi(p) = p - 1),一个下界是 nloglogn\frac{n}{loglogn}
  • 欧拉函数前 nn 项和 S(n)=3π2n2+O(nlogn)S(n) = \frac{3}{\pi^2}n^2 + O(nlogn),所以欧拉函数的平均阶是 Θ(n)\Theta(n)
  • 欧拉函数的二阶复合 ϕ(ϕ(n))\phi(\phi(n)),上界阶是 nloglognlogn\frac{nloglogn}{logn},下界阶是 n(loglogn)2\frac{n}{(loglogn)^2},平均阶 nlogn\frac{n}{logn}

1.2.6 本质不同素因数个数函数 ω(n)\omega(n)

本质不同素因数个数函数 ω(n)\omega(n) 定义为 ω(n)=pn1(p)\omega(n) = \sum_{p \mid n} 1(p),即在标准素因数分解下 n>1,n=p1α1p2α2...psαsn > 1, n = p_1^{\alpha_1} p_2^{\alpha_2} … p_s^{\alpha_s},其中 pip_i 互不相同,则

ω(n)={sif n>10if n=1\omega(n) = \begin{cases} s & \text{if } n > 1\ 0 & \text{if } n = 1 \end{cases}

数论性质
  • ω(n)\omega(n) 是加性函数而不是积性函数
  • ω(ab)=ω(a)+ω(b)ω(gcd(a,b))=ω(lcm(a,b))\omega(ab) = \omega(a) + \omega(b) - \omega(gcd(a, b)) = \omega(lcm(a, b))
分析学性质
  • ω(n)\omega(n) 的下界显然是 11,且存在无穷项接近它(ω(p)=1\omega(p) = 1),上界阶是 lognloglogn\frac{logn}{loglogn},平均阶是 loglognloglogn

1.2.7 全部素因子个数函数 Ω(n)\Omega(n)

在标准素因数分解下 n>1,n=p1α1p2α2...psαsn > 1, n = p_1^{\alpha_1} p_2^{\alpha_2} … p_s^{\alpha_s},其中 pip_i 互不相同,全部素因子个数函数 Ω(n)\Omega(n) 定义为

Ω(n)={i=1sαiif n>11if n=1\Omega(n) = \begin{cases} \sum_{i = 1}^s \alpha_i & \text{if } n > 1 \ 1 & \text{if } n = 1 \end{cases}

数论性质
  • Ω(n)\Omega(n) 是完全加性函数

1.2.8 莫比乌斯函数 μ(n)\mu(n)

莫比乌斯函数 μ(n)\mu(n) 定义为 μ(n)={1if n=10if d>1s.t.d2n(1)w(n)others\mu(n) = \begin{cases} 1 & \text{if } n = 1 \ 0 & \text{if } \exists d > 1 ::s.t.:d^2 \mid n \ (-1)^{w(n)} & \text{others} \end{cases}

即假设 n2n \geq 2 的标准素因数分解 n=p1α1p2α2...psαsn = p_1^{\alpha_1} p_2^{\alpha_2} … p_s^{\alpha_s},其中 pip_i 互不相同,若 1is\exists 1 \leq i \leq s 使得 αi2\alpha_i \geq 2 则函数值为 00,否则取值为 ±1\pm 1,具体与素因子个数奇偶有关

数论性质
  • μ(n)\mu(n) 是积性函数
  • μ(n)\mu(n) 在狄利克雷卷积里具有非常重要的地位
  • 其中最重要的容斥性质是 ε(n)=[n=1]=dnμ(d)\varepsilon(n) = [n = 1] = \sum_{d \mid n}\mu(d)
  • μ(n)2=d2nμ(d)\mu(n)^2 = \sum_{d^2 \mid n}\mu(d)
分析学性质
  • μ(n)\mu(n) 的值域只有 0,±1{0, \pm1}
  • limnμ(n)n=0\lim_{n \to \infty} \frac{\mu(n)}{n} = 0
  • μ(n)2=6π2n+O(n)\mu(n)^2 = \frac{6}{\pi^2}n + O(\sqrt n)

1.3 从狄利克雷Dirichlet卷积到莫比乌斯变换

1.3.1 狄利克雷Dirichlet卷积及数论函数代数结构

定义一种在数论函数集合 AA 上的二元运算 *fg(n)=dnf(d)g(nd)f * g(n) = \sum_{d \mid n}f(d)g(\frac{n}{d})

容易证明,Dirichlet卷积在数论函数集合 AA封闭,且具有交换律结合律和与数论函数加法的分配律 Dirichlet卷积意义下,数论函数的单位元元函数 ε\varepsilon 加法意义下,数论函数的零元O(n)=0O(n) = 0 于是容易发现: 全体数论函数 AA 在狄利克雷卷积下构成一个含幺交换环

但想要构成群还需要逆元的性质,下面这个定理告诉了我们数论函数关于卷积存在逆元的充要条件 定理:数论函数 ff 关于Dirichlet卷积可逆当且仅当 f(1)0f(1) \neq 0

于是可逆数论函数集合 U={fAf(1)0}U = {f \in A | f(1) \neq 0} 关于Dirichlet卷积构成阿贝尔群 又由于积性函数总是满足 f(1)=1f(1) = 1 故全体积性函数集合 MMUU 的一个子群

1.3.2 莫比乌斯变换及反演

我们定义一种变换 F(n)=dnf(d)F(n) = \sum_{d \mid n}f(d) 称为莫比乌斯变换 从狄利克雷Dirichlet卷积的更高观点来看,这只是 g(n)=1(n)g(n) = 1(n) 的卷积特例,但这个变换本身具有很重要的价值,这个变换的计算过程在算法意义下又称作狄利克雷前缀和,后面预处理部分我们会研究这种变换的计算方法 这个变换的逆变换就称为莫比乌斯逆变换莫比乌斯反演

定理:对于两个数论函数 ffFF,满足莫比乌斯变换 F(n)=dnf(d)F(n) = \sum_{d \mid n}f(d) 的充要条件是满足莫比乌斯反演公式 f(n)=dnμ(d)F(nd)f(n) = \sum_{d \mid n}\mu(d)F(\frac{n}{d})

  • 初等数论教材上是利用莫比乌斯的容斥性质给出的初等证明
  • 但若使用Dirichlet卷积在代数结构上的性质这个定理是显然的,高观点有高观点的好处

我们下面给出之前所有提及的数论函数的莫比乌斯变换及其反演

f(n)f(n) F(n)F(n) 公式 特性
μ(n)\mu(n) ε(n)\varepsilon(n) ε(n)=dnμ(d)\varepsilon(n) = \sum_{d \mid n} \mu(d)μ(n)=dnμ(d)ε(nd)\mu(n) = \sum_{d \mid n} \mu(d)\varepsilon(\frac{n}{d}) 莫比乌斯函数 μ(n)\mu(n) 是常数函数逆元
ϕ(n)\phi(n) I(n)I(n) I(n)=dnϕ(d)I(n) = \sum_{d \mid n} \phi(d)ϕ(n)=dnμ(d)I(nd)\phi(n) = \sum_{d \mid n} \mu(d)I(\frac{n}{d}) 这个公式也常常称为欧拉变换及反演
1(n)1(n) d(n)d(n) d(n)=dn1(d)d(n) = \sum_{d \mid n} 1(d)1(n)=dnμ(d)d(nd)1(n) = \sum_{d \mid n} \mu(d)d(\frac{n}{d}) 因数个数函数定义
I(n)I(n) σ(n)\sigma(n) σ(n)=dnI(d)\sigma(n) = \sum_{d \mid n} I(d)I(n)=dnμ(d)σ(nd)I(n) = \sum_{d \mid n} \mu(d)\sigma(\frac{n}{d}) 因数和函数定义
μ(n)n\frac{\mu(n)}{n} ϕ(n)n\frac{\phi(n)}{n} ϕ(n)n=dnμ(d)d\frac{\phi(n)}{n} = \sum_{d \mid n}\frac{\mu(d)}{d}μ(n)=dndμ(d)ϕ(nd)\mu(n) = \sum_{d \mid n} d \mu(d)\phi(\frac{n}{d}) 是欧拉变换的推论

2. 数论函数预处理

接下来的部分研究数论函数的预处理,包括单点求值,前缀求值,单点前缀和,前缀和数组,狄利克雷前缀和等,是处理数论函数问题的武器库

2.1 素数及素数筛

预处理素数是一个重要的基本问题,后面很多的预处理依赖于素数筛这个基本的技能

2.1.1 单个数的素数检测

试除法

由于最小质因子是小于 n\sqrt n 的,我们只需要判断 22n\sqrt n 的数字是否是 nn 的因数,时间复杂度 O(n)O(\sqrt n)

bool is_prime(ll n)
{
    if(n <= 1)
        return 0;
    for(ll i = 2; i * i <= n; ++i)
    {
        if(n % i == 0)
            return 0;
    }
    return 1;
}
cpp

另外有一种卡常版本的这里不展示

Miller_Rabin素性检测

一种基于随机化的素数检测方法,时间复杂度 O(k(logn)3)O(k(logn)^3)kk 是检测次数,至少在 101810^{18} 内不担心错误性

vector<ll> p = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};

bool Miller_Rabin(lll n)
{
    if (n < 3 || n % 2 == 0)
        return n == 2;
    lll u = n - 1, t = 0;
    while (u % 2 == 0)
        u /= 2, ++t;
    for (auto a : p)
    {
        if (n == a)
            return 1;
        if (n % a == 0)
            return 0;
        lll v = qpow(a, u, n);
        if (v == 1)
            continue;
        lll s = 1;
        for (; s <= t; ++s)
        {
            if (v == n - 1)
                break;
            v = v * v % n;
        }
        if (s > t)
            return 0;
    }
    return 1;
}
cpp

2.2.2 素数筛

更加有价值的问题是如何获得从 11nn 的素数 一个基本知识点是素数定理,定义 π(x)\pi(x) 是小于等于 xx 的素数个数函数,则 π(x)xlnx\pi(x) \sim \frac{x}{lnx}

大概可以估计素数密度和个数

首先这个问题如果我们使用单个素数判别方法,时间复杂度是 O(nn)O(n\sqrt n) 是不优的,因为考查连续数字素性时有很多新性质

朴素筛法

筛法引申出的核心思想就是调和级数枚举,最简单常见的例子是我们枚举小于等于 nn 的每个数及其倍数,这样的时间复杂度是 O(n(11+12+...+1n))=O(nlogn)O(n(\frac{1}{1} + \frac{1}{2} + … + \frac{1}{n})) = O(nlogn)

所以每个数除了自己以外的所有倍数自然就可以被我们筛掉了,时间复杂度 O(nlogn)O(nlogn)

埃式筛Eratosthenes

我们注意到一个数若没有被更小的数筛过,那么这个数就是素数了,所以在从小到大的遍历时,我们天然的就得到了素数,而事实上,我们只用素数去枚举倍数就把所有的合数筛掉了,就是一个小优化

时间复杂度中调和级数的项就变成了全体素数倒数和,我们不加证明的给出这个级数的阶大概是 loglognloglogn,于是总时间复杂度就是 O(nloglogn)O(nloglogn)

bitset<N> vis;//判断素性的数组
vector<ll> pr;//质数数组

void EratosthenesSieve(ll n)
{
	vis.reset();
	pr.clear();
	vis[0] = vis[1] = 1;
	for (int i = 2; i <= n; ++i)//不在这里使用专门存储素数的数组的话到sqrt(n)即可
	{
		if (vis[i])
			continue;
		pr.push_back(i);
		for (int j = i; j <= n; j += i)
			vis[j] = 1;
	}
}
cpp

欧拉筛EulerSieve

为什么埃式筛不是线性?因为多素因子的合数被多次筛掉了,欧拉筛通过控制每个合数都被最小的质因子筛掉实现了 O(n)O(n) 的线性预处理质数

bitset<N> vis;
vector<int> pr;

void EulerSieve(ll n)//欧拉筛
{
	pr.clear();
	vis.reset();
	vis[0] = vis[1] = 1;
	for (ll i = 2; i <= n; ++i)
	{
		if (!vis[i])
			pr.push_back(i);
		for (ll j = 0; j < pr.size() && i * pr[j] <= n; ++j)
		{
			vis[i * pr[j]] = 1;
			if (i % pr[j] == 0)
				break;
		}
	}
}
cpp

2.2 数论函数朴素预处理

在朴素预处理部分,主要运用暴力计算,分解因数质因数,调和级数枚举,素数筛,组合数等基本的数学方法做计算

其中需要单独拿出来的方法性知识是欧拉筛一般化预处理积性函数 我们知道对于积性函数 f(n)f(n),首先 f(1)=1f(1) = 1

n>1n > 1 时标准素因数分解 n=p1α1p2α2...psαsn = p_1^{\alpha_1} p_2^{\alpha_2} … p_s^{\alpha_s}f(n)=i=1sf(piαi)f(n) = \prod_{i = 1}^s f(p_i^{\alpha_i})

所以在素数筛预处理 f(n)f(n) 时需要同步维护 g(n)=f(pα)g(n) = f(p^{\alpha}),其中 pp 是最小质因子,α\alpha 是最小质因子对于指数,如果这一项是可以在素数筛同步计算的,那么这个积性函数就可以 O(n)O(n) 预处理了

2.2.0 元函数 ε(n)\varepsilon(n) 和幂函数 Ik(n)I_k(n)

求值是比较显然的,即使是高次幂也可以快速幂或者欧拉公式降幂

  • 对于预处理 1n1-nIkI_k 也可以使用欧拉筛处理,由于 IkI_k 是完全积性函数,在素数处使用快速幂,其余数字使用乘积计算,可以做到 O(n)O(n) 的预处理(素数快速幂部分两个对数可以消掉)

下面主要针对前缀和,常见函数主要有以下求和式 i=1nε(i)=I0(n)=1i=1nI0(i)=I1(n)=ni=1nI1(i)=i=1ni=n(n+1)2i=1nI2(i)=i=1ni2=n(n+1)(2n+1)6i=1nI3(i)=i=1ni3=n2(n+1)24\sum_{i = 1}^n\varepsilon(i) = I_0(n) = 1 \ \sum_{i = 1}^nI_0(i) = I_1(n) = n \ \sum_{i = 1}^nI_1(i) = \sum_{i = 1}^ni = \frac{n(n + 1)}{2} \ \sum_{i = 1}^nI_2(i) = \sum_{i = 1}^ni^2 = \frac{n(n + 1)(2n + 1)}{6} \ \sum_{i = 1}^nI_3(i) = \sum_{i = 1}^ni^3 = \frac{n^2(n + 1)^2}{4}

对于一般化的 i=1nIk(i)=i=1nik\sum_{i = 1}^nI_k(i) = \sum_{i = 1}^ni^k

可以使用拉格朗日插值获取各项系数或值

2.2.1 最小质因子函数 minp(n)minp(n) 及其扩展

欧拉筛获取即可,时间复杂度 O(n)O(n)

  • 在可素数筛范围内预处理 minpminp,可以实现 O(logn)O(logn) 的分解质因数(注意可处理值域只到大概 10810^8

扩展:

  1. minexp(n)minexp(n) 最小质因子指数函数可以类似处理,将在预处理 d(n)d(n) 时发挥作用
  2. minsum(n)minsum(n) 最小质因子等比和,将在预处理 σ(n)\sigma(n) 时发挥作用 类似这样的扩展

2.2.2 因数个数函数 d(n)d(n)

  • 可以使用因数分解在 O(n)O(\sqrt n) 的时间复杂度内获得单个数 nnd(n)d(n)
  • 通过调和级数枚举在 O(nlogn)O(nlogn) 的时间获取 1n1 - nd(n)d(n) 函数值
ll d[N];
ll n = 1e6;
for(int i = 1; i <= n; ++i)
    for(int j = i; j <= n; j += i)
        d[j]++;
cpp
  • 通过欧拉筛同时预处理最小质因子指数 minexp(n)minexp(n) 可以 O(n)O(n) 得到 1n1 - nd(n)d(n) 函数值
void get_d(ll n)//欧拉筛
{
	pr.clear();
	vis.reset();
	vis[0] = vis[1] = 1;
    d[1] = 1;
	for (ll i = 2; i <= n; ++i)
	{
		if (!vis[i])
        {
            pr.push_back(i);
            d[i] = 2;
            minexp[i] = 1;
        }
		for (ll j = 0; j < pr.size() && i * pr[j] <= n; ++j)
		{
			vis[i * pr[j]] = 1;
			if (i % pr[j] == 0)
            {
                minexp[i * pr[j]] = minexp[i] + 1;
                d[i * pr[j]] = d[i] / (minexp[i] + 1) * (minexp[i * pr[j]] + 1);
                break;
            }
			d[i * pr[j]] = d[i] * 2;
            minexp[i * pr[j]] = 1;
		}
	}
}
cpp

2.2.3 因数和函数 σ(n)\sigma(n)

  • 可以使用因数分解在 O(n)O(\sqrt n) 的时间复杂度内获得单个数 nnσ(n)\sigma(n)
  • 通过调和级数枚举在 O(nlogn)O(nlogn) 的时间获取 1n1 - nσ(n)\sigma(n) 函数值
ll sigma[N];
ll n = 1e6;
for(int i = 1; i <= n; ++i)
    for(int j = i; j <= n; j += i)
        sigma[j] += i;
cpp
  • 通过欧拉筛同时预处理 minsum(n)minsum(n) 可以 O(n)O(n) 得到 1n1 - nσ(n)\sigma(n) 函数值

2.2.4 欧拉函数 ϕ(n)\phi(n)

  • 基于公式 ϕ(n)=npn(1p1)\phi(n) = n\prod_{p \mid n}(1 - p^{-1}) 可以通过分解质因数做到在 O(n)O(\sqrt n) 的时间得到单点欧拉函数值
  • 通过欧拉筛可以在 O(n)O(n) 时间复杂度预处理 1n1 - nϕ(n)\phi(n) 函数值
void get_phi(ll n)//欧拉筛
{
	pr.clear();
	vis.reset();
	vis[0] = vis[1] = 1, phi[1] = 1;
	for (ll i = 2; i <= n; ++i)
	{
		if (!vis[i])
        {
            pr.push_back(i);
            phi[i] = i - 1;
        }
		for (ll j = 0; j < pr.size() && i * pr[j] <= n; ++j)
		{
			vis[i * pr[j]] = 1;
			if (i % pr[j] == 0)
            {
                phi[i * pr[j]] = phi[i] * pr[j];
                break;
            }
			phi[i * pr[j]] = phi[i] * phi[pr[j]];
		}
	}
}
cpp

2.2.5 本质不同素因子个数函数 ω(n)\omega(n)

  • 基于质因数分解可以在 O(n)O(\sqrt n) 的时间计算单点 ω(n)\omega(n)
  • 欧拉筛可以在 O(n)O(n) 时间复杂度预处理 1n1 - nω(n)\omega(n) 函数值

2.2.6 全部素因子个数函数 Ω(n)\Omega(n)

  • 同 2.2.5 ω(n)\omega(n)

2.2.7 莫比乌斯函数 μ(n)\mu(n)

  • 质因数分解可以在 O(n)O(\sqrt n) 的时间计算单点 μ(n)\mu(n)
  • 作为积性函数可以使用欧拉筛 O(n)O(n) 预处理 1n1 - nμ(n)\mu(n) 函数值
void get_mu(ll n)
{
    vis.reset();
    pr.clear();
    vis[0] = vis[1] = 1, mu[1] = 1;
    for(ll i = 2; i <= n; ++i)
    {
        if(!vis[i])
        {
            pr.push_back(i);
            mu[i] = -1;
        }
        for(int j = 0; j < pr.size() && i * pr[j] <= n; ++j)
        {
            vis[i * pr[j]] = 1;
            if(i % pr[j] == 0)
            {
                mu[i * pr[j]] = 0;
                break;
            }
            mu[i * pr[j]] = -mu[i];
        }
    }
}
cpp

2.3 进阶武器

2.3.1 狄利克雷前后缀和与差分

狄利克雷前缀和

在已知数论函数 f(n)f(n) 的基础上计算数论函数 g(n)=dnf(d)g(n) = \sum_{d \mid n} f(d)

的各点取值,容易发现这里 g(n)g(n)f(n)f(n) 的莫比乌斯变换 但是,要注意狄利克雷前缀和是一种操作,从含义上不等同于作为数论语言的莫比乌斯变换

狄利克雷前缀和应当和普通前缀和类比着看,普通前缀和可以线性递推 O(n)O(n) 获得,那如何计算狄利克雷前缀和 形式上,狄利克雷前缀和是每个下标的因数下标的总和,那自然可以枚举每个点的因数,O(nn)O(n\sqrt n) 的时间完成 反过来想枚举每个数的倍数,容易发现我们可以调和级数枚举做到 O(nlogn)O(nlogn) 的时间复杂度

这里我们将展示一种 O(nloglogn)O(nloglogn) 的方法,这个方法基于素数筛,我们首先获取素数数组 prpr,接着枚举素数及其倍数

for(int i = 1; i <= n; ++i)
    b[i] = a[i];
for(const auto & p : pr)
    for(ll i = 1; i * p <= n; ++i)
        b[i * p] += b[i];
cpp
for(const auto & p : pr)
    for(ll i = 1; i * p <= n; ++i)
        a[i * p] += a[i];
cpp

这样的做法同时也能够把非素因子加到倍数下标上,具体的分析可以阅读下面的博客

算法学习笔记(28):狄利克雷前/后缀和、差分

狄利克雷后缀和

狄利克雷前缀和解决了因数和,后缀和则是解决倍数和的问题

在已知数论函数 f(n)f(n) 的基础上计算数论函数 g(d)=dnf(n)g(d) = \sum_{d \mid n} f(n)

类似的我们也可以做到 O(nloglogn)O(nloglogn) 的处理

for(int i = 1; i <= n; ++i)
    b[i] = a[i];
for(const auto & p : pr)
    for(ll i = n / p; i * p >= 1; ++i)
        b[i] += b[i * p];
cpp
for(const auto & p : pr)
    for(ll i = n / p; i * p >= 1; ++i)
        a[i] += a[i * p];
cpp
狄利克雷前缀差分(因数差分)

前缀差分则是狄利克雷前缀和的逆运算

在已知数论函数 g(n)g(n) 的基础上计算 f(n)f(n),且有 g(n)=dnf(d)          f(n)=dnμ(nd)g(d)g(n) = \sum_{d \mid n} f(d);;;;;f(n) = \sum_{d \mid n} \mu(\frac{n}{d})g(d)

我们也可以做到 O(nloglogn)O(nloglogn) 的处理

for(int i = 1; i <= n; ++i)
    b[i] = a[i];
for(const auto & p : pr)
    for(ll i = n / p; i * p >= 1; ++i)
        b[i * p] -= b[i];
cpp
for(const auto & p : pr)
    for(ll i = n / p; i * p >= 1; ++i)
        a[i * p] -= a[i];
cpp
狄利克雷后缀差分(倍数差分)

前缀差分则是狄利克雷后缀和的逆运算

在已知数论函数 g(n)g(n) 的基础上计算 f(n)f(n),且有 g(d)=dnf(n)g(d) = \sum_{d \mid n} f(n)

我们也可以做到 O(nloglogn)O(nloglogn) 的处理

for(int i = 1; i <= n; ++i)
    b[i] = a[i];
for(const auto & p : pr)
    for(ll i = 1; i * p <= n; ++i)
        b[i] -= b[i * p];
cpp
for(const auto & p : pr)
    for(ll i = 1; i * p <= n; ++i)
        a[i] -= a[i * p];
cpp
Dirichlet前缀和及其差分与Dirichlet卷积的关系

我们前面注意到了狄利克雷前缀和与前缀差分给出了莫比乌斯变换及其反演的一种计算方法 而莫比乌斯变换又是一种特殊的Dirichlet卷积(与单位函数做卷积) 所以利用狄利克雷卷积的关系,使用其前缀和与差分也能够预处理很多数论函数

例如

对于欧拉函数 ϕ(i)\phi(i),我们知道 ϕ1=I\phi * 1 = I

于是也可以使用狄利克雷前缀差分做到 O(nloglogn)O(nloglogn) 的预处理 ϕ(n)\phi(n)

for(int i = 1; i <= n; ++i)
    a[i] = i;
for(const auto & p : pr)
    for(int i = n / p; i * p >= 1; --i)
        a[i * p] -= a[i];
cpp

对于莫比乌斯函数 μ(n)\mu(n),我们知道 μ1=ε\mu * 1 = \varepsilon

于是也可以使用狄利克雷前缀差分做到 O(nloglogn)O(nloglogn) 的预处理 μ(n)\mu(n)

a[1] = 1;
for(const auto & p : pr)
    for(int i = n / p; i * p >= 1; --i)
        a[i * p] -= a[i];
cpp

对于因数个数函数 d(n)d(n),我们知道 11=d(n)1 * 1 = d(n)

于是也可以使用狄利克雷前缀和做到 O(nloglogn)O(nloglogn) 的预处理 d(n)d(n)

for(int i = 1; i <= n; ++i)
        a[i] = 1;
for(const auto & p : pr)
    for(int i = 1; i * p <= n; ++i)
        a[i * p] += a[i];
cpp

同理对于因数函数 σ(n)\sigma(n),我们知道 I1=σI * 1 = \sigma

也可以快速处理

你可能会认为这些常见的数论函数都可以使用欧拉筛 O(n)O(n) 直接得到了,这个没什么价值,但对于题目出现的一般数论函数,我们可能不太容易研究其素数筛怎么处理,直接使用狄利克雷前缀和或者差分是更为快捷的 例如 dndμ(d),dndϕ(d)\sum_{d \mid n}d\mu(d), \sum_{d \mid n}d\phi(d)

2.3.2 杜教筛

杜教筛给出了一类计算数论函数 f(n)f(n) 前缀和 F(n)F(n) 的方法 具体想法是寻找另一个函数 g(n)g(n)h(n)=fg(n)h(n) = f * g(n) 建立三个函数前缀和的关系

具体变形如下

H(n)=i=1nh(i)=i=1ndng(d)f(id)=d=1ng(d)di,inf(id)=d=1ng(d)i=1ndf(i)=d=1ng(d)F(nd)=d=2ng(d)F(nd)+g(1)F(n)\begin{aligned} H(n) &= \sum_{i = 1}^nh(i) \ &= \sum_{i = 1}^n\sum_{d \mid n}g(d)f(\frac{i}{d}) \ &= \sum_{d = 1}^ng(d)\sum_{d \mid i, i \leq n}f(\frac{i}{d}) \ &= \sum_{d = 1}^ng(d)\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}f(i) \ &= \sum_{d = 1}^ng(d)F(\lfloor\frac{n}{d}\rfloor) \ &= \sum_{d = 2}^ng(d)F(\lfloor\frac{n}{d}\rfloor) + g(1)F(n) \end{aligned}

最后一步将 d=1d = 1 抽取出来,就得到了计算 F(n)F(n) 的表达式,注意到这个表达式本质上是一个递推式,这意味着我们事实上需要递归的计算 F(nd)F(\lfloor\frac{n}{d}\rfloor)

在递归计算时,注意到 d=2ng(d)F(nd)\sum_{d = 2}^ng(d)F(\lfloor\frac{n}{d}\rfloor) 是可以分块计算的,所以我们同时需要计算 G(n)G(n)H(n)H(n) 的快捷方法

在假设 G(n)G(n)H(n)H(n) 可以 O(1)O(1) 计算的情况下,计算出 F(n)F(n) 的时间复杂度应该如何分析呢? 计算 F(n)F(n) 时,根据整数分块的理论,我们知道 knk \leq \sqrt nknk \geq \sqrt n 计算的 F(k)F(k) 个数相当,所以我们只考虑 knk \geq \sqrt n 的部分即可 另外对于单独的 F(k)F(k) 我们递归计算的次数大概是 O(k)O(\sqrt k),从而计算 F(n)F(n) 的时间复杂度 T(n)T(n) 可做如下估计 T(n)=i=2nO(ni)O(0nnxdx)=O(n34)T(n) = \sum_{i = 2}^{\sqrt n}O(\sqrt{\frac{n}{i}}) \sim O(\int_0^{\sqrt n}\sqrt\frac{n}{x},dx) = O(n^{\frac{3}{4}})

事实上还可以对上述过程进行优化,对于较小的取值我们可以提前使用线性筛预处理,假设对于前 knk \geq \sqrt n 的部分我们已经获得,那么需要递归计算的就只有 k<ni<nk < \lfloor\frac{n}{i}\rfloor < n,于是 2ink2 \leq i \leq \frac{n}{k} 那么此时的复杂度可做下列估计(注意不要忽略线性筛) T(n)=O(k)+i=2nkO(ni)O(0nknxdx)=O(k+nk)T(n) = O(k) + \sum_{i = 2}^{\frac{n}{k}}O(\sqrt{\frac{n}{i}}) \sim O(\int_0^{\frac{n}{k}}\sqrt\frac{n}{x},dx) = O(k + \frac{n}{\sqrt k})

观察上式我们可以得到知道当 k=n23k = n^{\frac{2}{3}} 时复杂度取最小值 O(n23)O(n^\frac{2}{3})

从而满足上述条件的数论函数前缀和可以处理到 101110^{11} 级别 下面我们列举计算 μ(n),ϕ(n),n2ϕ(n)\mu(n), \phi(n), n^2\phi(n) 的前缀和,恰好对应下面两个洛谷的题目

LG P4213 【模板】杜教筛 LG P3768 简单的数学题 (这个题还需要一些莫反变形)

f=μ,g=1,h=εf = \mu, g = 1, h = \varepsilon 可以获得计算 μ(n)\mu(n) 的前缀和 S(n)S(n) 表达式 i=1nμ(i)=1i=2nj=1niμ(j)S(n)=1i=2nS(ni)\sum_{i = 1}^n\mu(i) = 1 - \sum_{i = 2}^n\sum_{j = 1}^{\lfloor\frac{n}{i}\rfloor}\mu(j) \ S(n) = 1 - \sum_{i = 2}^nS(\lfloor\frac{n}{i}\rfloor)

f=ϕ,g=1,h=If = \phi, g = 1, h = I 可以获得计算 ϕ(n)\phi(n) 的前缀和 S(n)S(n) 表达式 i=1nϕ(i)=n(n+1)2i=2nj=1niϕ(j)S(n)=n(n+1)2i=2nS(ni)\sum_{i = 1}^n\phi(i) = \frac{n(n + 1)}{2} - \sum_{i = 2}^n\sum_{j = 1}^{\lfloor\frac{n}{i}\rfloor}\phi(j) \ S(n) = \frac{n(n + 1)}{2} - \sum_{i = 2}^nS(\lfloor\frac{n}{i}\rfloor)

f=n2ϕ,g=I2,h=I3f = n^2\phi, g = I_2, h = I_3 可以获得计算 n2ϕ(n)n^2\phi(n) 的前缀和 S(n)S(n) 表达式 i=1ni2ϕ(i)=n2(n+1)24i=2ni2j=1nij2ϕ(j)S(n)=n2(n+1)24i=2ni2S(ni)\sum_{i = 1}^ni^2\phi(i) = \frac{n^2(n + 1)^2}{4} - \sum_{i = 2}^ni^2\sum_{j = 1}^{\lfloor\frac{n}{i}\rfloor}j^2\phi(j) \ S(n) = \frac{n^2(n + 1)^2}{4} - \sum_{i = 2}^ni^2S(\lfloor\frac{n}{i}\rfloor)

2.3.3 Dirichlet双曲线法

我们希望计算数论函数 f(n)f(n) 的前缀和 F(n)F(n),我们类比杜教筛的想法,但是我们构造 f(n)=(gh)(n)f(n) = (g * h)(n)

于是 F(n)=ijng(i)h(j)=iu,ijng(i)h(j)+u<in,ijng(i)h(j)=iu,ijng(i)h(j)+vj,ijng(i)h(j)iu,jvg(i)h(j)=iug(i)jnih(j)+jvh(j)injg(i)iug(i)jvh(j)=i=1ug(i)H(ni)+j=1vh(j)G(nj)G(u)H(v)(u,vZ)\begin{aligned} F(n) &= \sum_{ij \leq n}g(i)h(j) \ &= \sum_{i \leq u, ij \leq n}g(i)h(j) + \sum_{u < i \leq n, ij \leq n}g(i)h(j) \ &= \sum_{i \leq u, ij \leq n}g(i)h(j) + \sum_{v \leq j, ij \leq n}g(i)h(j) - \sum_{i \leq u, j \leq v}g(i)h(j) \ &= \sum_{i \leq u}g(i)\sum_{j \leq \frac{n}{i}}h(j) + \sum_{j \leq v}h(j)\sum_{i \leq \frac{n}{j}}g(i) - \sum_{i \leq u}g(i)\sum_{j \leq v}h(j) \ &= \sum_{i = 1}^ug(i)H(\lfloor \frac{n}{i} \rfloor) + \sum_{j = 1}^vh(j)G(\lfloor \frac{n}{j} \rfloor) - G(u)H(v) (当 u, v \in \mathbb{Z}) \end{aligned}

注意在连续条件下应该有 uv=nuv= n,但由于整数的离散性,我们只需要在当 u,vZu, v \in \mathbb{Z} 下,u=nvu = \lfloor \frac{n}{v} \rfloor,没有这样的严格。为了平衡复杂度,通常取 u=v=nu = v = \lfloor \sqrt n \rfloor,于是

F(n)=i=1ng(i)H(ni)+j=1nh(j)G(nj)G(n)H(n)F(n) = \sum_{i = 1}^{\lfloor \sqrt n \rfloor}g(i)H(\lfloor \frac{n}{i} \rfloor) + \sum_{j = 1}^{\lfloor \sqrt n \rfloor}h(j)G(\lfloor \frac{n}{j} \rfloor) - G(\lfloor \sqrt n \rfloor)H(\lfloor \sqrt n \rfloor)

具体而言我们怎么使用这个表达式呢,首先注意到我们需要 H,GH, Gn\sqrt nnn 的取值,所以局限性是非常大的(也就是说为了计算 ff 的单点前缀和,我们需要另外两个函数同范围的大量的前缀和取值)

我们举一些例子

g=h=1g = h = 1,则 G=H=I,f=dG = H = I, f = d,于是有 i=1nd(n)=2i=1nnin2\sum_{i = 1}^nd(n) = 2\sum_{i = 1}^{\lfloor \sqrt n \rfloor}\lfloor \frac{n}{i} \rfloor - \lfloor \sqrt n \rfloor^2

g=1,h=Ig = 1, h = I,则 G=I,H(n)=n(n+1)2,f=σG = I, H(n) = \frac{n(n + 1)}{2},f = \sigma,于是有 i=1nσ(n)=i=1niH(ni)+i=1nininH(n)\sum_{i = 1}^n\sigma(n) = \sum_{i = 1}^{\lfloor \sqrt n \rfloor}iH(\lfloor\frac{n}{i}\rfloor) + \sum_{i = 1}^{\lfloor \sqrt n \rfloor}i\lfloor\frac{n}{i}\rfloor - \lfloor \sqrt n \rfloor H(\lfloor \sqrt n \rfloor)

值得注意的是以上两个公式通过拆贡献本身就有 O(n)O(\sqrt n) 的整除分块做法,所以这个公式本质无用

3. 数论分块

Liuguang_Ji的算法笔记