跳转至

狄利克雷卷积

本文介绍 Dirichlet 卷积和 Dirichlet 生成函数.

Dirichlet 卷积

数论函数 f(n)g(n)Dirichlet 卷积(Dirichlet convolution),记作 fg,定义为数论函数

(fg)(n)=knf(k)g(nk)=k=nf(k)g().

Dirichlet 卷积是数论函数的重要运算.数论函数的许多性质都是通过这个运算挖掘出来的.

例子
  1. 单位函数 ε 是莫比乌斯函数 μ 和常数函数 1 的 Dirichlet 卷积:
ε=μ1ε(n)=dnμ(d).
  1. 除数个数函数 τ 是常数函数 1 和它自身的 Dirichlet 卷积:
τ=11τ(n)=dn1.
  1. 除数和函数 σ 是恒等函数 id 和常数函数 1 的 Dirichlet 卷积:
σ=id1σ(n)=dnd.
  1. 欧拉函数 φ 是恒等函数 id 和莫比乌斯函数 μ 的 Dirichlet 卷积:
φ=idμφ(n)=dndμ(nd).

莫比乌斯反演 就是利用 ε=μ1 对数论函数恒等式进行变形.

性质

Dirichlet 卷积具有一系列代数性质.

定理

f,g,h 都是数论函数.那么,有:

  1. 交换律fg=gf
  2. 结合律(fg)h=f(gh)
  3. 分配律(f+g)h=fh+gh
  4. 单位元fε=εf=f,其中,ε(n)=[n=1] 是卷积单位元,[] 是 Iverson 括号.
  5. 逆元:当且仅当 f(1)0 时,存在 g 使得 fg=gf=ε,且 g 称为 fDirichlet 逆元(Dirichlet inverse),可以记作 f1.而且,逆元 g 满足递推公式
g(n)=ε(n)k=n, k1f(k)g()f(1).
证明

为验证交换律,直接计算可知

(fg)(n)=k=nf(k)g()=(gf)(n).

为验证结合律,直接计算可知

((fg)h)(n)=km=nf(k)g()h(m)=(f(gh))(n).

为验证分配律,直接计算可知

((f+g)h)(n)=k=n(f(k)+g(k))h()=k=nf(k)h()+k=ng(k)h()=(fh+gh)(n).

为验证 ε(n) 是单位元,直接计算可知

(fε)(n)=k=nf(k)ε()=f(n).

第二个等号是因为 ε() 仅在 =1k=n 时取得非零值.

最后,需要证明 f1 存在,当且仅当 f(1)0.对于任意 f,假设存在 g 使得 fg=ε.这意味着

(fg)(n)=k=nf(k)g()=ε(n).

这实际上给出了一系列关于 g(n) 取值的方程组,从中可以直接求出 g(n).特别地,当 n=1 时,等式变为 f(1)g(1)=1,所以 g 存在,至少要求 f(1)0.而只要 f(1)0,可以直接解出

g(n)=ε(n)k=n, k1f(k)g()f(1).

它可以用于递归计算 g(n) 的取值.因此,逆元 g 存在,当且仅当 f(1)0

用抽象代数的语言说,这些代数性质说明,全体数论函数在(逐点)加法运算和 Dirichlet 卷积运算下构成 交换环,且它的全体可逆元就是那些在 n=1 处取非零值的函数.这个环称为 Dirichlet 环(Dirichlet ring).

积性函数是一类特殊的数论函数.它对于 Dirichlet 卷积和 Dirichlet 逆都是封闭的.

定理

f,g 是积性函数,那么,fg 也是积性函数.而且,逆元 f1 一定存在,它也是积性函数.

证明

对于第一点,设 h=fg,直接验证可知,对于 n1n2,都有

h(n1)h(n2)=(k11=n1f(k1)g(1))(k22=n2f(k2)g(2))=k11=n1, k22=n2f(k1)f(k2)g(1)g(2)=k=n1n2f(k)g()=h(n1n2).

其中,第三个等号改变求和顺序的逻辑是:当 k 遍历 n1n2 的因数时,k 的素因子可以根据它是 n1 还是 n2 的素因子分为两类,将两类中的素因子(计重复)分别乘起来得到 k1k2,它们将分别遍历 n1n2 的因数;反过来,根据 n1n2 的因数 k1k2,总是可以得到 n1n2 的因数 k=k1k2

对于第二点,设 g=f1,考虑应用数学归纳法.首先,g(1)=1/f(1)=1.此时,逆元的递归公式可以写作

g(n)=ε(n)k=n, k1f(k)g().

所以,对于 n1n2n1n2>1,有

g(n1n2)=k=n1n2, k1f(k)g()=k11=n1, k22=n2, k1k21f(k1)f(k2)g(1)g(2)=f(1)f(1)g(n1)g(n2)k11=n1, k22=n2f(k1)f(k2)g(1)g(2)=g(n1)g(n2)(k11=n1f(k1)g(1))(k22=n2f(k2)g(2))=g(n1)g(n2)ε(n1)ε(n2)=g(n1)g(n2).

其中,第二个等号用到了归纳假设,即对于 12<n1n212,条件 g(12)=g(1)g(2) 成立.

用抽象代数的语言说,全体积性函数在 Dirichlet 卷积运算下构成 Dirichlet 环乘法群的 子群

更为特殊的是完全积性函数.

定理

α 是完全积性函数,f,g 是数论函数.那么,有:

  1. 分配律:(αf)(αg)=α(fg)
  2. 逆元:(αf)1=αf1,只要 f1 存在.
  3. 积性函数 f 是完全积性函数,当且仅当 f1=μf,其中,μ莫比乌斯函数
证明

对于第一条,直接验证可知

((αf)(αg))(n)=k=n(αf)(k)(αg)()=k=nα(k)f(k)α()g()=k=nα(n)f(k)g()=α(n)(fg)(n).

其中,第三个等号用到了完全积性函数的性质:α(n)=α(k)α() 对所有 n=k 都成立.

对于第二条,利用第一条就有

(αf)(αf1)=α(ff1)=αε=ε.

其中,最后一个等号只利用了 α(1)=1.由逆元定义,(αf)1=αf1

对于第三条,利用第二条和 11=μ 可知,如果 f 是完全积性函数,那么

f1=(1f)1=11f=μf.

其中,1 是常数函数.反过来,如果 f 是积性函数且 f1=μf,那么只需要证明对于所有素数 peN+,都有 f(pe)=f(p)e 成立,就能证明 f 是完全积性函数.为此,对 eN+ 应用数学归纳法.归纳起点 e=1 处命题显然成立.对于任意 e>1,应用逆元递推公式,都有

f1(pe)=i=1ef(pi)f1(pei)=i=1ef(pi)μ(pei)f(pei)=f(pe)f(1)μ(1)f(pe1)μ(p)f(p)=f(pe)+f(p)e.

其中,最后一个等号用到了归纳假设 f(pe1)=f(p)e1.应用 f1=μf,就得到

f1(pe)=μ(pe)f(pe)=0.

代入前式,就得到

f(pe)=f(p)e.

所以,归纳步骤成立.原命题得证.

用抽象代数的语言说,如果 α 是完全积性函数,映射 fαf 是 Dirichlet 环的 自同态

Dirichlet 生成函数

与 Dirichlet 卷积紧密相关的是 Dirichlet 生成函数.

数论函数 f(n)——也就是数列 {f(n)}——对应的 Dirichlet 生成函数(Dirichlet series generating function,DGF)定义为形式 Dirichlet 级数(formal Dirichlet series):

F(s)=n=1f(n)ns.

级数中的 s 是形式变元.常见的 Dirichlet 生成函数中,s 往往可以看作是复变量,进而讨论 Dirichlet 级数的解析性质,但这超出了算法竞赛的范围.

Dirichlet 生成函数的乘积对应着相应的数论函数的 Dirichlet 卷积:

定理

对于数论函数 f,g 及其 Dirichlet 生成函数 F,G,它们的 Dirichlet 卷积 fg 的生成函数等于 FG

证明

直接验证:

F(s)G(s)=(k=1f(k)ks)(=1g()s)=k=1=1f(k)g()(k)s=n=1k=nf(k)g()ns=n=1(fg)(n)ns.

利用 Dirichlet 卷积和 Dirichlet 生成函数乘积之间的对应关系,可以从 Dirichlet 生成函数的角度理解 Dirichlet 卷积的性质.由于形式 Dirichlet 级数的乘法运算满足交换律、结合律、对加法的分配律,数论函数的 Dirichlet 卷积运算满足同样的代数性质.

Euler 乘积

积性函数的特殊性同样反映在 Dirichlet 生成函数上.由于整数有 唯一分解定理,积性函数 f(n) 的生成函数 F(s) 可写成如下形式:

F(s)=n=1f(n)ns=n=1pPf(pe)pes=pPe=0f(pe)pes=pP(1+f(p)ps+f(p2)p2s+f(p3)p3s+).

这意味着,F(s) 可以分解为若干 Fp(s) 的乘积,且每个 Fp(s) 对应的数论函数都只在 p 的幂次处可能取非零值.这一无穷乘积也称为 Euler 乘积(Euler product).如果 F(s)G(s) 都能分解成类似的形式,那么它们的乘积同样如此;将这一观察对应到数论函数上,就是积性函数的 Dirichlet 卷积仍然是积性函数.

进一步地,如果 f(n) 还是完全积性函数,那么 f(pe)=f(p)e,上式可以继续简化:

F(s)=pPe=0f(p)epes=pP(1f(p)ps)1.

与积性函数不同,完全积性函数的 Dirichlet 生成函数的形式在乘法运算下并不具有封闭性,因此,完全积性函数的 Dirichlet 卷积和 Dirichlet 逆都未必是完全积性函数,但一定是积性函数.

例子
  1. 单位函数 ε(n) 是完全积性函数.它的 Dirichlet 生成函数是关于不定元 s 的常值函数

    E(s)=n=1ε(n)ns=1.
  2. 常数函数 1(n) 是完全积性函数.它的 Dirichlet 生成函数是 Riemann 函数

    I(s)=n=11ns=pP11ps=ζ(s).
  3. 莫比乌斯函数 μ(n) 是常数函数的 Dirichlet 逆.它的 Dirichlet 生成函数是 ζ(s) 的倒数:

    M(s)=n=1μ(n)ns=pP(1ps)=1ζ(s).
  4. 幂函数 idk(n)=nk 是完全积性函数.特别地,当 k=0 时,它就是常数函数;当 k=1 时,它就是恒等函数.它的 Dirichlet 生成函数是

    Ik(s)=n=1nkns=pP11pks=ζ(sk).
  5. 欧拉函数 φ(n) 是积性函数.它的 Dirichlet 生成函数是

    Φ(s)=pP(1+p1ps+p(p1)p2s+p2(p1)p3s+)=pP(11p1s1ps11p1s)=pP1ps1p1s=ζ(s1)ζ(s).

    结合幂函数的 Dirichlet 函数表达式,就得到 id=φ1

  6. 约数函数 σk(n)=dndk 是积性函数.它的 Dirichlet 生成函数是

    Σk(s)=pP(1+1+pkps+1+pk+p2kp2s+1+pk+p2k+p3kp3s+)=pP11pk((1pk)+1p2kps+1p3kp2s+1p4kp3k+)=pP11pk(11pspk1pks)=pP1(1ps)(1pks)=ζ(sk)ζ(s).

    结合幂函数的 Dirichlet 表达式,就得到 σk=idk1.这正是 σk 的定义式.

  7. 无平方因子数的指示函数 u(n)=|μ(n)| 是积性函数.它的 Dirichlet 生成函数是

    U(s)=pP(1+ps)=pP1p2s1ps=ζ(s)ζ(2s).

应用

Dirichlet 生成函数可以用于将积性函数表示为 Dirichlet 卷积.

例如在杜教筛的过程中,要计算积性函数 f 的前缀和,需要找到另一个积性函数 g 使得 fgg 都可以快速求前缀和.可以利用 Dirichlet 生成函数推导这一过程.

以杜教筛一节的例题 Luogu P3768 简单的数学题 为例,需要对 f(n)=n2φ(n) 构造满足上述条件的数论函数 g(n).由于 f 是积性函数,它的 Dirichlet 生成函数为

F(s)=pP(1+k=1p3k1(p1)pks)=pP1p2s1p3s=ζ(s3)ζ(s2).

对比幂函数的 Dirichlet 生成函数可知,只要取 g=id2,就有 fg=id3.两者都是可以快速计算前缀和的.

Dirichlet 卷积的计算

本节讨论 Dirichlet 卷积的计算问题,即给定序列 {f(k)}k=1n{g(k)}k=1n,求解 Dirichlet 卷积 h=fg 的前若干项 {h(k)}k=1n 的问题.根据涉及到的函数性质,算法的复杂度也略有不同.

一般情形

如果 f,g,h 都没有特殊性质,那么 Dirichlet 卷积的计算只能利用其定义:

h(n)=k=nf(k)g().

枚举 k,将贡献 f(k)g() 累加到 h(k) 上即可.枚举复杂度为

O(k=1nnk)=O(nlogn).

参考实现如下:

参考实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
// Compute the Dirichlet convolution h = f * g.
auto dirichlet_convolute(const std::vector<int>& f, const std::vector<int>& g) {
  int n = f.size() - 1;
  std::vector<int> h(n + 1);
  for (int k = 1; k <= n; ++k) {
    for (int d = 1; k * d <= n; ++d) {
      h[k * d] += f[k] * g[d];
    }
  }
  return h;
}

与积性函数卷积的情形

如果 g 是积性函数,那么可以利用 Euler 乘积加速 Dirichlet 卷积的计算.计算 h 相当于计算它的 Dirichlet 生成函数 H 中各项的系数.由于

H(s)=F(s)G(s)=F(s)pPGp(s).

其中,Gp(s)G(s) 的 Euler 乘积分解中的因式,它只包含 p 的幂次处的系数:

Gp(s)=pknf(pk)pks=1+f(p)ps+f(p2)p2s+.

那么,从 F(s) 开始,遍历所有不超过 n 的素数 p,将 Gp(s) 逐一乘上去,同样可以得到最终结果 H(s).将 Gp(s) 乘上去时,直接应用一般情形中的暴力枚举算法即可.总枚举次数

pP, pnk=1npkpP, pnnp1pP, pn2npO(nloglogn).

最后一步复杂度的估计与 Eratosthenes 筛法 复杂度的证明一致.所以,本算法的时间复杂度为 O(nloglogn)

参考实现如下:

参考实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
// Compute the Dirichlet convolution h = f * g.
// Assume that g is multiplicative.
auto dirichlet_convolute(const std::vector<int>& f, const std::vector<int>& g) {
  int n = f.size() - 1;
  std::vector<int> h(f);
  std::vector<bool> vis(n + 1);
  for (int i = 2; i <= n; ++i) {
    if (vis[i]) continue;
    // Reverse the order for in-place computation.
    for (int k = n / i * i; k; k -= i) {
      vis[k] = true;
      int d = k;
      while (true) {
        d /= i;
        h[k] += h[d] * g[k / d];
        if (d % i) break;
      }
    }
  }
  return h;
}

特别地,当积性函数 g 是完全积性函数或其 Dirichlet 逆时,例如当 g=1g=μ 时,那么算法可以进一步简化。此时,Dirichlet 卷积 h=fg 的计算可以采用常数更小的 Dirichlet 前缀和/差分 算法,但是算法时间复杂度仍为 O(nloglogn)

结果为积性函数的情形

最后,考虑 h 是积性函数的情形.特别地,当 f,g 都是积性函数时,h=fg 就是积性函数.要计算 h,只需要确定它在素数幂处的取值,就可以通过 线性筛O(n) 时间内计算.而对于素数幂 pe 处的取值 h(pe) 直接暴力计算即可:

h(pe)=i=0ef(pi)g(pei).

这些暴力计算需要的枚举次数为

pP, pne=1logpn(e+1)pP, pnlogpn2+pP, n<pn1n(log2n)2+nO(n).

因此,这一算法的总时间复杂度为 O(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
// Compute the Dirichlet convolution h = f * g.
// Assume that h is multiplicative.
auto dirichlet_convolute(const std::vector<int>& f, const std::vector<int>& g) {
  int n = f.size() - 1;
  std::vector<int> h(n + 1), primes, rem(n + 1), lpf(n + 1);
  std::vector<bool> vis(n + 1);
  h[1] = 1;
  for (int x = 2; x <= n; ++x) {
    if (!vis[x]) {
      primes.push_back(x);
      rem[x] = 1;
      lpf[x] = x;
    }
    for (int p : primes) {
      if (x * p > n) break;
      vis[x * p] = true;
      rem[x * p] = x % p ? x : rem[x];
      lpf[x * p] = p;
      if (x % p == 0) break;
    }
    if (rem[x] == 1) {  // prime powers.
      for (int k = x; k; k /= lpf[x]) {
        h[x] += f[k] * g[x / k];
      }
    } else {  // other cases.
      h[x] = h[rem[x]] * h[x / rem[x]];
    }
  }
  return h;
}

参考资料与注释