泡沫

博客园 首页 联系 订阅 管理

如何快速求出区间范围[0,x]内的素数?

Brute Force 总是那么的让人觉得亲近。

枚举每一个数,判断该数是否为素数?

bool isPrime( int x ) {
    if(x < 2) return 0;
    if(x == 2) return 1;
    for ( int i=2; i*i<=x; ++i )
      if(0 == x%i) return 0;
    return 1;
} 
// 用prime[0]存素数个数
for ( int i=0; i<=x; ++i ) if( isPrime(i) ) prime[++prime[0] ] = i;

暴力解决问题在很多时候都会存在效率低下这个问题。那么对于求素数这个问题如何提高效率呢?

其实很多数都是有公共因子的(公约数)。如果直接把这个公约数拿出来去踢出其倍数,被踢出的数都不会是素数的。例如:

2 可以踢出: 4 6 8 10 12 。。。2*x (x=2,3,4,5,6,7,8....)

3 可以踢出: 6 9 12 15 18。。。3*x(x=2,3,4,5,6,7,8....)

5 可以踢出: 10 15 20 25。。。5*x(x=2,3,4,5,6,7,8.....)

......... 依次做下去,剩余的数就是素数了。

int prime[3438];
bool isPrime[40007];

for ( int i=1; i<=x; ++i ) // 初始化所有的数都是素数
  isPrime[i] = 1;
for ( int i=2; i<=x; ++i ) if(isPrime[i] ) { // 该数是素数
    for ( int j=i+i; j<x; j+=i ) //踢出不是素数的数
      isPrime[j] = 0;
}
for ( int i=2; i<=x; ++i ) // 存取素数
  if(isPrime[i] ) prime[++prime[0] ] = i;

这样看貌似效率是提高了不少额。可是仔细分析上面的踢出非素数的过程会发现,有些数被踢出了好几次,为什么会出现这个问题呢?

x 可以踢出: x*2 x*3 x*4  x*(x-1) x*x x*(x+1).....

不难发现其实x*2 x*3 x*4....x*(x-1)这些被更小的因子给踢出了。也就是说对于x应该踢出的数是从

x*x开始的。那么是不是这就已经优化到顶了呢? 没有哦。

开始考虑x*x是偶数,还是奇数。如果是偶数的话,x必然是偶数的。这中情况可以直接用素因子2直接把这些数踢出掉。

x*x是奇数,x必然是奇数了。 x*x+x是什么数呢? 是偶数。奇数+奇数=偶数。因为奇数%2=1,两个1相加就是偶数了。这样的话又会造成部分偶数被踢出多次。如何解决这个问题。呵呵,只需对x*x+x+x就好了,因为x*x+x是偶数 x*x+x+x,就是奇数了, 接下来的奇数就是x*x+2*(x+x)了。优化到此结束。

int prime[3438];
bool isPrime[40007];

void getPrime( int x ){
   for ( int i=1; i<x; i+=2 )
      isPrime[i] = 1, isPrime[i-1] = 0; // 首先踢出2的倍数
   prime[prime[0]=1 ] = 2; // 存取第一个素数
   for ( int i=3; ; i+=2 ) if(isPrime[i] ) { // 只考虑奇数了
       int j = i*i, k = i+i;
       if(j >= x) break;
       while(j < x ) {
           isPrime[j] = 0;  j += k;
       }
   }// 存素数
   for ( int i=3; i<x; i += 2 )
     if(isPrime[i]) prime[++prime[0] ] = i;
}

对次求解素数已经完毕了。

 接下来是将x分解成素因子的乘积

求x = p[1]^cnt[1]*p[2]^cnt[2]*p[3]^cnt[3]。。。。p[...]^cnt[...]  p[]是素因子

int p[3438], cnt[3438];

void getPrimeDivisor( int x ) {
     p[0] = cnt[0] = 0; int t;
     for ( int i=1; prime[i]*prime[i]<=x  && i<=prime[0]; ++i ) {
         t = 0;
         while( x%prime[i] == 0 ) {
              ++t; x /= prime[i];
         }
         if( t ) p[++p[0] ] = prime[i], cnt[++cnt[0] ] = t;
     }
     if(x > 1) p[++p[0] ] = x, cnt[++cnt[0] ] = 1;
// 此时x必定是个素数。 因为在sqrt(x)内都没有找到其因子 };

此时x的因子个数=(cnt[1]+1)*(cnt[2]+1)*(cnt[3]+1)*(cnt[4]+1)....

额,只知道因子个数而已。不是很好玩ing。能不能求出x的所有约数呢? 恭喜你,这是可以滴。

利用乘法原理对约数集合乘以一个约数得到新的集合,

int divisor[1500];

void getDivisor( int x ) { // 或得x的所有因子
     getPrimeDivisor(x); // 对x进行素因子分解
     divisor[0] = 1; // 存取因子个数
     divisor[1] = 1;
     for ( int i=1; i<=p[0]; ++i ) {
         int nowNum = divisor[0];
         int base = 1;
         for ( int j=1; j<=cnt[i]; ++j )
         {
             base *= p[i];
             for ( int k=1; k<=divisor[0]; ++k )
               divisor[++nowNum ] = divisor[k]*base;
         }
         divisor[0] = nowNum;
     }
}

 

下面通过一道题来测试代码的正确性:

hdoj 4432 http://acm.hdu.edu.cn/showproblem.php?pid=4432

题意是求n的所有约数,并把每个约数转化为m进制,把每个进制为的平方加起来。再把最后的结果转化为m进制。简单吧。

View Code
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

int prime[3438];
bool isPrime[40007];

void getPrime( int x ){
   for ( int i=1; i<x; i+=2 )
      isPrime[i] = 1, isPrime[i-1] = 0;
   prime[prime[0]=1 ] = 2;
   for ( int i=3; ; i+=2 ) if(isPrime[i] ) {
       int j = i*i, k = i+i;
       if(j >= x) break;
       while(j < x ) {
           isPrime[j] = 0;  j += k;
       }
   }
   for ( int i=3; i<x; i+=2 )
     if(isPrime[i]) prime[++prime[0] ] = i;
}

int p[3438], cnt[3438];

void getPrimeDivisor( int x ) {
     p[0] = cnt[0] = 0; int t;
     for ( int i=1; prime[i]*prime[i]<=x  && i<=prime[0]; ++i ) {
         t = 0;
         while( x%prime[i] == 0 ) {
              ++t; x /= prime[i];
         }
         if(t ) p[++p[0] ] = prime[i], cnt[++cnt[0] ] = t;
     }
     if(x > 1) p[++p[0] ] = x, cnt[++cnt[0] ] = 1;
};

int divisor[1500];

void getDivisor( int x ) {
     getPrimeDivisor(x);
     divisor[0] = 1;
     divisor[1] = 1;
     for ( int i=1; i<=p[0]; ++i ) {
         int nowNum = divisor[0];
         int base = 1;
         for ( int j=1; j<=cnt[i]; ++j )
         {
             base *= p[i];
             for ( int k=1; k<=divisor[0]; ++k )
               divisor[++nowNum ] = divisor[k]*base;
         }
         divisor[0] = nowNum;
     }
}

int get( int x, int bit ){
    int ret = 0;
    while(x ) {
        ret += (x%bit)*(x%bit);
        x /= bit;
    }
    return ret;
};

void print( int x, int bit ) {
    int a[100]={0};
    while( x ) {
       a[++a[0] ] = x%bit;
       x /= bit;
    }
    for ( int i=a[0]; i>0; --i)
      if(a[i] > 9) printf("%c", a[i]-10+'A');
      else printf("%d", a[i]); puts("");
}

int main(){
    getPrime(32000);
    int n, m;
    while( ~scanf("%d%d", &n, &m) ) {
         getDivisor(n);
         int ans = 0;
         for ( int i=1; i<=divisor[0]; ++i ) {
             ans += get(divisor[i], m);
         }
         print(ans, m);
    }
}

 

 

关于素数的运用还有很多的,比如求数x的所有约数的和。初等数论还是挺有意思的,可惜没时间整了,悲剧啊。

posted on 2013-04-23 20:08  木-天空  阅读(3604)  评论(0编辑  收藏  举报