Re^2: Mean of 10^9 consecutive primes?
cino hilliard
hillcino368 at hotmail.com
Fri Sep 29 20:38:11 CEST 2006
Hi,
>From: Richard Mathar <mathar at strw.leidenuniv.nl>
>To: seqfan at ext.jussieu.fr
>Subject: Re^2: Mean of 10^9 consecutive primes?
>Date: Thu, 28 Sep 2006 16:28:29 +0200
>
>zs> From zakseidov at yahoo.com Thu Sep 28 11:45:15 2006
>zs> Date: Thu, 28 Sep 2006 02:45:11 -0700 (PDT)
>zs> Subject: Re: Mean of 10^9 consecutive primes?
>zs> To: Richard Mathar <mathar at strw.leidenuniv.nl>, seqfan at ext.jussieu.fr
>zs>
>zs> Right now I'm trying to calculate
>zs> sum of first 10^9 primes... using Mathematica... spent
>zs> already hours...
>zs> n=106333207,s=112600741085770645
>zs> (n may be correct up to =/- 1 because
>zs> interrupting Mathematica may cause it(?))
>
>The sum of primes starting at p=2 =prime(1), accumulating 106333207 terms
>as if to compute A007504(106333207),
>we get the following PARI sequence of partial sums:
>
>#primes sum=A007504
>106333207 112600743259812234
>
>This is based on PARI's nextprime(), which might smuggle
>some pseudo-primes in, the probability of which I cannot estimate
>right now. So far the results are compatible.
>
The 106333207-th prime is 2174041589. Using the following excellent memory
efficient sieve
program sumprimes.c based on Achim Flammenkamp' program. I modified it
somewhat to sum
the primes and to call Pari for higi precision arithmetic to get the correct
output.
f:\eratosieve>sumprimes 2174041589
Sqrt stop = 46626
Extending sieve size to 2097152
Memory Used 6881352
# {(null) <= Primes <= 2174041589} = 106333207
Sum estimate = stop^2/(2*log(stop)-1) - start^2/(2*log(start)-1)
Sum of primes = 112600743259812234
Sum Estimate = 112535470144235294
Relative Error = 0.00057968636518090921
log(SumPrimes/Pi(x)) = 20.78053652748839251854305866
Sec = 64.296249
This agrees with your result.
****A very effecient Gcc sieve program based on Achim Flammenkamp org
program******
//Sumprimes uses Achim Flammenkamp's algorithm below. Sum the primes from
//start and stop up to 2^64. The routine splits the sums into smaller parts
// for a Pari call to handle big numbers whose sum > 10^16. In addition the
// the estimated sum = stop^2/(2*log(stop)-1)- start^2/(2*log(start)-1) is
//computed for comparison. Comments below with Cino show some of my
additions.
// Cino Hilliard
// July 2006
// http://groups.yahoo.com/group/seqfun/files/
#include <stdio.h> /* Achim Flammenkamp, Bielefeld, D,
achim at uni-bielefeld.de */
#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)
//Cino save partial sums to be added by Pari
#define use(x) do {\
anz++;\
x1=x;\
s1+=x1/10000000;\
x1= x1%10000000;\
s2+=x1/1000000;\
x1= x1%1000000;\
s3+=x1/100000;\
x1= x1%100000;\
s4+=x1/10000;\
x1= x1%10000;\
s5+=x1/1000;\
s6+=x1%1000;\
} while(0);
//Cino add the timer
#define timer GetTickCount()/1000.0
unsigned long long s1=0,s2=0,s3=0,s4=0,s5=0,s6=0,x1;
static FILE* fp1;
static FILE* fp2;
static FILE* fp3;
LONG anz,pi2;
int main(unsigned argc, char *argv[])
{
LONG anz,pi2;
float t;
LONG register size, hi, h, i, ji, js;
uns32 k, hj, c, cj;
LONG j;
unsigned char *sieve, *prime_diff;
unsigned char bi, b, m;
LONG stop, start, stlcm, *index, ii, jj, hh;
HALF psize;
int lim,ln;
if (argc <= 1)
return fprintf(stderr,"usage: %s stop [start [size] ]\nversion
2.0c(%u:%u)"
"
2003-06-04\n",argv[0],(unsigned)sizeof(LONG)*BITS,LCM), 0;
if (argc <= 3 || 1!=sscanf(argv[3],"%u",&k) || !(size= k&~(uns32)1))
size = SIEVE_SIZE;
if (argc <= 2 || 1!=sscanf(argv[2],UL,&start))
start= 0;
if (argc <= 1 || 1!=sscanf(argv[1],UL,&stop))
stop= 1;
anz=0;
pi2=0;
stop = atoll(argv[1]); //Cino- Necessary for it to work
start = atoll(argv[2]);
if (stop < start) goto finish;
t = timer;
j =(LONG)sqrt((long double)stop+0.1);
printf("Sqrt stop = %llu\n",j);
if (!(j&1)) j--;
if ((j-1)/2 >= (5-1)/2*(LCM/5)*size) /* to prevent overflow if j == 2^32-1
*/
// return fprintf(stderr,"error: sieve size must be at least %u\n",
// 1+(j-1)/(5-1)/(LCM/5)+!((j-1)/(5-1)/(LCM/5)&1)), 1;
// size = (j/48+1)*2;
size = 1+(j-1)/(5-1)/(LCM/5)+!((j-1)/(5-1)/(LCM/5)&1);
printf("Extending sieve size to %d\n",size);
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) use(2);
if (start <= 3) use(3);
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) use(5);
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 ddded break in place of goto
out0-out7.
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 ***/
char str[24],str1[24];
i=(start<=hh ? 0 : (start-hh)/LCM);
hh += LCM*i;
bi= sieve[i];
switch (start<=hh ? 0 : (BITS*(unsigned)(start-hh))/LCM)
{
for (;i<size;i++)
{ bi= sieve[i];
case 0:
if (!(bi&1))
{ ii= hh+1;
if (ii > stop) goto end;
use(ii);
}
case 1:
if (!(bi&2))
{ ii= hh+7;
if (ii > stop) goto end;
use(ii);
}
case 2:
if (!(bi&4))
{ ii= hh+11;
if (ii > stop) goto end;
use(ii);
}
case 3:
if (!(bi&8))
{ ii= hh+13;
if (ii > stop) goto end;
use(ii);
}
case 4:
if (!(bi&16))
{ ii= hh+17;
if (ii > stop) goto end;
use(ii);
}
case 5:
if (!(bi&32))
{ ii= hh+19;
if (ii > stop) goto end;
use(ii);
}
case 6:
if (!(bi&64))
{ ii= hh+23;
if (ii > stop) goto end;
use(ii);
}
case 7:
if (!(bi&128))
{ ii= hh+29;
if (ii > stop) goto end;
use(ii);
}
hh += LCM;
sieve[i]=0;
}
}
}
end:
free(index);
free(sieve);
finish:
t=timer-t;
//Cino - Estimate of sum of primes in an interval start stop
//Sum = stop^2/(2*log(stop)-1)- start^2/(2*log(start)-1)
printf("# {%s <= Primes <= %s} = %.18G\n",argv[2],argv[1],(double)anz);
printf("%s\n","Sum estimate = stop^2/(2*log(stop)-1) -
start^2/(2*log(start)-1)");
char sump[200],gp[100],logdata[100];
//Cino Here we combine thepartial sums and call Pari to get the total.
//This assumes gp.exe is in the path.
sprintf(sump,"%.18G0000000+%.18G000000+%.18G00000+%.18G0000+%.18G000+%.18G",
(double)s1,(double)s2,(double)s3,(double)s4,(double)s5,(double)s6);
// printf("%s\n %s\n" ,"String to fgp.gp for Pari to evaluate. ",sump);
if((fp1=fopen("fgp.gp","w"))==0)
{
fprintf(stderr,"Can't open file %s\n","fgp.gp");exit(1);
}
fprintf(fp1,"%s\n",sump);
if(fp1)
{
fflush(fp1);
fclose(fp1);
}
strcpy((char*)gp,"gp.exe -f -q < fgp.gp > pdata.txt");
system(gp);
if((fp2=fopen("pdata.txt","r"))==0)
{
fprintf(stderr,"Can't open file %s\n","pdata.txt");exit(1);
}
sump[0]=0;
fgets(sump,50,fp2);
ln = strlen(sump);
fseek(fp2,0,0);
fgets(sump,ln,fp2);
if(fp2)
{
fflush(fp2);
fclose(fp2);
}
printf("%s%s\n","Sum of primes = ",sump);
if((fp1=fopen("fgp.gp","w"))==0)
{
fprintf(stderr,"Can't open file %s\n","fgp.gp");exit(1);
}
sprintf(logdata,"round((%.18G/%.18G - %.18g/%.18G))\n",
(double)stop*stop,(double)(2*log(stop)-1),(double)start*start,(double)(2*log(start)-1));
fprintf(fp1,"%s\n",logdata);
if(fp1)
{
fflush(fp1);
fclose(fp1);
}
strcpy((char*)gp,"gp.exe -f -q < fgp.gp > pdata.txt");
system(gp);
if((fp2=fopen("pdata.txt","r"))==0)
{
fprintf(stderr,"Can't open file %s\n","pdata.txt");exit(1);
}
logdata[0]=0;
fgets(logdata,50,fp2);
ln = strlen(logdata);
fseek(fp2,0,0);
fgets(logdata,ln,fp2);
if(fp2)
{
fflush(fp2);
fclose(fp2);
}
char * pEnd;
double logd,sumd;
printf("Sum Estimate = %s%s\n",logdata);
logd = strtod(logdata,&pEnd);
sumd = strtod(sump,&pEnd);
printf("Relative Error = %.18G\n",1-logd/sumd);
// output the log(sum/pi(x))
if((fp1=fopen("fgp.gp","w"))==0)
{
fprintf(stderr,"Can't open file %s\n","fgp.gp");exit(1);
}
sprintf(logdata,"log(%s/%.18G)\n",sump,(double)anz);
fprintf(fp1,"%s\n",logdata);
if(fp1)
{
fflush(fp1);
fclose(fp1);
}
strcpy((char*)gp,"gp.exe -f -q < fgp.gp > pdata.txt");
system(gp);
if((fp2=fopen("pdata.txt","r"))==0)
{
fprintf(stderr,"Can't open file %s\n","pdata.txt");exit(1);
}
logdata[0]=0;
fgets(logdata,50,fp2);
ln = strlen(logdata);
fseek(fp2,0,0);
fgets(logdata,ln,fp2);
if(fp2)
{
fflush(fp2);
fclose(fp2);
}
printf("%s%s\n","log(SumPrimes/Pi(x)) = ",logdata);
printf("Sec = %.6f\n",t);
return anz;
}
Compiled as follows with gc2.bat file
g++ %1.c -o %1.exe -O3 -s -mtune=pentium4
g++ -dumpversion
compilation output in the superb ConText editor.
>Executing: c:\context\ConExec.exe "gc2.bat" sumprimes
F:\eratosieve>g++ sumprimes.c -o sumprimes.exe -O3 -s -mtune=pentium4
sumprimes.c: In function `int main(unsigned int, char**)':
sumprimes.c:573: warning: unreachable code at beginning of switch statement
F:\eratosieve>g++ -dumpversion
3.4.4
>Execution finished.
Maybe someone can improve upon the warning.
Here is some output for the 10^n-th prime with counts,sums,estimates, and
timing. I will
take about 11 days to compute the sum of the 10^12 primes on 2.53 GHZ win
xp pro.
The advantage of the program is that input can be in ranges. Eg.,
f:\eratosieve>sumprimes 1000000100000000 1000000000000000
Sqrt stop = 31622778
Extending sieve size to 2097152
Memory Used 6881312
# {1000000000000000 <= Primes <= 1000000100000000} = 2893937
Sum estimate = stop^2/(2*log(stop)-1) - start^2/(2*log(start)-1)
Sum of primes = 2893937144655970518295
Sum Estimate = 2894671936480109822287
Relative Error = -0.00025390731982422032
log(SumPrimes/Pi(x)) = 34.53877644489655810508907688
Sec = 6.999750
Here we have determined the sum of the 2893937 primes between 10^15 and
10^15+10^6. The Relative error is an indicator of the accuracy of this
count.
*********************Some output to estimate timing****************
f:\eratosieve>sumprimes 22801763489
Sqrt stop = 151002
Extending sieve size to 2097152
Memory Used 6881352
# {(null) <= Primes <= 22801763489} = 1000000000
Sum estimate = stop^2/(2*log(stop)-1) - start^2/(2*log(start)-1)
Sum of primes = 11138479445180240497
Sum Estimate = 11133150082192420760
Relative Error = 0.00047846413992580494
log(SumPrimes/Pi(x)) = 23.13367156708259542967846274
Sec = 676.515015
f:\eratosieve>sumprimes 252097800623
Sqrt stop = 502093
Extending sieve size to 2097152
Memory Used 6881352
# {(null) <= Primes <= 252097800623} = 10000000000
Sum estimate = stop^2/(2*log(stop)-1) - start^2/(2*log(start)-1)
Sum of primes = 1234379338586942892505
Sum Estimate = 1233896951484385578900
Relative Error = 0.00039079324116808589
log(SumPrimes/Pi(x)) = 25.53900430684028721415942634
Sec = 7890.264648
Currently I am doing primes < 2760727302517 or the first 100 billion primes.
This should
take a day. My guess the sum of the first 1 trillion primes will take 11
days.
Have fun,
Cino
More information about the SeqFan
mailing list