高效率的算法 计算100亿内的素数的个数 3.01G的cpu运算只用8.64秒！！！！

 计算100亿内的素数的个数 3.01G的cpu运算只用8.64秒！！！！

#include <windows.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <time.h>

__int64 maxn = 10000000L;
const long maxlen = 100000;
long m = 5;
const long dots = 80;


BYTE     *mask;
long     primes[maxlen+2], mprimes[maxlen+2], ans[maxlen+2];
long     bits[256];
long     len, r; 

  
void solve(long a, long b, long c, long &x, long &y)
{
    if(a == 0)
    {
        x = 0; 
        y = c/b;
    }
    else
    {
        solve(b%a, a, c, y, x);
        x -= b/a*y;
    }
}

// primes数组从下标1开始存放素数，例如primes[1]为2
void init()
{
    long         p, k;
    
    k = 2;
    len = 0;
    
    // 求待求范围的平凡根范围内的素数
    do
    {
        ++len;
        primes[len] = k;
        
        do
        {
            ++k;
            p = 1;
            while((primes[p]*primes[p] < k) && (k%primes[p] != 0))
            {
                ++p;
            }
        }while(k%primes[p] == 0);
    }while(1.0*k*k <= maxn);
    
    bits[0] = 0;
    for(k=1; k<=255; ++k)
    {
        bits[k] = bits[k>>1] + (k & 1);    
    }
}

void makeprimes()
{
    long i, k, n, size, kmax, min;
    long last, p, x, dot;

    // 以2*3*5*7*11*13为基，循环完后n为30030
    n = 1;
    for(i=1; i<=m; ++i)
        n *= primes[i];

    for(i=m+1; i<=len; ++i)
    {
        p = primes[i];
        solve(p, -n, 1, x, mprimes[i]);
        mprimes[i]  = (mprimes[i]%p+p)%p;
    }
    
    r = len;
    size = long(maxn/n+7)>>3;
    last = 0;
    dot = 0;
    
    memset(ans, 0, sizeof(ans));
    mask = new BYTE[size];
    assert(mask);

    for(i=1; i<n; ++i)
    {
        k = 1;
        while(k < m && (i%primes[k] != 0))
            ++k;
        if(i%primes[k] == 0)    // 最小值即为2，3，5，7，11，13的倍数的段可以直接排除
            continue;

        p = i-last;
        last = i;
        kmax = (long)((maxn-i)/n);
        for(k=m+1; k<=len; ++k)
            ans[k] = (ans[k] + mprimes[k]*p)%primes[k];

        // 重置mask数组为0
        memset(mask, 0, sizeof(BYTE)*size);
        for(k=m+1; k<=len; ++k)
        {
            p = primes[k];
            if(ans[k] == 0)
                min = p-1;
            else
                min = ans[k]-1;

//            if((i>(maxn-maxn/n*n) && (min > maxn/n-2)) || (i<=(maxn-maxn/n*n) && (min > maxn/n-1)))
//            {
//                printf("min>size*8: i=%d, k = %d, min=%d, size=%d,prime=%d\n", i, k, min,size, primes[k]);
//                continue;
//            }

            // 下面的汇编代码用来给所有的合数设置标志位
            __asm
            {
                mov eax, mask
                mov ebx, kmax
                mov ecx, min
                mov edx, p
            $001:
                bts [eax], ecx    // 置mask数组的第min位为1
                add ecx, edx    // min位置指向下一个合数
                cmp ecx, ebx    // 是否到了kmax位置
                jl $001            // 不到则转$001位置
            }
        }

        // 去掉本段内所有的合数的个数
        for(k=0; k<size; ++k)
        {
            assert(k <= size);
            kmax -= bits[mask[k]];  
        }
        // 加上本段内的素数个数
        r += kmax;
        if(i/n > dot/dots)
        {
            ++dot;
            putchar('.');
        }
    }
    
    delete []mask;
    
    putchar('\n');    
}

    
int main()
{
    long    start, end, t;

    printf("Input a range: ");
    scanf("%I64d", &maxn);


    start = clock();
    init();
    end = clock();
    printf("Init time=%lf\n", (double)(end-start)/1000);

    start = clock();
    makeprimes();
    end = clock();

    t = end-start;
    printf("sum=%ld, cacl time=%lf\n", r, (double)t/1000);
    
    return 0;
}