求素数

发表于

码一波求1-n之间素数的办法。。。。

硬刚解法–n,n/2

1
2
3
4
5
6
7
8
9
10
int main()
{
  int i,n;
  while(~scanf("%d",&n))
  {   for(i=2;i<n;i++)
         if(n%i==0)    break; 
      if(i==n)    printf("YES\n");
      else           printf("NO\n");
  }
}

这是理所当然的想法，按照素数的定义，除了1和它本身没有其他的因数，就是素数。
这种解法的缺点就是循环规模n稍微大点，运行时间。。。
只需要简单思索就能考虑到，n的因数是成对出现的，如(2,n/2),(3,n/3)等等，因此只枚举i最大到sqr(n)即可退出循环。这里需要注意浮点误差，一般我们用：

1
for(i=2;i<=(sqrt(n)+0.5);i++)

硬刚解法估计只能在入门题中使用。。。也就是对于某些数据量小且n小的题目。。还是有奇效的//我就是在说你校的比赛题
emmm…但不是所有题都对菜鸡这么友好的。。我们还可以在时间上改进算法。

普通筛选法–埃拉托色尼筛法

基本思想：素数的倍数一定不是素数
实现方法：用一个长度为N+1的bool数组保存信息（0表示素数，1表示非素数），先假设所有的数都是素数（初始化为0），从第一个素数2开始，把2的倍数都标记为非素数（置为1），一直到大于N；然后进行下一趟，找到2后面的下一个素数3，进行同样的处理，直到最后，数组中依然为0的数即为素数。
说明：整数1特殊处理即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
bool judge[n+1]={0};
int tot = 0;
for (int i = 2; i <= n; ++i)
{
  if (!judge[i])
  {
    prime[tot++] = i;
  }
  for (int j = i+i; j <= n; j += i)
  {
    judge[j] = 1;
  }
}

当然打表打一次就够了。。。判断n以内一个数是不是指数直接judge[]就可以了。
但是我们大概算一下就会发现时间复杂度也不少。。。比如1e8,就是1e8*(1/2+1/3+1/5+1/7+1/11+…..)循环次数，前三个数加起来就以及1e8了，还不算判断时间。。。emmm 当然慢也很好解释，比如说一个合数30030，它分别在(2,3,5,7,9,11,13)时筛了一共是七次。。。emmm….更大的就更不用说了。。。

能不能再给力点呢？

线性筛法–欧拉筛法

原理：每个合数只会被它的最小质因数筛去，因此每个数只会被标记一次，所以时间复杂度近似是O(n)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
const int N=1000000;
bool a[N+1];
vector<int> prime;
void get_prime()
{
    memset(a,true,sizeof(a));
  for(int i=2;i<=N;i++)
  {    
      if(a[i]) prime.push_back(i);
    for(int j=0;j<prime.size()&&i*prime[j]<=N;j++)
    { 
        a[i*prime[j]]=false;
      if(!(i%prime[j])) break;
    }
  }
}

节省空间的写法(少一个bool数组)：

1
2
3
4
5
6
7
8
9
10
11
const int MAXN = 1e7 + 5;
void getPrime(){
  memset(prime, 0, sizeof(prime));
  for (int i = 2;i <= MAXN;i++) {
    if (!prime[i])prime[++prime[0]] = i;
    for (int j = 1;j <= prime[0] && prime[j] <= MAXN / i;j++) {
      prime[prime[j] * i] = 1;
      if (i%prime[j] == 0) break;
    }
  }
}

容斥原理

1e8的情况是筛选法完全无法满足的，但是还是考虑a * b = c的情况，1e8只需要考虑10000以内的素数p[10000]，然后每次先减去n / p[i]，再加上n / (p[i] * p[j])再减去n / (p[i] * p[j] * p[k])以此类推…于是就可以得到正确结果了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
#include <cmath>
#include <cstdio>
using namespace std;
 
const int maxn = 10005;
int sqrn, n, ans = 0;
bool vis[maxn];
int pri[1500] = {0};
void init(){
    vis[1] = true;
    int k = 0;
    for(int i = 2; i < maxn; i++){
        if(!vis[i]) pri[k++] = i;
        for(int j = 0; j < k && pri[j] * i < maxn; j++){
            vis[pri[j] * i] = true;
            if(i % pri[j] == 0) break;
        }
    }
}
void dfs(int num, int res, int index){
    for(int i = index; pri[i] <= sqrn; i++){
        if(1LL * res * pri[i] > n){
            return;
        }
        dfs(num + 1, res * pri[i], i+1);
        if(num % 2 == 1){
            ans -= n / (res * pri[i]);
        }else{
            ans += n / (res * pri[i]);
        }
 
        if(num == 1) ans++;
    }
}
int main(){
    init();
    while(~scanf("%d",&n) && n){
        ans = n;
        sqrn = sqrt((double)n);
        dfs(1,1,0);
        printf("%d\n",ans-1);
    }
    return 0;
}

Meissel-Lehmer算法

最后介绍的这个算法可以说是黑科技级别的，它能够在3s内求出1e11之内的素数个数。据说还有在300ms内求出1e11的个数的。可以参考wiki里面原理。然后给出来自Codeforces 665F题目里面的代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#define MAXN 100    // pre-calc max n for phi(m, n)
#define MAXM 10010 // pre-calc max m for phi(m, n)
#define MAXP 40000 // max primes counter
#define MAX 400010    // max prime
#define setbit(ar, i) (((ar[(i) >> 6]) |= (1 << (((i) >> 1) & 31)))) 
#define chkbit(ar, i) (((ar[(i) >> 6]) & (1 << (((i) >> 1) & 31))))
#define isprime(x) (( (x) && ((x)&1) && (!chkbit(ar, (x)))) || ((x) == 2))
namespace pcf {
  long long dp[MAXN][MAXM];
  unsigned int ar[(MAX >> 6) + 5] = { 0 };
  int len = 0, primes[MAXP], counter[MAX];
  void Sieve() {
    setbit(ar, 0), setbit(ar, 1);
    for (int i = 3; (i * i) < MAX; i++, i++) {
      if (!chkbit(ar, i)) {
        int k = i << 1;
        for (int j = (i * i); j < MAX; j += k) setbit(ar, j);
      }
    }
    for (int i = 1; i < MAX; i++) {
      counter[i] = counter[i - 1];
      if (isprime(i)) primes[len++] = i, counter[i]++;
    }
  }
  void init() {
    Sieve();
    for (int n = 0; n < MAXN; n++) {
      for (int m = 0; m < MAXM; m++) {
        if (!n) dp[n][m] = m;
        else dp[n][m] = dp[n - 1][m] - dp[n - 1][m / primes[n - 1]];
      }
    }
  }
  long long phi(long long m, int n) {
    if (n == 0) return m;
    if (primes[n - 1] >= m) return 1;
    if (m < MAXM && n < MAXN) return dp[n][m];
    return phi(m, n - 1) - phi(m / primes[n - 1], n - 1);
  }
  long long Lehmer(long long m) {
    if (m < MAX) return counter[m];
    long long w, res = 0;
    int i, a, s, c, x, y;
    s = sqrt(0.9 + m), y = c = cbrt(0.9 + m);
    a = counter[y], res = phi(m, a) + a - 1;
    for (i = a; primes[i] <= s; i++) res = res - Lehmer(m / primes[i]) + Lehmer(primes[i]) - 1;
    return res;
  }
}

引申–求欧拉函数

在数论，对正整数n，欧拉函数是小于n的正整数中与n互质的数的数目。此函数以其首名研究者欧拉命名，它又称为Euler’s totient function、φ函数、欧拉商数等。

//来自约三百年前的碾压

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#define MAXN 100005
#define MAXL 1299710
int prime[MAXN];
int check[MAXL];
int phi[MAXL];
int tot = 0;
phi[1] = 1;
memset(check, 0, sizeof(check));
for (int i = 2; i < MAXL; ++i)
{
  if (!check[i])
  {
    prime[tot++] = i;
    phi[i] = i - 1;
  }
  for (int j = 0; j < tot; ++j)
  {
    if (i * prime[j] > MAXL)
    {
      break;
    }
    check[i*prime[j]] = 1;
    if (i % prime[j] == 0)
    {
      phi[i*prime[j]] = phi[i] * prime[j];
      break;
    }else
    {
      phi[i*prime[j]] = phi[i] * (prime[j]-1);
    }
  }
}

参考：

https://www.cnblogs.com/grubbyskyer/p/3852421.html
http://blog.csdn.net/snow_me/article/details/52588819