Skip to main content

__debug's Home Keep it simple, stupid

筛法小结 (Eratosthenes/Euler)

Eratosthenes 筛法 (Sieve of Eratosthenes)

由于 思想 非常简单, 故只给出实现.

void eratosthenes_sieve(int n)
{
    totPrimes = 0;
    memset(flag, 0, sizeof(flag));

    int sqrtn = sqrt(n + 0.5);
    for (int i = 2; i <= sqrtn; i++) {
	if (!flag[i]) {
	    primes[totPrimes++] = i;
	    for (int j = i * i; j <= n; j += i) {
		flag[j] = true;
	    }
	}
    }
    for (int i = sqrtn + 1; i <= n; i++) {
	if (!flag[i])
	    primes[totPrimes++] = i;
    }
}

时间复杂度 \(O(n \log\log n)\).

Euler 筛法 (Sieve of Euler)

欧拉筛是一种线性算法, 并且同时可以计算出每个数的 \(\varphi\).

(好吧其实只要是积性函数并且容易求出 \(f(p^{k})\), 都是可以用欧拉筛算的)

做法

回顾经典的 Eratosthenes 筛法, 它可能对同一个质数筛去多次. 那么如果用某种方法使得每个合数只被筛去一次就变成是线性的了.

不妨规定每个合数只用其 最小的 一个质因数去筛, 这便是欧拉筛了.

不妨先看代码:

void euler_sieve(int n)
{
    totPrimes = 0;
    memset(flag, 0, sizeof(flag));

    for (int i = 2; i <= n; i++) {
	if (!flag[i])
	    primes[totPrimes++] = i;
	for (int j = 0; i * primes[j] <= n; j++) {
	    flag[i * primes[j]] = true;
	    if (i % primes[j] == 0)
		break;
	}
    }
}

请仔细体会 i % primes[j] == 0 的含义.

时间复杂度 \(O(n)\).

简单证明

这个看似很简单, 其实还是要注意一下细节的. 搞清了证明其他的问题也就清楚了.

证明分两部分. 首先证每个合数都会被筛到 (正确性), 其次证每个合数只会被筛到一次 (复杂度).

每个合数都会被筛到

设有一合数 \(n = p_1^{k_1}...p_m^{k_m}\) (\(p_i\) 为质数)

则 \(n\) 一定会在 \(i = p_1^{k_1-1}...p_m^{k_m}\) 时被筛去 (此时 \(primes[j] = p_1\)), 因为对于小于 \(p_1\) 的质数, 一定不会被 \(i\) 整除.

每个合数都只会被筛到一次

与上面一样, 还是设有一合数 \(n = p_1^{k_1}...p_m^{k_m}\) (\(p_i\) 为质数)

倘若存在一个质因子 \(p_x(x \ne 1)\) 也筛去了 \(n\), 那么此时 \(i = p_1^{k_1}...p_x^{k_i-1}...p_m^{k_m}\).

  • \(i > p_x\), 此时在内层循环中已经早早地 break 掉了, 因为 \(p_1 \mid i\)
  • \(i < p_x\), 此时 \(p_x\) 还没加进质数表 (顺便一提: 这种情况只有 可能 在 \(x = m\) 时发生)

难以理解的地方

为何不把 i % primes[j] == 0 放在前面?

你可以去试试…… 实践出真知.

放前面的话, 所有的“某个质因子的次数不为 $1$”的合数便会被当成质数. 至于为什么, 请看证明.

为何没有 j < totPrimes ?

  • 当 \(i\) 为质数时, 内层循环会在最后一个质数 (也就是 \(i\) 自己) 终止
  • 当 \(i\) 为合数时, 内层循环会在它的第一个质因数终止

当然加了也没有问题.

顺便把 \(\varphi\) 算出来?

其实这是极简单的.

主要基于以下事实: (很容易通过定义推出来,不妨自己试试)

  • \(n\) 为质数时, \(\varphi(n) = n - 1\)
  • \(p\) 为质数且 \(p \mid n\) 时, \(\varphi(np) = \varphi(n) \times p\)
  • \(p\) 为质数且 \(p \not\mid n\) 时, \(\varphi(np) = \varphi(n) \times (p - 1)\)

直接线性筛.

代码:

void euler_sieve_with_phi(int n)
{
    totPrimes = 0;
    phi[1] = 1;
    memset(flag, 0, sizeof(flag));

    for (int i = 2; i <= n; i++) {
	if (!flag[i]) {
	    primes[totPrimes++] = i;
	    phi[i] = i - 1;
	}
	for (int j = 0; i * primes[j] <= n; j++) {
	    flag[i * primes[j]] = true;
	    if (i % primes[j])
		phi[i * primes[j]] = phi[i] * (primes[j] - 1);
	    else {
		phi[i * primes[j]] = phi[i] * primes[j];
		break;
	    }
	}
    }
}

速度比较

你可能会觉得 \(\log\log n\) 微不足道, 那么你错了. (2333)

实测结果: (精确到微秒, 不打开优化开关)

when n = 10000
    eratosthenes_sieve(1229):     0(us)
    euler_sieve(1229):            0(us)
when n = 100000
    eratosthenes_sieve(9592):     999(us)
    euler_sieve(9592):            0(us)
when n = 1000000
    eratosthenes_sieve(78498):    13004(us)
    euler_sieve(78498):           7004(us)
when n = 10000000
    eratosthenes_sieve(664579):   185130(us)
    euler_sieve(664579):          79067(us)
when n = 100000000
    eratosthenes_sieve(5761455):  2363692(us)
    euler_sieve(5761455):         842592(us)
when n = 1000000000
    eratosthenes_sieve(50847534): 25535159(us)
    euler_sieve(50847534):        8987385(us)

差距还是挺大的呢.

Comments

Comments powered by Disqus