1. 程式人生 > >超高速計算n以內素數個數(百億內3毫秒解決)

超高速計算n以內素數個數(百億內3毫秒解決)

     判斷n以內素數個數有很多演算法,最簡單的是迴圈直接判斷,這個效率不用說,n稍大就不行了。最流行的是篩選法,原理就是定義一個素數標誌位表,初始為1,遇到一個數如果對應標誌位為1判斷這個數是不是素數,是將該為置1,不是放0,然後將他的倍數位置全部置0,然後繼續。。這個效率還是比較快的,但是計算到10^8時候需要3s左右了,對於一般要求基本夠了,但是對於ACM裡面對時間要求很嚴還是不夠。可以對帥選法進行優化,不如偶數直接跳過,以後直接加偶數倍,甚至加入移位運算判斷是不是3的倍數,5的倍數等等,最後基本勉強在ACM要求的時間之內。下來介紹一種逆天的演算法:MEISSEL-LEHMER,布吉島的可以百度下。不容易看懂。。。。。。

先看下效果絕對碉堡!!!!

程式碼如下:

#include    <stdio.h>
#include    <string.h>
#include    <stdlib.h>
#include    <time.h>
#include    <math.h>
 
__int64 *primarr, *v;
__int64 q = 1, p = 1;
 
//π(n)
__int64 pi(__int64 n, __int64 primarr[], __int64 len)
{
    __int64 i = 0, mark = 0;
    for (i = len - 1; i > 0; i--) {
        if (primarr[i] < n) {
            mark = 1;
            break;
        }
    }
    if (mark)
        return i + 1;
    return 0;
}
 
//Φ(x,a)
__int64 phi(__int64 x, __int64 a, __int64 m)
{
    if (a == m)
        return (x / q) * p + v[x % q];
    if (x < primarr[a - 1])
        return 1;
    return phi(x, a - 1, m) - phi(x / primarr[a - 1], a - 1, m);
}
 
__int64 prime(__int64 n)
{
    char *mark;
    __int64 mark_len;
    __int64 count = 0;
    __int64 i, j, m = 7;
    __int64 sum = 0, s = 0;
    __int64 len, len2, len3;
 
    mark_len = (n < 10000) ? 10002 : ((__int64)exp(2.0 / 3 * log(n)) + 1);
 
    //篩選n^(2/3)或n內的素數
    mark = (char *)malloc(sizeof(char) * mark_len);
    memset(mark, 0, sizeof(char) * mark_len);
    for (i = 2; i < (__int64)sqrt(mark_len); i++) {
        if (mark[i])
            continue;
        for (j = i + i; j < mark_len; j += i)
            mark[j] = 1;
    }
    mark[0] = mark[1] = 1;
 
    //統計素數數目
    for (i = 0; i < mark_len; i++)
        if (!mark[i])
            count++;
 
    //儲存素數
    primarr = (__int64 *)malloc(sizeof(__int64) * count);
    j = 0;
    for (i = 0; i < mark_len; i++)
        if (!mark[i])
            primarr[j++] = i;
 
    if (n < 10000)
        return pi(n, primarr, count);
 
    //n^(1/3)內的素數數目
    len = pi((__int64)exp(1.0 / 3 * log(n)), primarr, count);
    //n^(1/2)內的素數數目
    len2 = pi((__int64)sqrt(n), primarr, count);
    //n^(2/3)內的素數數目
    len3 = pi(mark_len - 1, primarr, count);
 
    //乘積個數
    j = mark_len - 2;
    for (i = (__int64)exp(1.0 / 3 * log(n)); i <= (__int64)sqrt(n); i++) {
        if (!mark[i]) {
            while (i * j > n) {
                if (!mark[j])
                    s++;
                j--;
            }
            sum += s;
        }
    }
    free(mark);
    sum = (len2 - len) * len3 - sum;
    sum += (len * (len - 1) - len2 * (len2 - 1)) / 2;
 
    //尤拉函式
    if (m > len)
        m = len;
    for (i = 0; i < m; i++) {
        q *= primarr[i];
        p *= primarr[i] - 1;
    }
    v = (__int64 *)malloc(sizeof(__int64) * q);
    for (i = 0; i < q; i++)
        v[i] = i;
    for (i = 0; i < m; i++)
        for (j = q - 1; j >= 0; j--)
            v[j] -= v[j / primarr[i]];
 
    sum = phi(n, len, m) - sum + len - 1;
    free(primarr);
    free(v);
    return sum;
}
 
int main()
{
    __int64 n;
    __int64 count;
	int h;
    clock_t start, end;
	while(scanf("%I64d", &n)!=EOF)
	{

		p=1;
		q=1;
		start = clock();
		count = prime(n);
		end = clock() - start;
		printf("%I64d(%d億)內的素數個數為%I64d\n",n,n/100000000,count);
		printf("用時%lf毫秒\n",(double)end/1000);
	}
    return 0;
}