[转]素数基础篇 之 素数的个数

   转自:http://hi.baidu.com/vfxupdpaipbcpuq/item/0cbe522c3c9e4bcaa5275a12

1.小区间素数个数,百万级别。

    简述:有的时候,我们需要知道某个特定区间的素数(区间大小较小,但数可能很大)。那么数组就开不下,这时候我们仍然可以使用筛法,只是所有的下标都进行了偏移。大家理解下面这段代码可以先用普通筛法写,然后数组下标集体移动即可。

const int maxn = 100000;
int PrimeList[maxn];
int PrimeNum;
bool IsNotPrime[maxn]; // IsNotPrime[i] = 1表示i + L这个数是素数.
void SegmentPrime(int L, int U)
{ // 求区间[L, U]中的素数.
     int i, j;
     int SU = sqrt(1.0 * U);
     int d = U - L + 1;
     for (i = 0; i < d; i++) IsNotPrime[i] = 0; // 一开始全是素数.
     for (i = (L % 2 != 0); i < d; i += 2) IsNotPrime[i] = 1; // 把偶数的直接去掉.
     for (i = 3; i <= SU; i += 2)
     {
           if (i > L && IsNotPrime[i - L]) continue; // IsNotPrime[i - L] == 1说明i不是素数.
           j = (L / i) * i; // j为i的倍数,且最接近L的数.
           if (j < L) j += i;
           if (j == i) j += i; // i为素数,j = i说明j也是素数,所以直接 + i.
          j = j - L;
          for (; j < d; j += i) IsNotPrime[j] = 1; // 说明j不是素数(IsNotPrime[j - L] = 1).
     }
     if (L <= 1) IsNotPrime[1 - L] = 1;
     if (L <= 2) IsNotPrime[2 - L] = 0;
     PrimeNum = 0;
     for (i = 0; i < d; i++) if (!IsNotPrime[i]) PrimeList[PrimeNum++] = i + L;
}

 

2.大区间素数个数,千万、亿级别。

    了解了基本的判断方法后,你是不是有个疑问:“我们能判断素数的个数吗?”总所周知,素数的个数是无限的,且没有固定的公式…但如果我们只要判断[a, b]区间(a, b范围为1到1亿)内的素数的个数呢?

     首先,我们可以想到,如果要求的素数个数区间[a, b],当区间长度比较小(10^6内),我们可以用筛法求出区间内的所有的素数,然后统计个数即可。
但如果区间长度很长或者要求询问的次数很多,那该怎么办呢? [a,b]区间内素数的个数 = [1, b]的个数 - [1, a - 1]的个数,所以我们这里只讨论求[1, a]区间内的素数。以下提供个人的两种方法,时限都是1s内产生结果。如果哪位大牛有更好的方法,大家一起交流下~~

     1. 我们可以扩展上面的思想,当区间小的时候,我们可以很好的求出素数的个数。那我们可以把大的区间划分成一块块小的区间,比如把一个长度为1亿的区间划分1,000个长度为100,000的区间。我们可以利用Miller-Rabin事先把[1, 100000], [100001, 200000], [200001, 300000]的区间内的素数个数统计好,然后存在一个数组中。
完成这步后,思路就比较清晰:对于区间[1, a],可以拆分为一个个长度为100000的小区间([1, 100000], [100001, 200000]…),加上尾部的小区间[c * 100000, a]。前面的小区间只要数组的值相加即可,而后面的小区间[c * 100000, a],长度在100000内,直接用区间的筛法求出素数,统计个数即可。
     代码:参见上一篇文章的Miller-Rabin,区间求素数的代码。
     该方法速度很快,主要时间都花在数组打表上,然后直接存在数组里,求1到1亿的素数个数时间为0.06s。
     评价:优点是方法速度快,且直接套模板即可。缺点是需要事先打表,且代码长度很长(因为要给长度为1000的数组赋初值)。

     2. 第二种方法涉及到容斥原理(inclusion-exclusion principle),容斥原理参见(http://en.wikipedia.org/wiki/Inclusion-exclusion_principle)。
      当一个数是合数,那么它可以分解成几个素数的乘积。如30 = 2 * 3 * 5。我们可以统计合数的个数,然后拿总数减它就是素数的个数(注意还要去掉1的)。我们可以利用类似筛法的原理,去除2的倍数(它们肯定是合数,不包括2),然后去除3的倍数,5的倍数,知道去除到Sqrt(a)的倍数为止。但你会发现6 = 2 * 3,被去除了2次,而这正是容斥原理解决的问题。合数的个数 = 1个素数筛完的合数个数 – 2个素数筛完的合数个数 + 3个素数筛完的合数个数...
而容斥原理的累加过程,即可用DFS来解决。你可能会认为sqrt(1亿) = 10000,其中素数有很多,DFS要跑很长时间。但我们只需要加一些简单的优化即可很大程度地提高程序的效率。

    首先,我们写筛法出1到sqrt(a)的素数表,然后从小到大DFS。
    如果当前的乘积 > a,那么直接退到上一层。
    如果该层的所有乘积不能使总数发生变化(即所有乘积都 > a),那么直接退回第一层。(因为是从小到大,该层下面的乘积必将 > a)
    如果是第一层的所有乘积不能使总数发生变化,那么程序运行结束。(原理同上)
    经过这样优化后,求1到1亿的时间为0.4秒,1到10亿的时间为3.5s。

     核心代码:

void Solve(int index, int lcm, int K)
{
    int i;
    int t, t_temp;
    if (K == 0)
    {
        temp += n / lcm;
        return ;
    }
    for (i = index; i < total - K + 1; i++)
    {
        t = lcm * primelist[i];
        t_temp = temp;
        if (t <= n)
        {
            Solve(i + 1, t, K - 1);
        }
        if (t_temp == temp) return ; // 剪枝:同样道理,说明以后的K - 1个不能组成我们想要的值
    }
}

main()中:
for (k = 1; k <= total; k++)
{ // 计算size中选k个的总数.
    temp = 0;
    Solve(0, 1, k);
    if (temp == 0) break; // 说明最小的k个乘积都大于n了,那么可以直接break了.
    if (k & 1) ans += temp;
    else ans -= temp;
}

    评价:该方法巧妙地使用了容斥原理来计数,且DFS应用于容斥原理的剪枝十分重要。
     
    提到容斥原理,推荐一道前不久做的题目SRM 453.5 DIV 1 1000(http://www.topcoder.com/stat?c=problem_statement&pm=10420&rd=14174),两题的容斥原理思想差不多,但剪枝方法不同,而且两题的方法交换都会产生超时...(有兴趣一起交流下~~)

 

posted @ 2013-11-04 10:56  zhsl  阅读(1117)  评论(0编辑  收藏  举报