Some Data and discussion on summing twin primes
Cino hilliard
The C program is at the end
From T.R Nicely page
Li_2(x) = integral(2*c_2/(ln t)^2, t, 2, x), the Hardy-Littlewood integral approximation for pi_2(x).
Although this is the traditional formula, note that a slightly more accurate (for small x) approximation
is produced by the asymptotically equivalent formula Li_2*(x) = integral(2*c_2/((ln(t+6))^2), t, 2, x) .
c_2 = Hardy-Littlewood constant for twins = 0.660161815846869573927812110014555778432623...
The kth Hardy-Littlewood constant (k > 1) is defined as c_k = prod((p^(k-1))*(p-k)/(p-1)^k, p; p prime, p > k) .
The value c_2 is also referred to as the twin-primes constant; however, some authors
use 2*c_2 as the twin-primes constant.
Table showing the Largest known number of twin primes < n p12(n) from T.R Nicely page
n Pi(10^n) Pi2(10^n) 2C2Li2(10^n) error SumTP(10^2n)/4
1, 4, 2, 5, 150.00
2, 25, 8, 14, 75.00
3, 168, 35, 46, 31.43 24676
4, 1229, 205, 214, 4.39
5, 9592, 1224, 1249, 2.04
6, 78498, 8169, 8248, 0.97
7, 664579, 58980, 58754, -0.38
8, 5761455, 440312, 440368, 0.013
9, 50847534, 3424506, 3425308, 0.023
10,455052511, 27412679, 27411417, -0.0046
11,4118054813, 224376048, 224368865, -0.0032 222430448
12,37607912018, 1870585220, 1870559867, -0.0013 1856091662
13,346065536839, 15834664872, 15834598305, -0.00042 15900559214
14,3204941750802, 135780321665, 135780264894, -0.000042 136419649148
15,29844570422669, 1177209242304, 1177208491861, -0.000064 1176717615908
16,279238341033925, 10304195697298, 10304192554496, -0.000031 10301443659233
2C2Li2(n) appears to be the best estimator for the number of twin primes < n.
The link between the sum of primes < n and the number of primes < n^2 is amazing.
Indeed, Pi(10^2n) ~ Sum(Pi(n)).
This link can be extended to sums over twin primes. The following table shows the
largest known sums of twin primes < n SumTP(n) and the number of twin primes < n or Pi2(n).
Analyzing the tables and inserting estimates we will find that SumTP(n) ~ 4*Pi2(n^2)
Using this we can estimate Pi2(10^n) for n=17 to 25. This is SumTP(n)/4.
For the odd values, we will take log10 of SumPI and curvefit that to a polynomial
of appropriate degree, that is, running combinations of terms that produces the
least srtandard error and evaluate that polynomial to find the estimated log of
est(i) i = 6.5 to 12.5 in increments of. Then we compute 10^est(i) to get the
interpolation for odd values of n. Actually, 6.5 is an extrapolation.
Relationship between number of twin primes < 10^2n and sum of twin primes < 10^n or SumTP(n)
n SumTP(10^n) 4*Pi2(10^2n) Error
1 20 32 -0.60000
2 488 820 -0.68032
3 24236 32676 -0.34824
4 1726412 1761248 -0.02017
5 109114568 109650716 -0.00491
6 7424366648 7482340880 -0.00780
7 545678596592 543121286660 0.00468
8 41205774636932 41216742789192 -0.00026
9 3234489739234676
10 260643410442091112
11 21446976192435396140
12 1795640886305783918948
13 152542601906447626814216
14 13119246582832293524505360
Actual Estimate Relative Li_2(10^n) Relative
n Pi2(n) SumTP(10^2n)/4 error error
1 2 5 -1.2500
2 8 5 0.3750000 14 -0.7500
3 35 35 46 -0.3142
4 205 122 0.4048780 214 0.0439
5 1224 1322 1249 0.0204
6 8169 6959 0.1481209 8258
7 58980 52631 58754
8 440312 431603 0.0197791 440368
9 3424506 3310013 3425308
10 27412679 27278642 0.0048895 27411417
11 224376048 222430448 224368865
12 1870585220 1856091662 0.0077481 1870559867
13 15834664872 5900559214 15834598305
14 135780321665 136419649148 -0.0047085 135780264894
15 1177209242304 1176717615908 1177208491860 0.000000637477
16 10304185697298 10301443659233 0.0002661 10304192554495 -0.000000665476
n Pi2(10^n) SumT(10^n)
1.5 5 140
2.5 20 5288
3.5 83 210524
4.5 487 13240052
5.5 3133 889721792
6.5 21900 63602236856
7.5 160118 4706870463632
8.5 1222534 363550236216980
9.5 9659034 28954546423236536
10.5 78226897 2358474621067059692
11.5 646502826 195852235960038407492
12.5 5432867356 16522508800595010454128
13.5
14.5
15.5
Estimates of Pi2(n) = Pari Li_2(x)=intnum(t=2,x,2*c_2/log(t)^2)
c_2=0.6601618158468695739278121100
n SumTP(10^(n/2))/4 Li_2(10^n) Error
17 90887559054245, 90948833260990 0.000673721746
18 808622434808669, 808675901493606 0.000066116332
19 7238636605809134, 7237518062753712 -0.000154547877
20 65160852610522778, 65154265428712141 -0.000101101313
21 589618655266764923, 589629958672165100 0.000019170337
22 5361744048108849035, 5361443463012475433 -0.000056064210
23 48963058990009601873, 48962239581826462374 -0.000016735512
24 448910221576445979737, 448905221840875225106 -0.000011137619
25 4130627200148752613532, 4130654449164037217040 0.000006596779
26 38135650476611906703554, 38135381367225646263727 -0.000007056685
27
28 3279811645708073381126340, 3279812193898996337181839 0.000000167140
First 10^n lesser twin primes and some sums using c:\twinprimes\SumtwinsFirstN
a b c
n pi2x(10^n) pi2(a) pi2x(a) SumT(10^n) SumT(c)
1 107 10 4091 908 383912
2 3821 100 402341 328184 1271574140
3 79559 1000 14155367 69004076 1042100237012
4 1260989 10000 327329087 11556327260 388297227355664
5 18409199 100000 6448616177 1707243198956 64048524111768392
6 252427601 1000000 237064232862404
7 3285916169 31153163750203064
8 41375648687 3947120494191630260
9 507575862527 486665774050923191336
10 6100479510551 58727077924545848315800
11
Cino's
n pi2x(10^n) twinx(10^n) error
1 107, 100 0.0654205
2 3821, 3380 0.1154148
3 79559, 75610 0.0496361
4 1260989, 1257632 0.0026621
5 18409199, 18456351 -0.0025613
6 252427601, 252177334 0.0009914
7 3285916169, 3285912624 0.000001078
8 41375648687, 41374714817 0.000022570
9 507575862527, 507584081641 -0.000016192
10 6100479510551, 6100475249386 0.000000698
For reference pi(x)
n Pi(10^n) Rg(10^n) Relative Error
1, 4, 5, -0.25000
2, 25, 26, -0.04000
3, 168, 168, 0.00000
4, 1229, 1227, 0.16273 E-3
5, 9592, 9587, 0.52126 E-4
6, 78498, 78527, -0.36943 E-4
7, 664579, 664667, -0.13241 E-4
8, 5761455, 5761552, -0.16836 E-5
9, 50847534, 50847455, 0.15536 E-6
10, 455052511, 455050683, 0.40171 E-7
11, 4118054813, 4118052495, 0.56288 E-7
12, 37607912018, 37607910542, 0.92471 E-8
13, 346065536839, 346065531066, 0.16681 E-8
14, 3204941750802, 3204941731602, 0.59907 E-9
15, 29844570422669, 29844570495887, -0.24533 E-9
16, 279238341033925, 279238341360977, -0.11712 E-9
17, 2623557157654233, 2623557157055978, 2.28032 E-10
18, 24739954287740860, 24739954284239494, 1.41527 E-10
19, 234057667276344607, 234057667300228940, -1.02045 E-10
20, 2220819602560918840, 2220819602556027015, 2.20271 E-12
21, 21127269486018731928, 21127269485932299724, 4.09102 E-12
22, 201467286689315906290, 201467286689188773625, 6.31034 E-13
4*22 783964159847056303858, 783964159848153692021, -1.39979 E-12
/* ************************************************gcc 4.1.2 + gmp ************************************************************************
sumtwinsLessN.c uses Modification Achim Flammenkamp's algorithm in the c code below.
Achim Flammenkamp, Bielefeld, D, achim@uni-bielefeld.de
The program sums twin primes < n from start and stop up to 2^64. The routine uses gmp to add the
primes to the required precision. Comments below with Cino show some of my additions.
Cino Hilliard
Oct 2006 - Nov 2008
*/
#include <stdio.h>
#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<<20)
#define fact .2565
#define chr64(a) ulltoa(a,str,10)
//Cino save partial sums using gmp
#define use(x) do {\
anz++;\
mpz_set_str(x2,chr64(x),10);\
mpz_add(sum,sum,x2);\
} while(0);
/*
#define use(x) do {\
anz++;\
x1=x;\
mpz_add_ui(sum,sum,x1);\
} 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,x,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,tmp,x2;
printf(" Sieve, sum of twin primes <= n in ranges between 1 and 10^18\n");
printf(" Cino Hilliard\n");
printf(" Nov, 2008\n");
printf(" Input start and stop 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=3999,sq=0;
mpz_init(sum);
mpz_init(tmp);
mpz_init(x2);
printf("start => ");
gets(start1);
start = strtoull(start1,&pEnd,10);
printf("stop => ");
gets(stop1);
stop = strtoull(stop1,&pEnd,10)+3;
system("date/t&time/t");
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(3);
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(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 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==stop)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&&ii<=stop) //absolutely necessary!
{
if(ii2+2==ii) //If it is a twin, add it.
{
mpz_set_str(x2,chr64(ii2),10); //Must convert to string. no long long in gmp
mpz_add(sum,sum,x2);
anz++;
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 ",chr64(ii3));
//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 <= Twin primes <= %s} = %.18G\n",start1,stop1,(double)anz);
//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,anz);
mpz_add_ui(sum,sum,anz);
mpz_out_str(stdout,10,sum);
printf("\n");
printf("Sec = %.6f\n",t);
}//wend(1)
return anz;
}
sumprimesGmp.c
|
