The secret of being a bore is to tell everything

0%

数论筛法

欧拉筛/线性筛

之前的素数判定中,我们曾使用埃拉托斯特尼筛法(Sieve of Eratosthenes),进行素数的筛选。但是这个算法的时间复杂度是$O(n\log{\log{n}})$的,因为每个数都被它筛了它的素因子个数那么多次。

核心思想

欧拉筛法(Euler’s sieve)就是减少这种重复的一种筛法,埃筛的方法是筛去素数的倍数,而欧拉筛的核心思想是筛去最小素因子为它的数。这样一对比下来,因为每个数只有一个最小素因子,所以时间复杂度显然是$O(n)$。

然而,欧拉筛法的实现是有一点技巧的,我当初尝试理解的时候还是费了一番功夫。假设我们有一个数$x$,有唯一质因数分解:$\prod_{k} p_k^{\alpha_k} \ (p_{k-1} < p_k)$。此时欧拉筛会筛除一个$p_i*x$,$p_i \leq p_1$,$p_1$就是最小素因子,也就是说筛去的这个数的最小素因子会小于等于$x$的最小素因子,且是$x$的倍数。现在我们就要对于所有的$x$都筛掉这些倍数。

那么为什么每个数只会被筛一次呢?根据唯一分解定理,任何一个数的质因数排序后的序列一定是唯一的,那么假设$a\neq b$,并且他们删掉最小素因子后的数相同,那么它们的最小素因子一定不同。那么他们会被那个除掉最小素因子的数筛掉,然后就不会再被触及。

现在又有一个问题,怎么保证$x\in[2,p-1]$的合数全部被筛掉,没有遗漏呢?还是那个思路,因为小于$p$的数的最小素因子一定小于$p$,所以它一定会被小于$x/最小素因子$的数筛掉,这就保证了正确性。

实现的时候,我们筛到一个素数就要把它放到一个数组里面,此时素数的大小顺序非常重要,因为我们对于$x$需要查找所有小于等于$p_1$的素数,同时在等于$p_1$的时候停止查找。

代码:

1
bool nop[MAXN];
2
int primes[MAXN / ln(MAXN)];
3
void sieve() {
4
    int tot = 0;
5
    for (int i = 2; i <= n; i++) {
6
        if (!nop[i]) primes[tot++] = i;
7
        for (int j = 0; j < tot && i * primes[j] <= n; j++) {
8
            nop[i * primes[j]] = true;
9
            // i 整除 p_j 代表此时p_j已经是i的最小素因子
10
            if (!(i % primes[j])) break;
11
        }
12
    }
13
}

扩展应用

欧拉筛在数论筛法里面用途非常广泛,尤其是对于积性函数,以欧拉函数筛为例,我们令$\varphi(x)=n\prod_k(\frac{p_k-1}{p_k})$,如果$x$是质数,那么$\varphi(x)=x-1$。同时$\varphi(pq) = \varphi(p) * \varphi(q) \ (p \perp q)$,也就是积性函数性质。

于是我们可以在筛到质数的时候把它的$\varphi(x)$设为$x-1$,否则,我们就需要去计算筛到合数的时候$\varphi(x)$值的变化。因为每次我们只加入一个质因数$p$,如果这个质因数不在原数的质因数分解中,那么显然答案要乘以$\varphi(p)$,否则就乘以$p$本身。

代码:

1
bool nop[MAXN];
2
int primes[MAXN / ln(MAXN)], phi[MAXN];
3
void sieve() {
4
    int tot = 0;
5
    phi[1] = 1;
6
    for (int i = 2; i <= n; i++) {
7
        if (!nop[i]) primes[tot++] = i, phi[i] = i - 1;
8
        for (int j = 0; j < tot && i * primes[j] <= n; j++) {
9
            nop[i * primes[j]] = true;
10
            if (!(i % primes[j])) {
11
                phi[i * primes[j]] = primes[j] * phi[i];
12
                break;
13
            } else {
14
                phi[i * primes[j]] = phi[primes[j]] * phi[i];
15
            }
16
        }
17
    }
18
}