/**********************************************************************/ /* */ /* */ /* Some fast programs for generating prime numbers, calculating */ /* ============================================================ */ /* the number of twin primes and Bruns constant */ /* ============================================ */ /* */ /* */ /* Author: */ /* ------- */ /* */ /* Stefan Spaennare, Lund, Sweden */ /* */ /* E-mail: stefan@spaennare.se */ /* */ /* Home page: http://www.spaennare.se/index.html */ /* Program page: http://www.spaennare.se/ssprog.html */ /* */ /* First version: June 1998 */ /* Latest update: 2007-08-13 */ /* */ /* */ /* Purpose: */ /* -------- */ /* */ /* This header is common for the following programs: */ /* */ /* esieve.c: Generating, counting and printing prime numbers. */ /* */ /* ptwins.c: Prime number and twin primes statistics. */ /* */ /* cbrun.c, cbrun2.c and cbrun3.c: Calculating prime number and */ /* twins primes statistics and Bruns constant (double precision */ /* versions). */ /* */ /* cbrunl.c, cbrun2l.c and cbrun3l.c: Calculating prime number and */ /* twin primes statistics and Bruns constant (long double precision */ /* versions). Recommended versions for large n (>=10^12)! */ /* */ /* */ /* Input variables: */ /* ---------------- */ /* */ /* Input variables ending with "i" are integers and "r" are real */ /* i.e. double or long double precision. */ /* */ /* out_file: Output file for the results. */ /* */ /* n(>=prim-max)r: The program calculates and use all primes <= n. */ /* */ /* base_bytesi: Tells how large part of n that fits into the CPUs L2 */ /* cache. A good starting value is about 80 % of the size of the L2 */ /* cache rounded to nearest 100000 bytes. */ /* */ /* The programs cbrun2.c, cbrun3.c, cbrun2l.c and cbrun3l.c use a */ /* bit-wise handling of the cache memory making it 8 and 16 times */ /* larger effectively. This makes these programs significantly */ /* faster than cbrun.c and cbrunl.c for large n (i.e. n >= 10^12). */ /* The programs esieve.c and ptwins.c do not have this feature and */ /* are slower. */ /* */ /* */ /* Algorithm, numbers and precision: */ /* --------------------------------- */ /* */ /* The programs use a fast implementation of the segmented */ /* Eratosthenes sieve for calculating prime numbers. */ /* */ /* Theoretical time complexity: O(n*loglog(n)) */ /* */ /* Only the programs cbrunl.c, cbrun2l.c and cbrun3l.c have the */ /* long double feature. */ /* */ /* At maximum n <= 2^53-1 (about 9.0*10^15) if double is used or */ /* n <= 2^64-1 (about 1.8*10^19) if long double is used. However */ /* calculations for n >= 1.0*10^12 will take long time also on */ /* very fast computers. */ /* */ /* Large numbers are stored as double (64 bit) or long double */ /* (80 bit) numbers depending on the flag LDOUBLE (0 or 1). */ /* A double precision floating point number is supposed to have */ /* 53 bit in the mantissa and 64 bit for long double. */ /* */ /* Important note! It is highly recommended to run these program with */ /* long double precision. Otherwise there will be precision loss in */ /* Bruns constant and the Twin prime constant for large n (>= 10^12). */ /* */ /* */ /* Memory requirements: */ /* -------------------- */ /* */ /* Memory requirements = 2*sqrt(n)+base bytes */ /* */ /* */ /* Compiling: */ /* ---------- */ /* */ /* >gcc -O3 -o program program.c -lm */ /* */ /* */ /* References: */ /* ----------- */ /* */ /* The web page "Mathematical Constants and Computation" by */ /* Xavier Gourdon and Pascal Sebah: */ /* */ /* http://numbers.computation.free.fr/Constants/constants.html */ /* */ /* */ /* Benchmarks: */ /* ----------- */ /* */ /* The tests were run on an Intel Celeron (Tualatin core) computer */ /* at 1400 MHz with 256 kbyte L2 cache memory (Advanced Transfer) */ /* and 576 Mbyte primary memory (PC133 memory running at 100 MHz). */ /* The operating system was Fedora Core 3 Linux. Compilers tested */ /* were GNU GCC 3.4.3 (gcc) and Intel C/C++ 8.1 (icc) for Linux. */ /* The load was 1.00 (i.e. no other programs running). */ /* */ /* In these benchmark tests base_bytes = 200000 (a quite optimum */ /* value for this computer). */ /* */ /* Double precision results: */ /* ------------------------- */ /* */ /* Program: n: CPU-time: Compiler options: */ /* ================================================================== */ /* */ /* esieve.c 10^9 38.87 s gcc -O3 */ /* ptwins.c 10^9 46.87 s gcc -O3 */ /* */ /* cbrun.c 10^9/10^12 140.22 s gcc -O3 */ /* cbrun2.c 10^9/10^12 88.90 s gcc -O3 */ /* cbrun3.c 10^9/10^12 79.56 s gcc -O3 */ /* */ /* */ /* Long double precision results: */ /* ------------------------------ */ /* */ /* Program: n: CPU-time: Compiler options: */ /* ================================================================== */ /* */ /* cbrunl.c 10^9/10^12 170.39 s gcc -O3 */ /* cbrun2l.c 10^9/10^12 106.78 s gcc -O3 */ /* cbrun3l.c 10^9/10^12 76.65 s gcc -O3 */ /* */ /* cbrunl.c 10^9/10^14 800.44 s gcc -O3 */ /* cbrun2l.c 10^9/10^14 189.53 s gcc -O3 */ /* cbrun3l.c 10^9/10^14 121.47 s gcc -O3 */ /* */ /* cbrunl.c 10^9/10^16 --- gcc -O3 */ /* cbrun2l.c 10^9/10^16 910.69 s gcc -O3 */ /* cbrun3l.c 10^9/10^16 473.04 s gcc -O3 */ /* *************/ /* */ /* cbrunl.c 10^9/10^12 145.51 s icc -O2 -tpp6 -xK -ip -static -long_double */ /* cbrun2l.c 10^9/10^12 84.60 s icc -O2 -tpp6 -xK -ip -static -long_double */ /* cbrun3l.c 10^9/10^12 67.21 s icc -O2 -tpp6 -xK -ip -static -long_double */ /* */ /* cbrunl.c 10^9/10^14 805.14 s icc -O2 -tpp6 -xK -ip -static -long_double */ /* cbrun2l.c 10^9/10^14 195.65 s icc -O2 -tpp6 -xK -ip -static -long_double */ /* cbrun3l.c 10^9/10^14 132.33 s icc -O2 -tpp6 -xK -ip -static -long_double */ /* */ /* cbrunl.c 10^9/10^16 --- icc -O2 -tpp6 -xK -ip -static -long_double */ /* cbrun2l.c 10^9/10^16 876.91 s icc -O2 -tpp6 -xK -ip -static -long_double */ /* cbrun3l.c 10^9/10^16 464.33 s icc -O2 -tpp6 -xK -ip -static -long_double */ /* */ /**********************************************************************************/ /********************************************************************************/ /* */ /* Notice */ /* ====== */ /* */ /* I make no warranties that this program is (1) free of errors, (2) consistent */ /* with any standard merchantability, or (3) meeting the requirements of a */ /* particular application. This software shall not, partly or as a whole, */ /* participate in a process, whose outcome can result in injury to a person or */ /* loss of property. It is solely designed for analytical work. Permission to */ /* use, copy and distribute is hereby granted without fee, providing that the */ /* header above including this notice appears in all copies. */ /* */ /* Stefan Spaennare */ /* */ /********************************************************************************/ #include #include #include #include #include double timetic=(double)(CLOCKS_PER_SEC); int main(argc,argv) int argc; char *argv[]; { FILE *outfile; unsigned long i,j,ii,nsqrt,base; unsigned long off,j0,flag,flag2; unsigned long nsqrt2,base1,nmod,nmod1; unsigned long c1,c2; double nd,nnd,kd,based,ndivd,nmodd; double psumd,nsumd; double twinsd,poldd,stepd,nlimd; unsigned char *pflag; unsigned long *oldj; clock_t time1,timediff1; if (argc != 4) { printf("Usage: %s out_file n(>=prim-max)r base_bytesi\n",argv[0]); exit(0); } /* if */ nd=atof(argv[2]); base=(unsigned long)(atoi(argv[3])); nd=floor(nd+0.5); printf("n = %19.0f\n",nd); outfile=fopen(argv[1],"w"); based=(double)(base); if ((nd+0.5) < based) { based=nd; base=(unsigned long)(based+0.5); } /* if */ ndivd=floor((nd+0.5)/based); nmodd=nd-ndivd*based; nnd=ndivd*based; pflag=(unsigned char *)calloc(base+3,sizeof(unsigned char)); nsqrt=(unsigned long)(sqrt(nd))+1; nsqrt2=nsqrt/2; oldj=(unsigned long *)calloc(nsqrt2+3,sizeof(unsigned long)); time1=clock(); nsumd=2.0; psumd=1.0; twinsd=0.0; off=4; base1=base+1; nmod=(unsigned long)(nmodd+0.5); nmod1=nmod+1; kd=0.0; flag=0; nlimd=3.0; c1=3; c2=100; stepd=1.0; poldd=0.0; fprintf(outfile," n primes <= n prime-twins <= n\n"); fprintf(outfile,"-----------------------------------------------------\n"); fprintf(outfile,"%17.0f %17.0f %17.0f\n",nsumd,psumd,twinsd); fflush(NULL); ii=1; for (i=3; i<=nsqrt; i=i+2) { oldj[ii]=2*i+1; ii++; } /* for i */ do { if (fabs(kd-nnd) < 0.5) { base=nmod; base1=nmod1; based=nmodd; flag=1; } /* if */ for (j=off; j<=base1; j=j+2) { pflag[j]=0; } /* for j */ ii=1; for (i=3; i<=nsqrt; i=i+2) { j0=oldj[ii]; for (j=j0; j<=base1; j=j+i) { pflag[j]=1; } /* for j */ oldj[ii]=j-base; ii++; } /* for i */ for (j=off; j<=base1; j++) { nsumd=nsumd+1.0; if ((j % 2) == 0) { if (pflag[j]==0) { psumd=psumd+1.0; if (fabs(nsumd-poldd) < 2.5) { twinsd=twinsd+1.0; } /* if */ poldd=nsumd; } /* if */ } /* if */ flag2=0; if (fabs(nlimd-nsumd) < 0.5) { fprintf(outfile,"%17.0f %17.0f %17.0f\n",nsumd,psumd,twinsd); fflush(NULL); flag2=1; } /* if */ if (c1==c2) { stepd=10.0*stepd; nlimd=10.0*stepd; c1=1; c2=91; } /* if */ if (flag2==1) { nlimd=nlimd+stepd; c1++; } /* if */ } /* for j */ kd=kd+based; off=2; } while (flag==0); fclose(outfile); timediff1=clock()-time1; printf("primes = %14.0f\n",psumd); printf("twins = %14.0f\n",twinsd); printf("CPU time = %12.2f s\n",(double)(timediff1)/timetic); fflush(NULL); } /* End */