本程序改编自intfree写的pascal程序。

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

__int64 maxn = 1000000000L;
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～255的每一个数中位为1的位的个数，
 // 例如7有3位为1，则bits[7] = 3
 bits[0] = 0;
 for(k=1; k<=255; ++k)
 {
  bits[k] = bits[k>>1] + (k & 1); 
 }
}

// 算法思想如下：
// 以2*3*5*7*11*13（即30030）为基，把整个范围（假设为1000000）分成30030段，
// 每段数据为：
// 第0段：0, 30030, 60060, 90090, ..., 990990 // 本段显然可以排除
// 第1段：1, 30031, 60061, 90091, ..., 990991
// 第2段：2, 30032, 60062, 90092, ... // 本段也可以排除（但丢掉了2）
// 第3段：3, 30033, 60063, 90093, ...  // 本段也可以排除（但丢掉了3）
// 第4段：4, 30034, 60064, 90094, ...  // 本段也可以排除
// 第5段：5, 30035, 60065, 90095, ...  // 本段也可以排除（但丢掉了5）
// ......
// 第13段：13, 30043, 60073, 90103, ...  // 本段也可以排除（但丢掉了13）
// ......
// 第16段：16, 30046, 60076, 90106, ...  // 本段也可以排除
// 第17段：17, 30047, 60077, 90107, ...  // 本段需要考察（但排除17）
// ......
// 第9010段：9010, 39040, 69070, ..., 969970, 1000000 //
// 第9011段：9011, 39041, 69071, ..., 969971
// ......
// 第30029段：30029, 60059, 90089, ..., 990989
// 从上面分析可以看出，丢掉的素数正好是范围的平方根范围内的素数
// 所以需要加上这个素数个数
 
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];

 // 找到13以后的素数的最小的倍数的偏移位置
 // 以17为例，加上最小的倍数的偏移为2（此时数据为60061），则偏移19时肯定还是合数
 // 即60061＋30030*17还是17的倍数，即为合数。
 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);

 if(maxn > 3500000000L)
  m = 6;
 else if(maxn > 10000000)
  m = 5;
 else if(maxn > 50000)
  m = 4;
 else if(maxn > 1000)
  m = 3;
 else if(maxn > 40)
  m = 2;
 else
  m = 1;

 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;
}

评论