Appriximating the nth prime number and nth lower twin prime number
By Cino Hilliard
Part 1
A very good approximation of prime(x), the xth prime number, using the formula Pi(prime(x)) = x.
We choose brackets r1 and r2 such that r1 < prime(x) <= r2. we then perform trial and error
using the best approximation for Pi(x) ~ Rg(x) or the Riemann Gram approximation for
Pi(x). Since Pi(prime(n)) = n, we seek an estimates of prime(n) such that Pi(prime(n))
converges to n. This turns out to be est = Rg(r) where r is the average of the previous
bounds, r1 and r2. Then if est <= n, we let the lower bound r1 = r orherwise we let
the upper bound r2 = r. Thus if we start with appropriate bounds r1 and r2, we will get a
convergence Rg(r) --> n. This desired result implies that r ~ prime(n). The Pari script below
uses bounds r1 = n*log(n) and r2 = r1*2.
*************************************Pari script *********************************
primex(n) =
{
local(r1,r2,r,est);
if(n==1,return(2));
r1 = n*log(n);
r2 = r1*2;
for(x=1,50,
r=(r1+r2)/2.;
est = Rg(r); /* Use Riemann Gram estimator of Pi(n) */
print(r);
if(est <= n,r1=r,r2=r);
);
r;
}
Rg(x) =
{
local(n=1,L,s=1,r);
L=r=log(x);
while(s<10^30*r,
s=s+r/(zeta(n+1)*n);
n=n+1;
r=r*L/n);
(s)
}
This is kinda like using a calculator to compute r = sqrt(2). We know 1 < r < 2. So we
set our lower and upper bounds r1=1 and r2 = 2 and go half way or take the average
1. and test the result 2 as is shown below.
1. r = (r1+r2)/2 = 1.5.
2. r^2 = 1.5^2 = 2.25 > 2.
This tells us that the root r must be less than 1.5. So we set the upper bound r2 to 1.5
and repeat 1. and 2.
r = (1+1.5)/2 = 1.25
r^2 = 1.5625
This tells us that the root r must be greater than 1 and less than 1.5. So we set the lower
bound r1 to 1.25 and repeat 1. and 2.
r = (1.25+1.5)/2 = 1.375
r^2 =1.890625
Bacically the rule is simple. If the test result r^2 > 2 then the lower upper bound r2 = r
orherwise the lower bound r1 = r.
{ r1 if r^2<2
Continuing this process using the rule r = {
{r2 if r^2>2,r2
r = (1.375+1.5)/2 = 1.4375
r^2= 2.06640625
r = (1.375+1.4375)/2 = 1.40625
r^2 = 1.977539063
r = (1.40625+1.4375)/2 = 1.421875
r^2 = 2.021728516
r = (1.40625+1.421875)/2 = 1.4140625
r^2 = 1.999572754
So 7 iterations gives sqrt(2) accurate to 3 places.
Part 2
A very good approximation of twin(x), the xth twin prime, using the formula Pi2(twin(x)) = x.
We choose brackets or bounds r1 and r2 such that r1 < twin(x) <= r2. we then perform trial and error
using the best approximation for Pi2(x) ~ Li_2(x). Then twinx(n) approximates the nth lower twin prime number
to an accuracy of log10(n)/2 places. This requires the Li_2() estimator for Pi2(x) or the count of twins < x
twinx(n) =
{
local(r1,r2,r,est);
if(n==1,return(2));
r1 = n;
r2 = n*n/log(n);
for(x=1,100,
r=(r1+r2)/2.;
est = Li_2(r); /* Use Riemann Gram estimator of Pi(n) */
if(est <= n,r1=r,r2=r);
);
r;
}
****************************Gcc+ Gmp code to count and sum twin primes *******************
//SumtwinsFirstN uses Achim Flammenkamp's algorithm below. Sum first n lower twin
//primes from start and stop up to 2^64. The routine uses gmp to add the primes
//Also 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
// Oct 2006 - Oct 2008
//Feb 2007 compile with c:\cygwin\bin\g++ but must use sprintf. ulltoa not allowed
#include <stdio.h> /* Achim Flammenkamp, Bielefeld, D, achim@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 "gmp.h" /* extern double sqrt(double), log(double), floor(double); */
#include <string.h>
#include <windows.h>
#define true_64bit_words 1 //Cino This is necessary
#define uns64 unsigned long long
#define uns32 unsigned long
#define LONG uns64
#define UL "%llu"
#define HALF uns32
#define BITS 8
#define LCM (2*3*5)
#define SIEVE_SIZE (16<<18)
#define fact .2565
#define chr64(a) ulltoa(a,str,10)
//Cino save partial sums using gmp.Use str to over come gmp no long long
#define use(x) do {\
anz++;\
mpz_set_str(x2,chr64(x),10);\
mpz_add(sum,sum,x2);\
} while(0);
//Cino add the timer
#define timer GetTickCount()/1000.0
unsigned long long j,s0,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,x1,ct,sq;
static FILE* fp1;
static FILE* fp2;
static FILE* fp3;
static char *pEnd;
char buffer[200];
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;
static long r1[]={0,3,5};
char str[24],str1[24];
char *binary(long);
char start1[20],stop1[20];
unsigned long long anz,pi2,start,stop,stop2,len;
int main(unsigned argc, char *argv[])
{
FILE *fp1;
float t;
uns64 k, hj, c, cj;
unsigned char *sieve, *prime_diff;
unsigned char bi, b, m;
unsigned long long stlcm, *index,ii,ii2,ii3,jj, hh;
unsigned long long size, hi, h, i, ji, js;
HALF psize;
int lim,ln;
mpz_t sum,x2;
mpz_init(sum);
printf(" Sieve, sum of lower twin primes in a range beg from start max = 1 to 10^18\n");
printf(" Cino Hilliard\n");
printf(" Oct, 2008\n");
printf(" Input start and Range as prompted below\n");
printf("\n");
while(1)
{
j=0;s0=0;s1=0;s2=0;s3=0;s4=0;s5=0;s6=0;s7=0;s8=0;s9=0;s10=0;x1=0;
anz=0;pi2=0;*start1=0;*stop1=0;start=0;stop=0;size=0;ct=1;
mpz_set_ui(sum,0LL);
mpz_init(x2);
printf("start => ");
gets(start1);
start = strtoull(start1,&pEnd,10);
printf("Range => ");
gets(stop1);
stop2 = strtoull(stop1,&pEnd,10);
stop = 1000000000000000LL;
t = timer;
j =(LONG)sqrt((long double)stop+0.1)+10;
//printf("Sqrt stop = %llu\n",j);
if (!(j&1)) j--;
size = (j/48+1)*2;
//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 %llu\n",index);
prime_diff = (unsigned char*)(index+psize);
if (start <= 3) use(3LL);
if (start%2 == 0) start += 1;
if (start%3 == 0) start += 2;
stlcm = start/LCM;
hh = size*(stlcm/size);
/*** sifting the sieve primes ***/
ii2=3;
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(5LL);
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 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. Modified by Cino Loop ***/
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(anz==stop2)goto end;
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) //absolutely necessary!
{
if(ii2+2==ii) //If it is a twin, add it.
{
use(ii2);
ii3=ii2; //save the last twin
sq = round(sqrtf(ii2-1));
if(sq*sq==ii2-1)
{
ct++;
printf("%s ",chr64(ct));
printf("%s\n",chr64(ii2));
}
}
ii2=ii;
}
}
}
hh += LCM;
sieve[i]=0;
}
}
end:
printf("Last lower twin prime %s \n",chr64(ii3));
printf("Start %s \n",chr64(start));
printf("Stop %s \n",chr64(start+stop2-1));
//printf("%d %s\n",r,chr64(ii));
free(index);
free(sieve);
finish:
t=timer-t;
//Cino - Estimate of sum of primes in an interval start stop
//printf("%s\n","Sum estimate = stop2^2/(2*log(stop2)-1) - start^2/(2*log(start)-1)");
printf("Sum of lower twin primes ");
mpz_out_str(stdout,10,sum);
printf("\n");
printf("Sum of both twin primes ");
mpz_add(sum,sum,sum);
mpz_add_ui(sum,sum,2*anz);
mpz_out_str(stdout,10,sum);
printf("\n");
printf("Sec = %.6f\n",t);
}//wend(1)
return anz;
}
sumprimesGmp.c
|
