/**********************************************************************/ /* */ /* 64-bit version 2009-01-07 by Stefan Spaennare */ /* ============================================= */ /* */ /* Note, this program has about the same speed as the 32-bit */ /* version on 64-bit computers. It is significantly slower */ /* on 32-bit machines. */ /* */ /**********************************************************************/ /**********************************************************************/ /* */ /* */ /* Benchmarks added by Stefan Spaennare 2010-02-24 */ /* =============================================== */ /* */ /* Contents */ /* -------- */ /* */ /* 1. Intel Core 2 Duo at 2400 MHz */ /* 2. AMD Athlon 64 at 1800 MHz */ /* 3. Intel Celeron at 1400 MHz (in old header part below) */ /* */ /* */ /* 1. Intel Core 2 Duo */ /* ------------------- */ /* */ /* The tests were run on an Intel Core 2 Duo computer at 2400 MHz */ /* with 2048 kbyte L2 cache memory and 4096 Mbyte primary memory */ /* (DDR2 memory running at 800 MHz). The operating system was */ /* Fedora 10 Linux (64-bit). The compiler was Intel icc C/C++ 11.1 */ /* (64-bit) for Linux. The load was 1.00 (i.e. no other programs */ /* running). */ /* */ /* In these benchmark tests optimum base_bytes = 500000. */ /* Compiler options were -O2 -msse3 -pc80 */ /* (-pc80 means long double precision 80-bit). */ /* */ /* Program: n: CPU-time: */ /* ================================================================== */ /* */ /* cbrun2l32.c 10^9/10^12 22.98 s */ /* cbrun2l32.c 10^9/10^14 38.36 s */ /* cbrun2l32.c 10^9/10^16 97.49 s */ /* */ /* cbrun2l64.c 10^9/10^12 22.45 s */ /* cbrun2l64.c 10^9/10^14 37.63 s */ /* cbrun2l64.c 10^9/10^16 96.47 s */ /* */ /* cbrun3l32.c 10^9/10^12 18.39 s */ /* cbrun3l32.c 10^9/10^14 28.55 s */ /* cbrun3l32.c 10^9/10^16 64.48 s */ /* */ /* cbrun3l64.c 10^9/10^12 18.33 s */ /* cbrun3l64.c 10^9/10^14 29.88 s */ /* cbrun3l64.c 10^9/10^16 65.81 s */ /* */ /* */ /* 2. AMD Athlon 64 */ /* ---------------- */ /* */ /* The tests were run on an AMD Athlon 64 2800+ computer at 1800 */ /* MHz with 512 kbyte L2 cache memory and 1024 Mbyte primary memory */ /* (DDR memory running at 166 MHz). The operating system was */ /* Fedora 10 Linux (32-bit). The compiler was Intel icc C/C++ 11.1 */ /* (32-bit) for Linux. The load was 1.00 (i.e. no other programs */ /* running). */ /* */ /* In these benchmark tests optimum base_bytes = 50000 for n=10^12 */ /* and 500000 for n=10^14 or 10^16. Compiler options were: */ /* -O2 -march=pentium3 -mtune=pentium3 -pc80 */ /* (-pc80 means long double precision 80-bit). For some reason this */ /* computer was faster with Pentium 3 optimizations than Pentium 4 */ /* optimizations. */ /* */ /* Program: n: CPU-time: */ /* ================================================================== */ /* */ /* cbrun2l32.c 10^9/10^12 29.12 s */ /* cbrun2l32.c 10^9/10^14 58.48 s */ /* cbrun2l32.c 10^9/10^16 120.18 s */ /* */ /* cbrun2l64.c 10^9/10^12 103.69 s */ /* cbrun2l64.c 10^9/10^14 171.71 s */ /* cbrun2l64.c 10^9/10^16 303.57 s */ /* */ /* cbrun3l32.c 10^9/10^12 33.10 s */ /* cbrun3l32.c 10^9/10^14 57.85 s */ /* cbrun3l32.c 10^9/10^16 73.31 s */ /* */ /* cbrun3l64.c 10^9/10^12 88.29 s */ /* cbrun3l64.c 10^9/10^14 126.79 s */ /* cbrun3l64.c 10^9/10^16 199.35 s */ /* */ /* */ /**********************************************************************/ /**********************************************************************/ /* */ /* */ /* 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-09 */ /* */ /* */ /* 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 #define LDOUBLE 1 #if LDOUBLE #define double long double #define sqrt sqrtl #define floor floorl #define fabs fabsl #define log logl #endif double timetic=(double)(CLOCKS_PER_SEC); double julian_day(year,month,day) int year,month,day; { double f,g,a,jdd; if (month>=3) { f=(double)(year); g=(double)(month); } else { f=(double)(year-1); g=(double)(month+12); } /* if */ a=2.0-floor(f/100.0)+floor(f/400.0); jdd=floor(365.25*f)+floor(30.6001*(g+1.0))+ (double)(day)+a+1720994.5; return(jdd); } /* julian_day */ double get_float_days(ttime) time_t *ttime; { int year,month,day,hour,min,sec; double fdays,fhours; struct tm *tt; tt=localtime(ttime); year=tt->tm_year+1900; month=tt->tm_mon+1; day=tt->tm_mday; hour=tt->tm_hour; min=tt->tm_min; sec=tt->tm_sec; fhours=(double)(hour)+(double)(min)/60.0+(double)(sec)/3600.0; fdays=julian_day(year,month,day)+fhours/24.0; return(fdays); } /* get_float_days */ void day_and_time(ttime) time_t *ttime; { int year,month,day,hour,min,sec; struct tm *tt; tt=localtime(ttime); year=tt->tm_year+1900; month=tt->tm_mon+1; day=tt->tm_mday; hour=tt->tm_hour; min=tt->tm_min; sec=tt->tm_sec; printf("%04d-%02d-%02d %02d:%02d:%02d",year,month,day,hour,min,sec); } /* day_and_time */ void fday_and_time(outfile,ttime) FILE *outfile; time_t *ttime; { int year,month,day,hour,min,sec; struct tm *tt; tt=localtime(ttime); year=tt->tm_year+1900; month=tt->tm_mon+1; day=tt->tm_mday; hour=tt->tm_hour; min=tt->tm_min; sec=tt->tm_sec; fprintf(outfile,"%04d-%02d-%02d %02d:%02d:%02d",year,month,day,hour,min,sec); } /* fday_and_time */ int main(argc,argv) int argc; char *argv[]; { FILE *outfile; unsigned long long i,j,ii,kk,nsqrt,base,base_bytes; unsigned long long off,j0,flag,ltemp,base4,base41; unsigned long long nsqrt2,base1,nmod,nmod1,hh,ll,bit; double nd,nnd,kd,based,ndivd,nmodd; double psumd,nsumd; double twinsd,poldd,nlimd; double brun,b2,c2,temp,ttime; double elapseddays,lt; unsigned long long *pflag; unsigned long long *oldj; clock_t time1,timediff1; time_t start_time,stop_time; #if LDOUBLE double strtold(); #else double atof(); #endif if (argc != 4) { printf("Usage: %s out_file n(>=prim-max)r base_bytesi\n",argv[0]); exit(0); } /* if */ outfile=fopen(argv[1],"w"); #if LDOUBLE nd=strtold(argv[2],NULL); #else nd=atof(argv[2]); #endif base_bytes=(unsigned long long)(atoi(argv[3])); nd=floor(nd+0.5); base4=base_bytes/8ull; base=64ull*base4; based=(double)(base); if ((nd+0.5) < based) { based=nd; base=(unsigned long long)(based+0.5); } /* if */ ndivd=floor((nd+0.5)/based); nmodd=nd-ndivd*based; nnd=ndivd*based; pflag=(unsigned long long *)calloc(base4+3ull,sizeof(unsigned long long)); nsqrt=(unsigned long long)(sqrt(nd))+1ull; nsqrt2=nsqrt/2ull; oldj=(unsigned long long *)calloc(nsqrt2+3ull,sizeof(unsigned long long)); time1=clock(); time(&start_time); printf("\n"); printf("Bruns constant and the Twin prime constant\n"); printf("==========================================\n"); printf("\n"); #if LDOUBLE printf("Long double precision (80 bit).\n"); printf("\n"); printf("n = %19.0Lf = %9.3LE; base (bytes) = %10Lu; n base eff. = %10Lu\n",nd,nd,base_bytes,base); #else printf("Double precision (64 bit).\n"); printf("\n"); printf("n = %19.0lf = %9.3lE; base (bytes) = %10Lu; n base eff. = %10Lu\n",nd,nd,base_bytes,base); #endif printf("\n"); printf("Start time: "); day_and_time(&start_time); printf("\n"); printf("\n"); printf(" s CPU lm m Primes <= m Prime twins <= m Bruns constant B2 Twin prime const.\n"); printf("-----------------------------------------------------------------------------------------------------------------------------------\n"); fprintf(outfile,"\n"); fprintf(outfile,"Bruns constant and the Twin prime constant\n"); fprintf(outfile,"==========================================\n"); fprintf(outfile,"\n"); #if LDOUBLE fprintf(outfile,"Long double precision (80 bit).\n"); fprintf(outfile,"\n"); fprintf(outfile,"n = %19.0Lf = %9.3LE; base (bytes) = %10Lu; n base eff. = %10Lu\n",nd,nd,base_bytes,base); #else fprintf(outfile,"Double precision (64 bit).\n"); fprintf(outfile,"\n"); fprintf(outfile,"n = %19.0lf = %9.3lE; base (bytes) = %10Lu; n base eff. = %10Lu\n",nd,nd,base_bytes,base); #endif fprintf(outfile,"\n"); fprintf(outfile,"Start time: "); fday_and_time(outfile,&start_time); fprintf(outfile,"\n"); fprintf(outfile,"\n"); fprintf(outfile," s CPU lm m Primes <= m Prime twins <= m Bruns constant B2 Twin prime const.\n"); fprintf(outfile,"-----------------------------------------------------------------------------------------------------------------------------------\n"); fflush(NULL); nsumd=2.0; psumd=1.0; twinsd=0.0; off=4ull; base1=base+1ull; base41=base4+1ull; nmod=(unsigned long long)(nmodd+0.5); nmod1=nmod+1ull; kd=0.0; flag=0ull; nlimd=10.0; poldd=0.0; b2=0.0; c2=1.0; brun=0.0; kk=1ull; lt=log(10.0); ii=1ull; for (i=3ull; i<=nsqrt; i=i+2ull) { oldj[ii]=2ull*i+1ull; ii++; } /* for i */ do { if (fabs(kd-nnd) < 0.5) { base=nmod; base1=nmod1; based=nmodd; flag=1ull; } /* if */ for (j=1ull; j<=base41; j++) { pflag[j]=0ull; } /* for j */ ii=1ull; for (i=3ull; i<=nsqrt; i=i+2ull) { j0=oldj[ii]; for (j=j0; j<=base1; j=j+i) { hh=j / 64ull + 1ull; ll=j % 64ull; bit= 1ull << ll; pflag[hh]=pflag[hh] | bit; } /* for j */ oldj[ii]=j-base; ii++; } /* for i */ for (j=off; j<=base1; j++) { nsumd=nsumd+1.0; if ((j % 2ull) == 0ull) { hh=j / 64ull + 1ull; ll=j % 64ull; bit= 1ull << ll; ltemp=pflag[hh] & bit; if (ltemp == 0ull) { psumd=psumd+1.0; temp=nsumd-1.0; c2=c2*(1.0-1.0/(temp*temp)); if (fabs(nsumd-poldd) < 2.5) { twinsd=twinsd+1.0; b2=b2+1.0/nsumd+1.0/poldd; } /* if */ poldd=nsumd; } /* if */ } /* if */ if (fabs(nlimd-nsumd) < 0.5) { brun=b2+4.0*c2/log(nsumd); timediff1=clock()-time1; ttime=(double)(timediff1)/timetic; #if LDOUBLE printf("%10.2Lf %2Lu %17.0Lf %17.0Lf %17.0Lf %18.16Lf %18.16Lf %18.16Lf\n",ttime,kk,nsumd,psumd,twinsd,brun,b2,c2); fprintf(outfile,"%10.2Lf %2Lu %17.0Lf %17.0Lf %17.0Lf %18.16Lf %18.16Lf %18.16Lf\n",ttime,kk,nsumd,psumd,twinsd,brun,b2,c2); #else printf("%10.2lf %2Lu %17.0lf %17.0lf %17.0lf %18.16lf %18.16lf %18.16lf\n",ttime,kk,nsumd,psumd,twinsd,brun,b2,c2); fprintf(outfile,"%10.2lf %2Lu %17.0lf %17.0lf %17.0lf %18.16lf %18.16lf %18.16lf\n",ttime,kk,nsumd,psumd,twinsd,brun,b2,c2); #endif fflush(NULL); nlimd=10.0*nlimd; kk++; } /* if */ } /* for j */ kd=kd+based; off=2ull; } while (flag == 0ull); brun=b2+4.0*c2/log(nsumd); timediff1=clock()-time1; ttime=(double)(timediff1)/timetic; #if LDOUBLE printf("%10.2Lf ln %17.0Lf %17.0Lf %17.0Lf %18.16Lf %18.16Lf %18.16Lf\n",ttime,nsumd,psumd,twinsd,brun,b2,c2); fprintf(outfile,"%10.2Lf ln %17.0Lf %17.0Lf %17.0Lf %18.16Lf %18.16Lf %18.16Lf\n",ttime,nsumd,psumd,twinsd,brun,b2,c2); #else printf("%10.2lf ln %17.0lf %17.0lf %17.0lf %18.16lf %18.16lf %18.16lf\n",ttime,nsumd,psumd,twinsd,brun,b2,c2); fprintf(outfile,"%10.2lf ln %17.0lf %17.0lf %17.0lf %18.16lf %18.16lf %18.16lf\n",ttime,nsumd,psumd,twinsd,brun,b2,c2); #endif printf("-----------------------------------------------------------------------------------------------------------------------------------\n"); fprintf(outfile,"-----------------------------------------------------------------------------------------------------------------------------------\n"); time(&stop_time); elapseddays=get_float_days(&stop_time)-get_float_days(&start_time); printf("\n"); printf("Stop time: "); day_and_time(&stop_time); printf("\n"); printf("\n"); #if LDOUBLE printf("Elapsed clock time: %18.8Lf days = %18.6Lf hours = %18.2Lf seconds\n", elapseddays,24.0*elapseddays,24.0*3600.0*elapseddays); #else printf("Elapsed clock time: %18.8lf days = %18.6lf hours = %18.2lf seconds\n", elapseddays,24.0*elapseddays,24.0*3600.0*elapseddays); #endif printf("\n"); fprintf(outfile,"\n"); fprintf(outfile,"Stop time: "); fday_and_time(outfile,&stop_time); fprintf(outfile,"\n"); fprintf(outfile,"\n"); #if LDOUBLE fprintf(outfile,"Elapsed clock time: %18.8Lf days = %18.6Lf hours = %18.2Lf seconds\n", elapseddays,24.0*elapseddays,24.0*3600.0*elapseddays); #else fprintf(outfile,"Elapsed clock time: %18.8lf days = %18.6lf hours = %18.2lf seconds\n", elapseddays,24.0*elapseddays,24.0*3600.0*elapseddays); #endif fprintf(outfile,"\n"); fflush(NULL); fclose(outfile); } /* End */