积性函数与筛法

【目录】

• 数论函数
• 积性函数
• 线性筛
• 狄利克雷卷积
• 杜教筛
• min_25筛

数论函数

定义

一个定义在正整数集上的实或复值函数\(f(n)\)叫做一个数论函数。

举例

• 数列\(\{a_n\}\)
• 阶乘\(n!\)
• 幂函数\(n^a\)

积性函数

定义

若数论函数\(f\)不恒等于\(0\),且当\(\gcd(m,n)=1\)时，有\(f(mn)=f(m)f(n)\)，则\(f\)叫做积性函数。若一个积性函数对于\(\forall m,n\)均有\(f(mn)=f(m)f(n)\)，则叫做完全积性函数。

性质

• 若\(f\)为积性函数，则\(f(1)=1\)
  证明： 设\(\forall n\in\mathbb{N+}\)，有\(\gcd(n,1)=1\)，故\(f(n)=f(n)f(1)\)，而\(f(n)\neq0\)，故\(f(1)=1\)
• 若\(f\)为积性函数，且\(n\)的唯一分解式为\(n=p_1^{c_1}p_2^{c_2}\cdots p_k^{c_k}\)，则\(f(n)=f(p_1^{c_1})f(p_2^{c_2})\cdots f(p_k^{c_k})\)
  这个就不证明了，很显然的。。。
• 若\(f,g\)为积性函数，则其狄利克雷卷积\(f*g\)为积性函数。
  这个下面会有证明的。

常见的积性函数

• \(\varphi(n)\)：欧拉函数
• \(\mu(n)\)：莫比乌斯函数
• \(\gcd(n,k)(k为定值时)\)
• \(d(n)=\sum\limits_{d|n}1\)：约数个数
• \(\sigma(n)=\sum\limits_{d|n}d\)：约数和
• \(\sigma_k(n)=\sum\limits_{d|n}d^k\)：约数\(k\)次幂和

常见的完全积性函数

• \(I(n)=1\)：常函数
• \(id(n)=n\)：单位函数
• \(id_k(n)=n^k\)：幂函数
• \(\epsilon(n)=[n=1]\)：元函数

线性筛

筛素数

最初，线性筛只是用来筛质数罢了。。。

void sieve(int n) {
  static int v[N], p[N], pr;
  // v[i] 表示 i 的最小质因子
  // p[N] 和 pr 用来存质数表
  for (int i = 2; i <= n; ++i) {
    if (v[i] == 0) v[i] = i, p[++pr] = i;
    // 没被筛，是质数
    for (int j = 1; j <= pr && i*p[j] <= n; ++j) {
      v[i*p[j]] = p[j];
      if (i % p[j] == 0) break;
      // 当 i 是 p[j] 的倍数时 break
    }
  }
}

代码简短，便于记忆。。。
然而背代码是不可能的，这辈子不可能背代码的。。。

正确性证明

• 在筛的过程中，上面的程序将\(i \cdot p[j]\)的最小质因子定为\(p[j]\)。
• 先不考虑其正确性，将\(i\)唯一分解，得到\(i=p_1^{c_1}\cdots p_k^{c_k}\)。
• 而此时一定有\(p[j]\leq p_1\)，因为当\(p[j]=p_1\)时就已经break了。
• 同样由于\(p[j]\leq p_1\)，\(i \cdot p[j]\)的最小质因子显然为\(p[j]\)。

线性时间复杂度证明

• 考虑每个数\(x\)被筛到的原因。
• 若\(x\)为质数，则\(x\)被筛是理所当然的。
• 若\(x\)非质数，则\(x\)只会被其最小质因子筛到。
• 综上述，每个数会且仅会被筛一次，故时间复杂度为\(O(n)\)。

筛积性函数

首先要完全理解线性筛质数的过程。
对于积性函数\(f(n)\)，首先可以断定\(f(1)=1\)；
因为对于\(\forall n \in \mathbb{N_+}\)，显然有\(\gcd(n,1)=1\)，所以\(f(n)=f(n)f(1)\)，故\(f(1)=1\)。
注：以下用\(\mathbb{P}\)表示质数集
若想线性筛出积性函数\(f(n)\)，就必须能快速求出下面两个值

• \(f(p)(p\in\mathbb{P})\)
• \(f(p^k)(p\in\mathbb{P},k\in\mathbb{N_+})\)

重新思考线性筛的过程

• 将\(i\)唯一分解，得到\(i=p_1^{c_1}\cdots p_k^{c_k}\)。
• 而此时一定有\(p[j]\leq p_1\)，因为当\(p[j]=p_1\)时就已经break了。
• 同样由于\(p[j]\leq p_1\)，\(i \cdot p[j]\)的最小质因子显然为\(p[j]\)。

考虑将其推广到求积性函数

• 若\(p[j]<p_1\)，则\(\gcd(i,p[j])=1\)，于是

\[f(i\cdot p[j])=f(i)f(p[j]) \]

• 若\(p[i]=p_1\)，则需要记录一个\(low_i=p_1^{c_1}\)，即\(i\)唯一分解式中最小质因子的幂。 可以发现\(\dfrac{i}{low_i}\)的最小质因子一定大于\(p_1\)，而\(p[j]=p_1\)，故

\[\gcd\left(\dfrac{i}{low_i},low_i\cdot p[j]\right)=1 \]

所以

\[f(i\cdot p[j])=f\left(\dfrac{i}{low_i}\right)f(low_i\cdot p[j]) \]

• 注意当\(low_i=i\)，即\(i=p^k(p\in\mathbb{P},k\in\mathbb{N_+})\)时，需要根据定义直接计算\(f(i\cdot p[j])\)

代码如下

void sieve(int n) {
  static int v[N], p[N], low[N], pr;
  f[1] = 1;
  for (int i = 2; i <= n; ++i) {
    if (!v[i]) {
      // i 为质数
      v[i] = i; // 最小质因子
      low[i] = i; // 唯一分解式中最小质因子的幂
      p[++pr] = i; // 加入质数表
      f[i] = ...; // 由定义直接计算
    }
    for (int j = 1; j <= pr && i*p[j] <= n; ++j) {
      v[i*p[j]] = p[j];
      if (i % p[j] == 0) {
        low[i*p[j]] = low[i] * p[j]; // 递推计算 low
        if (low[i] == i)
          f[i*p[j]] = ...; // 由定义直接计算
        else
          f[i*p[j]] = f[i/low[i]] * f[low[i]*p[j]];
        break;
      } else {
        low[i*p[j]] = p[j]; // 这是显然的
        f[i*p[j]] = f[i] * f[p[j]];
      }
    }
  }
}

狄利克雷卷积

我们知道

\[\varphi(n)=\sum\limits_{d|n}\mu(d)\dfrac{n}{d} \]

其右端和的形式，在数论中会经常出现。于是乎。。。

定义

设\(f,g\)是两个数论函数，则他们的狄利克雷卷积\(h(n)\)也是一个数论函数，由下式给出

\[h(n)=\sum\limits_{d|n}f(d)g\left(\dfrac{n}{d}\right) \]

简记为\(h=f*g\)。

性质

狄利克雷卷积有如下性质

• 交换律：\(f*g=g*f\)，这个结论还是很显然的。
• 结合律：\((f*g)*h=f*(g*h)\)。
  这是因为

\[\sum\limits_{(i*j)*k=n}(f(i)g(j))h(k)=\sum\limits_{i*(j*k)=n}f(i)(g(j)h(k)) \]

• \(\epsilon*f=f\)
• 对于满足\(f(1)=\not0\)的数论函数\(f\)，存在唯一的数论函数\(g\)使得\(f*g=\epsilon\)，\(g\)称为\(f\)的狄利克雷逆函数。
  逆函数的计算
  直接构造

\[g(n)=\dfrac{1}{f(1)}\left([n=1]-\sum\limits_{d|n,d\neq1f(d)}g\left(\dfrac{n}{d}\right)\right) \]

这样

\[\begin{aligned} &\sum\limits_{d|n}f(d)g\left(\dfrac{n}{d}\right) \\ =&f(1)g(n) + \sum\limits_{d|n,d\neq1}f(d)g\left(\dfrac{n}{d}\right)\\ =&[n=1] \end{aligned}\]

常用的狄利克雷卷积公式

• \(\mu*I=\epsilon\)
• \(\varphi*I=id\)
• \(\mu*id=\varphi\)

杜教筛

杜教筛可以通过构造狄利克雷卷积在非线性时间内求积性函数的前缀和。

前置知识

• 数论分块
• 积性函数
• 莫比乌斯反演
• 狄利克雷卷积

构造方法

给定积性函数\(f,g\)，设\(S(n)=\sum\limits_{i=1}^nf(i)\)，则有

\[\begin{aligned} \sum\limits_{i=1}^n(f*g)(i)&=\sum\limits_{i=1}^n\sum\limits_{d|i}f(d)g\left(\frac{i}{d}\right) \\ &=\sum\limits_{d=1}^ng(d)\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i) \\ &=\sum\limits_{d=1}^ng(d)S\left(\lfloor\frac{n}{d}\rfloor\right) \end{aligned}\]

由于\(g(1)=1\)，可以构造

\[S(n)=g(1)S(n)=\sum\limits_{d=1}^ng(d)S\left(\lfloor\frac{n}{d}\rfloor\right)-\sum\limits_{d=2}^ng(d)S\left(\lfloor\frac{n}{d}\rfloor\right) \]

而

\[\sum\limits_{d=1}^ng(d)S\left(\lfloor\frac{n}{d}\rfloor\right)=\sum\limits_{i=1}^n(f*g)(i) \]

故

\[S(n)=\sum\limits_{i=1}^n(f*g)(i)-\sum\limits_{d=2}^ng(d)S\left(\lfloor\frac{n}{d}\rfloor\right) \]

如果可以构造出合适的积性函数\(g\)，使\(\sum\limits_{i=1}^n(f*g)(i)\)可以快速计算，就可以使用数论分块递归地求出\(S(n)\)。
代码大概长这个样子

int S(int n) {
  int ret = ...; // f*g 的前缀和
  for (int l = 2, r; l <= n; l = r+1) {
    // 数论分块, 注意 l 从 2 开始
    r = n / (n/l);
    ret -= (sum_g[r]-sum_g[l-1]) * S(n/l);
    // sum_g 表示 g 的前缀和
  }
  return ret;
}

这样的时间复杂度为

\[T(n)=\sum\limits_{i=1}^{\sqrt n}O(\sqrt i)+O(\sqrt \frac{n}{i})=O(n^{\frac{3}{4}}) \]

时间复杂度的优化

线性筛与杜教筛结合使用。先用线性筛筛出前\(m\)个数的答案，之后再用杜教筛。
这样时间复杂度为

\[T(n)=\sum\limits_{i=1}^{\lfloor\frac{n}{m}\rfloor}\sqrt \frac{n}{i}=O\left(\frac{n}{\sqrt m}\right) \]

当\(m=n^{\frac{2}{3}}\)时，\(T(n)\)取最小值\(O(n^{\frac{2}{3}})\)
至于记忆化，可以手写哈希表，也可以用 map 或 unordered_map。
代码如下

unordered_map<int, int> rec;
// 这里使用 C++11 的 unordered_map
int S(n) {
  if (n <= m) return sum[n];
  // sum 为线性筛预处理的前缀和
  if (rec[n]) return rec[n];
  int ret = &rec[n];
  ret = ...; // f*g 的前缀和
  for (int l = 2, r; l <= n; l = r+1) {
    r = n / (n/l);
    ret -= (sum_g[r]-sum_g[l-1]) * S(n/l);
    // sum_g 表示 g 的前缀和
  }
  return ret;
}

构造举例

• \(S(n)=\sum\limits_{i=1}^n\varphi(i)\)
  因为\(\varphi*I=id\)，所以

\[\begin{aligned} S(n)=I(1)S(n)&=\sum\limits_{i=1}^n(\varphi*I)(i)-\sum\limits_{d=2}^nI(d)S\left(\lfloor\frac{n}{d}\rfloor\right)\\ &=\sum\limits_{i=1}^nid(i)-\sum\limits_{d=2}^nS\left(\lfloor\frac{n}{d}\rfloor\right)\\ &=\dfrac{n(n+1)}{2}-\sum\limits_{d=2}^nS\left(\lfloor\frac{n}{d}\rfloor\right) \end{aligned}\]

代码如下

unordered_map<int, int> rec;

int S(n) {
  if (n <= m) return sum[n];
  if (rec[n]) return rec[n];
  int ret = &rec[n];
  ret = n * (n+1) / 2;
  for (int l = 2, r; l <= n; l = r+1) {
    r = n / (n/l);
    ret -= (r-l+1) * S(n/l);
  }
  return ret;
}

• \(S(n)=\sum\limits_{i=1}^n\mu(i)\)
  因为\(\mu*I=\epsilon\)，所以

\[\begin{aligned} S(n)=I(1)S(n)&=\sum\limits_{i=1}^n(\mu*I)(i)-\sum\limits_{d=2}^nI(d)S\left(\lfloor\frac{n}{d}\rfloor\right)\\ &=\sum\limits_{i=1}^n\epsilon(i)-\sum\limits_{d=2}^nS\left(\lfloor\frac{n}{d}\rfloor\right)\\ &=1-\sum\limits_{d=2}^nS\left(\lfloor\frac{n}{d}\rfloor\right) \end{aligned}\]

代码如下

unordered_map<int, int> rec;

int S(n) {
  if (n <= m) return sum[n];
  if (rec[n]) return rec[n];
  int ret = &rec[n];
  ret = 1;
  for (int l = 2, r; l <= n; l = r+1) {
    r = n / (n/l);
    ret -= (r-l+1) * S(n/l);
  }
  return ret;
}

• \(S(n)=\sum\limits_{i=1}^n i\cdot\varphi(i)\)
  设\(f(n)=n\cdot\varphi(n)\)，即\(S(n)=\sum\limits_{i=1}^nf(i)\)
  构造\(g(n)=n\)，即\(g=id\)
  因为

\[\begin{aligned} (f*g)(n)&=\sum\limits_{d|n}f(d)g\left(\frac{n}{d}\right)\\ &=\sum\limits_{d|n}d\cdot\varphi(d)\cdot\dfrac{n}{d}\\ &=\sum\limits_{d|n}n\cdot\varphi(d)\\ &=n\sum\limits_{d|n}\varphi(d)\\ &=n^2 \end{aligned}\]

所以

\[\sum\limits_{i=1}^n(f*g)(i)=\sum\limits_{i=1}^ni^2=\dfrac{n(n+1)(2n+1)}{6} \]

故

\[\begin{aligned} S(n)=g(1)S(n)&=\sum\limits_{i=1}^n(f*g)(i)-\sum\limits_{d=2}^ng(d)S\left(\lfloor\frac{n}{d}\rfloor\right)\\ &=\dfrac{n(n+1)(2n+1)}{6}-\sum\limits_{d=2}^nd\cdot S\left(\lfloor\frac{n}{d}\rfloor\right) \end{aligned}\]

代码如下

unordered_map<int, int> rec;

int S(n) {
  if (n <= m) return sum[n];
  if (rec[n]) return rec[n];
  int ret = &rec[n];
  ret = n * (n+1) * (2*n+1) / 6;
  for (int l = 2, r; l <= n; l = r+1) {
    r = n / (n/l);
    ret -= (l+r) * (r-l+1) / 2 * S(n/l);
  }
  return ret;
}

min_25筛

留个坑吧，应该会填的。。。