（转载）线性时间内筛素数和欧拉函数

转载

mob60475705a319 2014-06-23 12:46:00

线性时间内筛素数和欧拉函数

最近学习了一种筛素数的方法，能够在 O(n) 的时间内完成。一开始不能理解其中的一句话，搜索了很久，大部分的结果都是一群人在csdn上面卖萌。现在终于理解了，记下来备忘。下面的内容不卖萌，但是也不严谨。

算法的实质

这个算法在1978年，由 David Gries 和 Jayadev Misra 在 Communication of the ACM上提出来，可以在ACM Digital Library看到，价格是 $15.00 。当然，我们只是学习用途，所以可以在找找，比如说：A Linear Sieve Algorithm for Finding Prime Numbers。因为ACM的版权限制还是比较严格的，这里就不放上pdf，如果上面的链接不起作用请自行搜索。

这个算法的核心思想是：每一个合数可以被唯一地表示成它的一个最小质因子和另外一个数的乘积。证明略。

伪代码：

S = {2, ..., n}
for (p = 2; p*p <= n; p = next(S, p)):
    print p
    for (i = p; i*p <= n; i = next(S, i)):
        for (j = i*p; j <= n; j = j*p):
            S = S - {j}

next(S, i) is a function defined only for i ∈ S such that 
there is an integer larger than i in S; it yields the next 
larger integer in S.

通俗的说，就是对于每一个质数p，枚举i，删除i*p^1, i*p^2, i*p^3, ...。

这个是筛素数的一个很原始的方案，也就是把质数的所有倍数删掉。看起来很暴力，但是很容易证明它的正确性，更容易证明它的时间复杂度是O(n)的：因为每一个元素最多被删除一次，并且没有新的元素插入S中！法

代码

下面代码就是带有计算欧拉函数的线性筛素数。代码原型的起源已经无从考证，可以作出一个合理的揣测，是某位搞OI或者ACM/ICPC的神牛第一次写出来的。

bool com[MAXN];
int primes, prime[MAXN], phi[MAXN];

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

最难理解的是这句话：

if (i % prime[j] == 0) break;

要理解这句话，（顺便不严谨地）证明这个算法的时间复杂度和正确性，要从下面两个方面：

每个数至少被访问一次
每个数至多被访问一次

每个数至少被访问一次

对于质数，一定会在i的循环中访问到，并确定为质数。

对于合数，因为每一个合数都可以表示成它最小的质因数和另一个数的乘积，而我们枚举了所有的另一个数（也就是i），所以它一定会被它的最小质因数筛掉。

每个数至多被访问一次

对于质数，不可能在j的循环中被访问到，因此仅会在i的循环中被访问到恰好一次。

对于合数，对于i = i1 = p * a，因为在i1 % prime[j1] == 0时break，所以不可能出现一个数x = i1 * prime[k] = p * a * prime[k] (k > j1)在i = i1, j = k的时候被筛掉一次，又在i = a * prime[k]的时候被p给筛掉的情况。

证毕

综上所述，每个数被访问一次且仅访问一次！因此整个算法的复杂度是O(n)的。

ps：我的pascal代码

 1 procedure getfai;
 2  var i,j,k:longint;
 3  begin
 4  tot:=0;
 5  fillchar(fai,sizeof(fai),0);
 6  for i:=2 to n do
 7   begin
 8   if fai[i]=0 then
 9    begin
10     inc(tot);
11     p[tot]:=i;
12     fai[i]:=i-1;
13    end;
14   for j:=1 to tot do
15    begin
16     k:=i*p[j];
17     if k>n then break;
18     if i mod p[j]=0 then
19      begin
20       fai[k]:=fai[i]*p[j];
21       break;
22      end
23     else
24       fai[k]:=fai[i]*(p[j]-1);
25    end;
26  end;
27 end;

View Code