【数学】各种积性函数的线性筛法

前置芝士：几种特殊的积性函数的定义及基本性质。

【资料图】

定义

积性函数：若函数\(f(x)\)满足\(f(x) = 1\)且\(\forall x,y \in N^+,gcd(x,y) = 1\) ，都有\(f(xy) = f(x)f(y)\)，则\(f(x)\)为积性函数。

完全积性函数：若函数满足\(f(x) = 1\)且\(\forall x,y \in N^+\)，都有\(f(xy) = f(x)f(y)\)，则\(f(x)\)为积性函数。

判定

特殊例子：（以下都是积性函数）

\(\epsilon(n) = [n = 1]\) （\(n = 1\)时为1，否则为0）

\(id_k(n) = n^k\)，当\(k = 1\)时简记为\(id(n)\)

\(1(n) = 1\)

\(\sigma_k(n) = \sum_{d|n}d^k\)，当\(k = 0\)时为因数个数，记为\(d(n)\)，\(k = 1\)时为因数和，记为\(\sigma(n)\)

\(\phi(n) = \sum_{i = 1}^n [gcd(i,n) = 1]\)

\[\mu(n) = \begin{cases}1 \ \ \ \ n = 1\\(-1)^k \ \ \ \ 所有质因数次数不超过1且有k种质因数\\0 \ \ \ \ otherwise\\\end{cases}\]

对于函数之间运算产生的函数，有以下规律：

若\(f(x)\)和\(g(x)\)都是积性函数，则以下函数也是积性函数：

\[h(x) = f(x^p)\\h(x) = f^p(x)\\h(x) = f(x)g(x)\\f(x)和g(x)的Dirichlet卷积，即：h(x) = \sum_{d|x}f(x)g(\frac gd)\]

有时在题目中，会要求筛出\(1 \to 10^6\)甚至\(1 \to 10^7\)的这些函数值，时间复杂度要求\(O(n)\)，我们如何计算这些东西呢？

我们发现一个普遍的性质，这些函数当自变量取值为质数时，都可以简单地\(O(1)\)求出，于是我们将这个过程代入筛素数的线性筛，再尝试每次用质数（创造互质条件）和积性函数的性质求出函数值。

不会线性筛质数的自行前往：P3383 【模板】线性筛素数 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

具体对于每一个函数将怎么做呢？

欧拉函数\(\phi\)

直接枚举比\(x\)小的数显然不可行，我们知道欧拉函数有一种计算方法：

\[\phi(x) = x\prod_{p \in prime\ \& \ p|x}\frac {p - 1}p\]

因为无论一个质数\(p\)在\(x\)的分解中出现多少次，它\(\frac {p - 1}p\)的贡献都只有一次，结合到质数筛中我们每次都用一个质数去乘以一个其他数，可以分类讨论：

\[\phi(x) = x - 1\ \ \ if \ x \in prime\]

以及(\(p\)是线性筛中每次枚举的质数)

\[\phi(x * p) = \begin{cases}\phi(x) * \phi(p)\ \ \ \ p \mid x\\\phi(x) * p\ \ \ \ p \nmid x\end{cases}\]

因为\(p \mid x\)时，满足\(gcd(p,x) = 1\)，所以可以直接相乘。

\(p \nmid x\)时，\(x\)中已有\(p\)，贡献就只有乘上了\(p\)倍。

由于线性筛一定会把每个数筛一次，所以保证范围内的数全都计算了\(\phi\)值。

（代码中的\(\phi(p)\)写作了\(p - 1\)，因为\(p\)是质数，两者是等价的，后面的函数也一样）

Code

inline void init(){    phi[1] = 1;    for(R int i = 2;i <= N - 1;i++)    {        if(!vis[i])        {            phi[i] = i - 1;            prime[++tot] = i;        }        for(R int j = 1;j <= tot && i * prime[j] < N;j++)        {            vis[prime[j] * i] = 1;            if(!(i % prime[j]))            {                phi[prime[j] * i] = (prime[j] - 1) * phi[i];                break;            }            phi[prime[j] * i] = phi[i] * prime[j];        }    }}

莫比乌斯函数\(\mu\)

同\(\phi\)的分析方法，我们不难发现，当自变量是一个质数的时候，其一定只有一个质因数。

即

\[\mu(x) = -1\ \ \ if\ x \in prime\]

考虑到线性筛中一个任意数和一个质数相乘的情况，如果\(p \mid i\)，那么\(p * i\)一定含有超过一个\(p\)，\(\mu\)为0，如果\(p \nmid i\)，至少说明乘上\(p\)这一操作让\(i\)多了一个质因数，不管之前\(\mu\)是否为0，将\(\mu\)取反，一定满足其定义。

\[\mu(i * p) = \begin{cases}\mu(i) * \mu(p) = -\mu(i)\ \ \ \ p \nmid i\\0\ \ \ \ p \mid i\\\end{cases}\]

Code

inline void init(){    miu[1] = 1;    for(R int i = 2;i <= N - 1;i++)    {        if(!vis[i])        {            miu[i] = -1;            prime[++tot] = i;        }        for(R int j = 1;j <= tot && i * prime[j] < N;j++)        {            vis[prime[j] * i] = 1;            if(!(i % prime[j]))            {                miu[prime[j] * i] = 0;                break;            }            miu[prime[j] * i] = -miu[i];        }    }}

因数和\(\sigma\)

精彩的部分来了，本人第一次看到时也十分震撼。

方法源自ldysy2102 - 博客园的博客线性筛约数个数、约数和的新方法 - ldysy2102 - 博客园 (cnblogs.com)

首先当\(p\)为质数时，因数只有\(1、p\)两个，因数和\(\sigma(p) = p + 1\)

由唯一分解定理，一个数可以被分解为：

\[x = \prod_{i = 1\ \& \ p_i \in prime}^kp_i^{c_i}\]

这个数的因数就是这些质因数任意次数的任意组合，每个质因数都可以选择\(0 \to c_i\)个。所以

每个因数就是质因数的任意次方乘起来，是：

\[\sigma(x) = (1 + p_1 + p_1^2 + .. + p_1^{c_1})(1 + p_2 + p_2^2 + .. + p_2^{c_2})...(1 + p_k + .. + p_k^{c_k})\]

令\((1 + p_2 + ..)..(1 + p_k + ..) = T\)，再次讨论线性筛中一个质数乘上一个任意数的情况，在\(p \mid i\)时，我们钦定乘上的质因数是\(p_1\)：

\[\sigma(x * p_1) = (1 + p_1 + .. + p_1^{c_1 + 1})T\\= \sigma(x) + p_1^{c_1 + 1}T\]\[\sigma(\frac x{p_1}) = (1 + p_1 + .. + p_1^{c_1 - 1})T\\= \sigma(x) - p_1^{c_1}T\]

两个柿子只有\(T\)前的系数不一样，将\(②\)式乘上\(p_1\)再加上\(①\)，神奇的事情发生了：

\[p_1\sigma(\frac x{p_1}) + \sigma(x * p_1) = (p_1 + 1)\sigma(x)\]\[\sigma(x * p_1) = (p_1 + 1)\sigma(x) - p_1\sigma(\frac x{p_1})\]

这样，在知道\(p_1\)和\(x\)之后我们就可以在枚举质数时完成计算。

由于这些函数都是积性函数，所以当\(p \nmid x\)时的处理方法都是一样的。

Code

inline void init(){    sig[1] = 1;    for(int i = 2;i < N;i++)    {        if(!vis[i])        {            prime[++tot] = i;            sig[i] = i + 1;        }        for(int j = 1;j <= tot && i * prime[j] < N;j+)        {            vis[i * prime[j]] = 1;            if(!(i % prime[j]))            {                sig[i * prime[j]] = ((prime[j] + 1) * sig[i] % MOD - prime[j] * sig[i / prime[j]] % MOD + MOD) % MOD;                break;            }            sig[i * prime[j]] = sig[i] * (prime[j] + 1);        }    }

因数个数d

同理，首先发现当\(p\)为质数时，因数个数\(d(p) = 2\)

对于\(x\)，唯一分解后得到

\[d(x) = (c_1 + 1)(c_2 + 1)...(c_k + 1)\]\[d(x) = (c_1 + 1)T\]\[d(x * p_1) = (c_1 + 2)T\]\[d(\frac x{p_1}) = c_1T\]

所以

\[d(x * p_1) + d(\frac x{p_1}) = 2d(x) \]

所以最终得到：

\[d(x * p_1) = 2d(x) - d(\frac x{p_1})\]

Code

inline void init(){    d[1] = 1;    for(int i = 2;i < N;i++)    {        if(!vis[i])        {            prime[++tot] = i;            d[i] = 2;        }        for(int j = 1;j <= tot && i * prime[j] < N;j+)        {            vis[i * prime[j]] = 1;            if(!(i % prime[j]))            {                d[i * prime[j]] = 2 * d[i] - d[i / prime[j]];                break;            }            d[i * prime[j]] = d[i] * 2;        }    }}

关键词：