seqfun
Yahoo Group

Count primes in a range up to 10^18

===============================================
Cino Hilliard     Message 1 of 1  Sep 23, 2007
-----------------------------------------------
//gcc 4.1.2
//Countprimes uses Achim Flammenkamp's algorithm below. 
//Count the primes from start and stop up to 2^64. 
//Comments below with Cino show some of my additions.
//                       Cino Hilliard
//                       Sept 2007
//      http://groups.yahoo.com/group/seqfun/files/
#include <stdio.h> 
/* Achim Flammenkamp, Bielefeld, D, achim@... */
#include <malloc.h>
/* ifdefs:  LONG=uns64  _ANSI_EXTENSION  true_64bit_words  */
#include <stdlib.h>
/* SUM_{p<x} 1/p= ln(ln(x))+0.5772156649015-0.315718451893 */
#include <math.h>  
/* extern double sqrt(double), log(double), floor(double); */
#include <string.h>
#include <windows.h>
#define  true_64bit_words 0 // Cino This is necessary
#define _ANSI_EXTENSION 0   // So is this
#ifndef  uns32              
/* Version 1.0a 1998-05-19  source availible     */
#define  uns32  unsigned    
/* Version 2.0a 1998-05-23  5 years in the WWW   */
#endif                      
/* Version 2.0b 2003-06-03  0-1-typo fixed Argh! */
#ifndef  uns64              
/* Version 2.0c 2003-06-04  works up to 2^64-1   */
# ifdef   _ANSI_EXTENSION   
/* uses < 1938 MB for minimal sieve size to 2^64 */
#  define  uns64  unsigned long long
#define LONG uns64
# else
#  define  uns64 unsigned long
# endif
#endif
#ifdef   LONG
# ifdef   _ANSI_EXTENSION
#  define  UL  "%llu"
# else
#  define  UL  "%lu"
# endif
# define  HALF  uns32
#else
# define  LONG  uns32
# define  UL  "%u"
# define  HALF  unsigned short
#endif
#define  BITS  8
#define  LCM  (2*3*5)
#define  SIEVE_SIZE  (16<<17)
#define chr64(a) _ui64toa(a,str,10)
#define ld long double
#define digits 40

//Cino add the timer
#define timer GetTickCount()/1000.0
long u,v,ai,j,aa[50],s[50],ln,lna,lnb;
unsigned long long p[50],cr,z,x,y,xx;
static char a1[50],a2[50],str[50];
static FILE*  fp1;
static FILE*  fp2;
LONG anz,pi2;
static long  a[]={1,2,4,8,16,32,64,128};
static long b1[]={1,7,11,13,17,19,23,29};
static long r;
long digits2;
static char *pEnd;
char stop1[30],start1[30];
int paricalc(void);
long double sum,est,xld;
unsigned long long p[50],s1[20],p1[9],x0;
LONG ii2;
LONG  stop, start, stlcm, ii, jj, hh;
int main(unsigned argc, char *argv[])
{
float t;
LONG size, hi, h, i, ji, js;
uns32  k, hj, c, cj;
LONG j,*index;
unsigned char  *sieve, *prime_diff;
unsigned char  bi, b, m;
HALF  psize;
int lim;
printf(" Count primes in a range from m to n up to 10^18 \n");
while(1)
{
sum=0;
anz=0;
pi2=0;
digits2=0;
for(u=0;u<digits;u++)s[u]=0;
printf("start  => ");
gets(start1);
start = strtoull(start1,&pEnd,10);
printf("stop   => ");
gets(stop1);
stop = strtoull(stop1,&pEnd,10);
//for(ai=ln;ai>=0;ai--)
//printf("%i",a1[ai]);

size = SIEVE_SIZE;
if (stop < start)  goto finish;
t = timer;
system("date/t");
system("time/t");
j =(LONG)sqrt((long double)stop+0.1);
if (!(j&1)) j--;
if ((j-1)/2 >= (5-1)/2*(LCM/5)*size)  
/* to prevent overflow if j == 2^32-1 */
size = (j/48+1)*2;
if (!(sieve = (unsigned char*)calloc(size, sizeof(unsigned char))))
return  fprintf(stderr,"error: can't get %u kB storage for sieve\n",
((unsigned)sizeof(unsigned char)*size+(1<<10)-1)>>10), 2;
psize = (int)((1.015*j)/(log((double)j)-1)); 
/* upper bound for Pi(stop^1/2) */
if (!(index = (LONG*) malloc((sizeof(LONG)+1)*psize)))
return fprintf(stderr,"error: can't get storage for %u 
primes\n",psize), 4;
//printf("Memory Used %d\n",index);
prime_diff = (unsigned char*)(index+psize);
if (start <= 2){anz++;} //{} are necessary!
if (start <= 3){anz++;}
if (start%2 == 0)start += 1;
if (start%3 == 0)start += 2;
stlcm = start/LCM;
hh = size*(stlcm/size);
/*** sifting the sieve primes ***/
k=0;
m=0;
h= i= cj= 0;
js=jj=0;
b=1;
c=1;
for(;;)
{  switch (c)
{  do
{  case 1:
i++;
if ((m+=1) >= LCM/3)  m -= LCM/3;
jj += h;
h++;
jj += h;
if (!(sieve[i]&b))
{  c= 2;
break;
}
case 2:
i++;
if (i == size)
{  i=0;
cj += size;
b += b;
if (!b)  break;
}
if ((m+=1) >= LCM/3)  m -= LCM/3;
jj += h;
if (!(sieve[i]&b))
{  c= 1;
break;
}
} while (1);
}
if (8*jj > stop/3)  break;
j = 3*(cj+i)+c;
bi= m - !!m - m/BITS;
prime_diff[k]= BITS*2*(j/LCM-js); /* difference of successive 
primes < 480 */
prime_diff[k] |= bi;
js= j/LCM;
index[k]= ((bi&1) != (bi>>2&1) ? 5 : 0);
ii= (8*jj)/(LCM/3);
if (ii < stlcm)
ii += j*((stlcm-ii)/j);
if (ii < stlcm)
{  hi = (index[k] ? 19 : 1);  
/* inverse zu index[k]  ---  0 <--> 1 Argh!! */
hj= js+js;
ji= 2*(3*m+c);
if (ji >= LCM)  ji-=LCM,  hj++;
switch (c)
{  do
{  case 1:  ii += 2*hj;
hi += 2*ji;
while (hi >= LCM)  hi -= LCM,  ii++;
if (ii >= stlcm  &&  hi!=5 && hi!=25)  break;
case 2:  ii += hj;
hi += ji;
while (hi >= LCM)  hi -= LCM,  ii++;
if (ii >= stlcm  &&  hi!=5 && hi!=25)  break;
}
while(1);
}
index[k] = (BITS*hi)/LCM;
}
index[k] |= BITS*(ii-hh);
k++;
if (jj >= size)
continue;
/*** sift with prime pre-sieve ***/
#ifdef  true_64bit_words
ii = 8*jj;
#define  hi ii
#else
hi = 8*jj;
#endif
ji= 2*(cj+i)+1;
hj= 2*(ji+c)-3;
bi= 1;
switch(m)
{   do
{  case 0:   while(hi >= size)
{  hi -= size;  bi+=bi;
if (!bi)  goto go_on;
}
sieve[hi] |= bi;  hi += 2*j;
case 2:   while(hi >= size)
{  hi -= size;  bi+=bi;
if (!bi)  goto go_on;
}
sieve[hi] |= bi;  hi += hj;
case 3:   while(hi >= size)
{  hi -= size;  bi+=bi;
if (!bi)  goto go_on;
}
sieve[hi] |= bi;  hi += ji;
case 4:   while(hi >= size)
{  hi -= size;  bi+=bi;
if (!bi)  goto go_on;
}
sieve[hi] |= bi;  hi += hj;
case 5:   while(hi >= size)
{  hi -= size;  bi+=bi;
if (!bi)  goto go_on;
}
sieve[hi] |= bi;  hi += ji;
case 6:   while(hi >= size)
{  hi -= size;  bi+=bi;
if (!bi)  goto go_on;
}
sieve[hi] |= bi;  hi += hj;
case 7:   while(hi >= size)
{  hi -= size;  bi+=bi;
if (!bi)  goto go_on;
}
sieve[hi] |= bi;  hi += 2*j;
case 9:   while(hi >= size)
{  hi -= size;  bi+=bi;
if (!bi)  goto go_on;
}
sieve[hi] |= bi;  hi += ji;
lim= size - LCM/3*j;
while ((int)hi < lim)
{  sieve[hi] |= bi;  hi += 2*j;
sieve[hi] |= bi;  hi += hj;
sieve[hi] |= bi;  hi += ji;
sieve[hi] |= bi;  hi += hj;
sieve[hi] |= bi;  hi += ji;
sieve[hi] |= bi;  hi += hj;
sieve[hi] |= bi;  hi += 2*j;
sieve[hi] |= bi;  hi += ji;
}
} while (1);
do {
case 1:   while(hi >= size)
{  hi -= size;  bi+=bi;
if (!bi)  goto go_on;
}
sieve[hi] |= bi;  hi += ji;
case 8:   while(hi >= size)
{  hi -= size;  bi+=bi;
if (!bi)  goto go_on;
}
sieve[hi] |= bi;  hi += hj;
lim= size - LCM/3* 5;
while ((int)hi < lim)
{  sieve[hi] |= bi;  hi += ji;
sieve[hi] |= bi;  hi += hj;
sieve[hi] |= bi;  hi += ji;
sieve[hi] |= bi;  hi += hj;
sieve[hi] |= bi;  hi += ji;
sieve[hi] |= bi;  hi += hj;
sieve[hi] |= bi;  hi += ji;
sieve[hi] |= bi;  hi += hj;
sieve[hi] |= bi;  hi += ji;
sieve[hi] |= bi;  hi += hj;
}
} while (1);
}
go_on: ;
#ifdef  true_64bit_words
#undef  hi
#endif
}
/****** main sifting starts *****/
for (i=size;i--;)
sieve[i]=0;
sieve[0] |= !hh;  /* 1 is not prime */
if (start <= 5){anz++;}
if (start%5 == 0)  start += 2*(3-start%3);
hh *= LCM;
while (1)
{  j= prime_diff[0]/BITS;
for (h=1;h<k;h++)   
/* sieve with next sieve prime (h=0 is prime 5) */
{  j += prime_diff[h]/BITS;
ii = index[h]/BITS;
if (ii >= size)
{  index[h] -= size*BITS;
continue;
}
hj= (size <= LCM/2*j ? 0 : size - LCM/2*j);
i=ji=js= j;
ji +=js;
i += ji;
#ifdef  true_64bit_words
#define  hi ii
#else
hi = ii;
#endif
//Cino - equated index[h] in place and added break in place of goto 

switch( prime_diff[h]%BITS )
{
case 0:   switch( index[h]%BITS )
{  do {
case 0:   if (hi >= size)  {index[h] = 0;  break;}
sieve[hi] |=  1;  hi += i;
case 1:   if (hi >= size)  {index[h] = 1;  break;}
sieve[hi] |=  2;  hi += ji;
case 2:   if (hi >= size)  {index[h] = 2;  break;}
sieve[hi] |=  4;  hi += js;
case 3:   if (hi >= size)  {index[h] = 3;  break;}
sieve[hi] |=  8;  hi += ji;
case 4:   if (hi >= size)  {index[h] = 4;  break;}
sieve[hi] |= 16;  hi += js;
case 5:   if (hi >= size)  {index[h] = 5;  break;}
sieve[hi] |= 32;  hi += ji;
case 6:   if (hi >= size)  {index[h] = 6;  break;}
sieve[hi] |= 64;  hi += i;
case 7:   if (hi >= size)  {index[h] = 7;  break;}
sieve[hi] |=128;  hi += js+1;
lim= hj-1;
while ((int)hi < lim)
{  sieve[hi] |=  1;  hi += i;
sieve[hi] |=  2;  hi += ji;
sieve[hi] |=  4;  hi += js;
sieve[hi] |=  8;  hi += ji;
sieve[hi] |= 16;  hi += js;
sieve[hi] |= 32;  hi += ji;
sieve[hi] |= 64;  hi += i;
sieve[hi] |=128;  hi += js+1;
}
} while (1);
}
case 1:   js+=1;  i+=1;  ji+=1;
switch( index[h]%BITS )
{  do {
case 5:   if (hi >= size)  {index[h] = 5;  break;}
sieve[hi] |= 32;  hi += ji;
case 4:   if (hi >= size)  {index[h] = 4;  break;}
sieve[hi] |= 16;  hi += js;
case 0:   if (hi >= size)  {index[h] = 0;  break;}
sieve[hi] |=  1;  hi += ji-1;
case 7:   if (hi >= size)  {index[h] = 7;  break;}
sieve[hi] |=128;  hi += js;
case 3:   if (hi >= size)  {index[h] = 3;  break;}
sieve[hi] |=  8;  hi += ji;
case 2:   if (hi >= size)  {index[h] = 2;  break;}
sieve[hi] |=  4;  hi += i;
case 6:   if (hi >= size)  {index[h] = 6;  break;}
sieve[hi] |= 64;  hi += js;
case 1:   if (hi >= size)  {index[h] = 1;  break;}
sieve[hi] |=  2;  hi += i;
lim= hj-7;
while ((int)hi < lim)
{  sieve[hi] |= 32;  hi += ji;
sieve[hi] |= 16;  hi += js;
sieve[hi] |=  1;  hi += ji-1;
sieve[hi] |=128;  hi += js;
sieve[hi] |=  8;  hi += ji;
sieve[hi] |=  4;  hi += i;
sieve[hi] |= 64;  hi += js;
sieve[hi] |=  2;  hi += i;
}
} while (1);
}
case 2:   i+=2;  ji+=2;
switch( index[h]%BITS )
{  do {
case 0:   if (hi >= size)  {index[h] = 0;  break;}
sieve[hi] |=  1;  hi += js;
case 6:   if (hi >= size)  {index[h] = 6;  break;}
sieve[hi] |= 64;  hi += ji;
case 1:   if (hi >= size)  {index[h] = 1;  break;}
sieve[hi] |=  2;  hi += js;
case 7:   if (hi >= size)  {index[h] = 7;  break;}
sieve[hi] |=128;  hi += ji;
case 3:   if (hi >= size)  {index[h] = 3;  break;}
sieve[hi] |=  8;  hi += i;
case 5:   if (hi >= size)  {index[h] = 5;  break;}
sieve[hi] |= 32;  hi += js+1;
case 2:   if (hi >= size)  {index[h] = 2;  break;}
sieve[hi] |=  4;  hi += i;
case 4:   if (hi >= size)  {index[h] = 4;  break;}
sieve[hi] |= 16;  hi += ji;
lim= hj-11;
while ((int)hi < lim)
{  sieve[hi] |=  1;  hi += js;
sieve[hi] |= 64;  hi += ji;
sieve[hi] |=  2;  hi += js;
sieve[hi] |=128;  hi += ji;
sieve[hi] |=  8;  hi += i;
sieve[hi] |= 32;  hi += js+1;
sieve[hi] |=  4;  hi += i;
sieve[hi] |= 16;  hi += ji;
}
} while (1);
}
case 3:   js+=1;  i+=3;  ji+=1;
switch( index[h]%BITS )
{  do {
case 5:   if (hi >= size)  {index[h] = 5;  break;}
sieve[hi] |= 32;  hi += ji+1;
case 2:   if (hi >= size)  {index[h] = 2;  break;}
sieve[hi] |=  4;  hi += js;
case 1:   if (hi >= size)  {index[h] = 1;  break;}
sieve[hi] |=  2;  hi += ji;
case 7:   if (hi >= size)  {index[h] = 7;  break;}
sieve[hi] |=128;  hi += i;
case 4:   if (hi >= size)  {index[h] = 4;  break;}
sieve[hi] |= 16;  hi += js;
case 3:   if (hi >= size)  {index[h] = 3;  break;}
sieve[hi] |=  8;  hi += i;
case 0:   if (hi >= size)  {index[h] = 0;  break;}
sieve[hi] |=  1;  hi += ji;
case 6:   if (hi >= size)  {index[h] = 6;  break;}
sieve[hi] |= 64;  hi += js;
lim= hj-13;
while ((int)hi < lim)
{  sieve[hi] |= 32;  hi += ji+1;
sieve[hi] |=  4;  hi += js;
sieve[hi] |=  2;  hi += ji;
sieve[hi] |=128;  hi += i;
sieve[hi] |= 16;  hi += js;
sieve[hi] |=  8;  hi += i;
sieve[hi] |=  1;  hi += ji;
sieve[hi] |= 64;  hi += js;
}
} while(1);
}
case 4:   js+=1;  i+=3;  ji+=3;
switch( index[h]%BITS )
{  do {
case 5:   if (hi >= size)  {index[h] = 5;  break;}
sieve[hi] |= 32;  hi += js;
case 6:   if (hi >= size)  {index[h] = 6;  break;}
sieve[hi] |= 64;  hi += ji;
case 0:   if (hi >= size)  {index[h] = 0;  break;}
sieve[hi] |=  1;  hi += i;
case 3:   if (hi >= size)  {index[h] = 3;  break;}
sieve[hi] |=  8;  hi += js;
case 4:   if (hi >= size)  {index[h] = 4;  break;}
sieve[hi] |= 16;  hi += i;
case 7:   if (hi >= size)  {index[h] = 7;  break;}
sieve[hi] |=128;  hi += ji;
case 1:   if (hi >= size)  {index[h] = 1;  break;}
sieve[hi] |=  2;  hi += js;
case 2:   if (hi >= size)  {index[h] = 2;  break;}
sieve[hi] |=  4;  hi += ji-1;
lim= hj-17;
while ((int)hi < lim)
{  sieve[hi] |= 32;  hi += js;
sieve[hi] |= 64;  hi += ji;
sieve[hi] |=  1;  hi += i;
sieve[hi] |=  8;  hi += js;
sieve[hi] |= 16;  hi += i;
sieve[hi] |=128;  hi += ji;
sieve[hi] |=  2;  hi += js;
sieve[hi] |=  4;  hi += ji-1;
}
} while (1);
}
case 5:   js+=2;  i+=4;  ji+=2;
switch( index[h]%BITS )
{  do {
case 0:   if (hi >= size)  {index[h] = 0;  break;}
sieve[hi] |=  1;  hi += ji;
case 4:   if (hi >= size)  {index[h] = 4;  break;}
sieve[hi] |= 16;  hi += i;
case 2:   if (hi >= size)  {index[h] = 2;  break;}
sieve[hi] |=  4;  hi += js-1;
case 5:   if (hi >= size)  {index[h] = 5;  break;}
sieve[hi] |= 32;  hi += i;
case 3:   if (hi >= size)  {index[h] = 3;  break;}
sieve[hi] |=  8;  hi += ji;
case 7:   if (hi >= size)  {index[h] = 7;  break;}
sieve[hi] |=128;  hi += js;
case 1:   if (hi >= size)  {index[h] = 1;  break;}
sieve[hi] |=  2;  hi += ji;
case 6:   if (hi >= size)  {index[h] = 6;  break;}
sieve[hi] |= 64;  hi += js;
lim= hj-19;
while ((int)hi < lim)
{  sieve[hi] |=  1;  hi += ji;
sieve[hi] |= 16;  hi += i;
sieve[hi] |=  4;  hi += js-1;
sieve[hi] |= 32;  hi += i;
sieve[hi] |=  8;  hi += ji;
sieve[hi] |=128;  hi += js;
sieve[hi] |=  2;  hi += ji;
sieve[hi] |= 64;  hi += js;
}
} while (1);
}
case 6:   js+=1;  i+=5;  ji+=3;
switch( index[h]%BITS )
{  do {
case 5:   if (hi >= size)  {index[h] = 5;  break;}
sieve[hi] |= 32;  hi += i;
case 1:   if (hi >= size)  {index[h] = 1;  break;}
sieve[hi] |=  2;  hi += js;
case 6:   if (hi >= size)  {index[h] = 6;  break;}
sieve[hi] |= 64;  hi += i;
case 2:   if (hi >= size)  {index[h] = 2;  break;}
sieve[hi] |=  4;  hi += ji;
case 3:   if (hi >= size)  {index[h] = 3;  break;}
sieve[hi] |=  8;  hi += js;
case 7:   if (hi >= size)  {index[h] = 7;  break;}
sieve[hi] |=128;  hi += ji+1;
case 0:   if (hi >= size)  {index[h] = 0;  break;}
sieve[hi] |=  1;  hi += js;
case 4:   if (hi >= size)  {index[h] = 4;  break;}
sieve[hi] |= 16;  hi += ji;
lim= hj-23;
while ((int)hi < lim)
{  sieve[hi] |= 32;  hi += i;
sieve[hi] |=  2;  hi += js;
sieve[hi] |= 64;  hi += i;
sieve[hi] |=  4;  hi += ji;
sieve[hi] |=  8;  hi += js;
sieve[hi] |=128;  hi += ji+1;
sieve[hi] |=  1;  hi += js;
sieve[hi] |= 16;  hi += ji;
}
} while (1);
}
case 7:   js+=2;  i+=6;  ji+=4;
switch( index[h]%BITS )
{  do {
case 0:   if (hi >= size)  {index[h] = 0;  break;}
sieve[hi] |=  1;  hi += js-1;
case 7:   if (hi >= size)  {index[h] = 7;  break;}
sieve[hi] |=128;  hi += i;
case 6:   if (hi >= size)  {index[h] = 6;  break;}
sieve[hi] |= 64;  hi += ji;
case 5:   if (hi >= size)  {index[h] = 5;  break;}
sieve[hi] |= 32;  hi += js;
case 4:   if (hi >= size)  {index[h] = 4;  break;}
sieve[hi] |= 16;  hi += ji;
case 3:   if (hi >= size)  {index[h] = 3;  break;}
sieve[hi] |=  8;  hi += js;
case 2:   if (hi >= size)  {index[h] = 2;  break;}
sieve[hi] |=  4;  hi += ji;
case 1:   if (hi >= size)  {index[h] = 1;  break;}
sieve[hi] |=  2;  hi += i;
lim= hj-29;
while ((int)hi < lim)
{  sieve[hi] |=  1;  hi += js-1;
sieve[hi] |=128;  hi += i;
sieve[hi] |= 64;  hi += ji;
sieve[hi] |= 32;  hi += js;
sieve[hi] |= 16;  hi += ji;
sieve[hi] |=  8;  hi += js;
sieve[hi] |=  4;  hi += ji;
sieve[hi] |=  2;  hi += i;
}
} while (1);
}
hi -= size;
index[h] |= BITS*hi;
}
}

/* output of remaining (prime) numbers 
cino version with no switch error*/
char str[24],str1[24];
i=(start<=hh ? 0 : (start-hh)/LCM);
hh += LCM*i;
bi= sieve[i];
for (;i<size;i++)
{
bi= sieve[i];
//This gets the next 8 numbers 30k+r, r=1,7,11,13,17,19,23,29
for(r=0;r<8;r++)
{
if(!(bi&a[r])) //test for 0 bit
{
ii=hh+b1[r]; //track the numbers in the wheel
if (ii > stop)goto end;
if(ii >= start&&ii<=stop) //absolutely necessary!
{
anz++;
}
}
}
hh += LCM;
sieve[i]=0;
}
}
end:
free(index);
free(sieve);
finish:
t=timer-t;
digits2=digits;
printf("# {%s <= Primes <= %s} = %s\n",start1,stop1,chr64(anz));
printf("Sec = %.06f\n",t);
}
return  anz;
}
===============================================
Cached by Georg Fischer at Nov 14 2019 12:47 with clean_yahoo.pl V1.4