浅谈积性函数的线性筛法

前置知识

数论函数及相关基本定义

素数的线性筛

线性筛

线性筛可以在严格$O(n)$的时间内筛出积性函数的值，

它有常见的套路

假设$n = p_1^{a_1} p_2^{a_2} \dots p_k^{a_k}$

如果我们能快速得出$f(p_i)$和$f(p_i^{k+1})$的取值，那么直接套板子就行了

在下文中如无特殊说明，默认$p_i$表示$n$质因数分解之后第$i$个质数，$a_i$表示$p_i$的指数

常见的有以下几种

线性筛素数

比较简单，这也是筛其他积性函数的基础

#include<cstdio>
const int MAXN = 1e4 + 10;
int N, prime[MAXN], vis[MAXN], tot;
void GetPrime(int N) {
    vis[1] = 1;
    for(int i = 2; i <= N; i++) {
        if(!vis[i]) prime[++tot] = i;
        for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j])) break;
        }
    }
}
int main() {
    GetPrime(1e3);
    for(int i = 1; i <= tot; i++)
        printf("%d ", prime[i]);
    return 0;
}

线性筛莫比乌斯函数

这个也是比较常见的

根据莫比乌斯函数的定义

$$\mu =\begin{cases}\left( -1\right) ^{k}\left( n=p_{1}p_{2}\ldots p_{k}\right) \\ 0\left( \exists P^{2}|n\right) \\ 1\left( n=1\right) \end{cases}$$

直接筛就可以了

#include<cstdio>
const int MAXN = 1e4 + 10;
int N, prime[MAXN], vis[MAXN], mu[MAXN], tot;
void GetMu(int N) {
    vis[1] = mu[1] = 1;
    for(int i = 2; i <= N; i++) {
        if(!vis[i]) prime[++tot] = i, mu[i] = -1;
        for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j])) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = mu[i] * mu[prime[j]];
            //根据莫比乌斯函数的定义，这里也可以写为
            //mu[i * prime[j]] = -mu[i];
        }
    }
}
int main() {
    GetMu(1e3);
    for(int i = 1; i <= tot; i++)
        printf("%d ", mu[i]);
    return 0;
}

 

线性筛欧拉函数

这个貌似更常用一点qwq。

我在以前的文章中也详细的讲过，这里只放一下代码

#include<cstdio>
const int MAXN = 1e4 + 10;
int N, prime[MAXN], vis[MAXN], phi[MAXN], tot;
void GetPhi(int N) {
    vis[1] = phi[1] = 1;
    for(int i = 2; i <= N; i++) {
        if(!vis[i]) prime[++tot] = i, phi[i] = i - 1;
        for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j])) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * phi[prime[j]];
        }
    }
}
int main() {
    GetPhi(1e3);
    for(int i = 1; i <= tot; i++)
        printf("%d ", phi[i]);
    return 0;
}

线性筛约数个数

这个就稍微高端一点了，按照上面的套路，我们只需要考虑最小的质因子对每个数的贡献

$n = p_1^{a_1} p_2^{a_2} \dots p_k^{a_k}$

设$d(i)$表示$i$的约数个数

那么根据约数定理

$d(i) = \prod_{i = 1}^k (a_i + 1) $

记$a(i)$表示$n$的最小的质因子($a_1$)的指数。

$d(p_i) = 2$

当$i % p_j = 0$时，考虑$i * p_j$，实际上也就是$a_i$的指数多了$1$

我们先除去原来的，再加上新的就OK了

#include<cstdio>
const int MAXN = 1e4 + 10;
int N, prime[MAXN], vis[MAXN], D[MAXN], a[MAXN], tot;
void GetD(int N) {
    vis[1] = D[1] = a[1] = 1;
    for(int i = 2; i <= N; i++) {
        if(!vis[i]) prime[++tot] = i, D[i] = 2, a[i] = 1;
        for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j])) {
                D[i * prime[j]] = D[i] / (a[i] + 1) * (a[i] + 2);
                a[i * prime[j]] = a[i] + 1;
                break;
            }
            D[i * prime[j]] = D[i] * D[prime[j]];
            a[i * prime[j]] = 1;
        }
    }
}
int main() {
    GetD(1e3);
    for(int i = 1; i <= tot; i++)
        printf("%d ", D[i]);
    return 0;
}

线性筛约数和

同样，按照上面的套路考虑

$n = p_1^{a_1} p_2^{a_2} \dots p_k^{a_k}$

设$SD(i)$表示$i$的约数和

$SD(n) = \prod_{i = 1}^k (\sum_{j = 1}^{a_i} p_i^j)$

$sum(i)$表示$i$最小的质因子的贡献，即$sum(i) = \sum_{i = 1}^{a_1}p_1^j$

$low(i)$表示$i$最小质因子的指数，$low(i) = a_1$

有了这三个我们就可以转移了

同样是考虑$i$的最小的因子对答案的贡献，应该比较好推，看代码吧

#include<cstdio>
const int MAXN = 1e4 + 10;
int N, prime[MAXN], vis[MAXN], SD[MAXN], sum[MAXN], low[MAXN], tot;
void GetSumD(int N) {
    vis[1] = SD[1] = low[1] = sum[1] = 1;
    for(int i = 2; i <= N; i++) {
        if(!vis[i]) prime[++tot] = i, sum[i] = SD[i] = i + 1, low[i] = i;
        for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j])) {
                low[i * prime[j]] = low[i] * prime[j];
                sum[i * prime[j]] = sum[i] + low[i * prime[j]];
                SD[i * prime[j]] = SD[i] / sum[i] * sum[i * prime[j]];
                break;
            }
            low[i * prime[j]] = prime[j];
            sum[i * prime[j]] = prime[j] + 1;
            //这里low和sum不是积性函数 
            SD[i * prime[j]] = SD[i] * SD[prime[j]];
        }
    }
}
int main() {
    GetSumD(1e3);
    for(int i = 1; i <= tot; i++)
        printf("%d ", SD[i]);
    return 0;
}

 

非常见积性函数的筛法

很多情况下我们会遇到求两个积性函数狄利克雷卷积的情况

很显然，这个函数也是积性函数，我们考虑如何求得

为了方便筛，我们需要把问题无限简化，

设$low(i)$表示$p_1^{a_1}$

考虑筛法中最关键的地方

$i \% p_j = 0$，

有了$low(i)$，此时我们需要分两种情况讨论

1. $low(i) = i$，此时$i$一定是某个素数的幂的形式（否则就会break掉）

这里就用到了我最开始说的那个套路

如果我们能快速的利用$f(p_i^{k})$更新出$f(p_i^{k + 1})$，那这个素数就很容易筛了

2. $low(i) \not = i$，那么$i / low(i)$一定与$low(i) * p_j$是互质的，我们可以直接利用积性函数的性质去更新

C++版的伪代码

vis[1] = low[1] = 1; H[1] = 初始化 
for(int i = 2; i <= N; i++) {
    if(!vis[i]) prime[++tot] = i, mu[i] = -1, H[i] = 质数的情况, low[i] = i;
    for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
        vis[i * prime[j]] = 1;
        if(!(i % prime[j])) {
            low[i * prime[j]] = (low[i] * prime[j]); 
            if(low[i] == i) H[i * prime[j]] = 特殊判断;
            else H[i * prime[j]] = H[i / low[i]] * H[prime[j] * low[i]];
            break;
        } 
        H[i * prime[j]] = H[i] * H[prime[j]];
        low[i * prime[j]] = prime[j];
    }
}

 

参考资料

积性函数与线性筛

线性筛约数个数和、约数和

线性筛,积性函数,狄利克雷卷积,常见积性函数的筛法