#include #include /*************************************************************************** * * * Program to verify the Riemann hypothesis up to Graham point END. The * * main routine uses the Riemann-Siegel formula for fast evaluation of * * the zeta function. * * * * This program counts the number of zeros of the Riemann zeta function * * between Graham points BEGIN and END. The program establishes the * * existence of a zero with complete certainty, and thus produces a lower * * bound on the number of zeros in the search range. There is the * * possibility that a zero will be missed in the case where |Z(t)| is * * less than the error tolerance, and hence the sign of Z(t) cannot be * * determined with certainty. In these situations, a warning is written * * to the RESULTFILE indicating that the questionable region should be * * examined with more accuracy. * * * * Output to the RESULTFILE consists of the count of roots found, the * * number of Graham blocks found, the total number of roots within these * * Graham blocks, and the longest and shortest Graham blocks found. * * * * on RedHat Linux compile with * * * * gcc -lm -o rootcount rootcount7.c * * * * and run with * * * * ./rootcount * * * ***************************************************************************/ #define BEGIN -1 /* initial graham point */ #define END 12193873 /* final graham point */ #define RESULTFILE "results.txt" /* output file for results */ /**************************************************************************/ int main() { FILE *results_file; /* output file for results */ long double Z(long double,int); /* Riemann-Siegel Z(t) function */ long double graham(int); /* graham point location function */ int grahamblock(int,int,FILE *); /* graham block root counting function */ int even(int); /* -1,+1 parity function */ int gpt_final; /* index of final graham point */ int gpt_indx; /* index of current graham point */ int tot_roots; /* root counter */ int gblock_start; /* current graham block lower index */ int gblock_end; /* current graham block upper index */ int gblock_roots; /* roots in current graham block */ int tot_gblocks; /* total graham blocks found */ int tot_gblock_roots; /* total graham block roots found */ int l_gblock_start; /* longest graham block left index */ int l_gblock_end; /* longest graham block right index */ gpt_indx = BEGIN; /* initializations */ gpt_final = END; /* */ tot_roots = 0; /* */ gblock_roots = 0; /* */ tot_gblocks = 0; /* */ tot_gblock_roots = 0; /* */ l_gblock_start = -1; /* */ l_gblock_end = -1; /* */ results_file = fopen(RESULTFILE, "w"); /* open output file */ while (gpt_indx < gpt_final) /* while not at the */ { /* last Graham point... */ ++gpt_indx; /* move to the next one */ if (even(gpt_indx)*Z(graham(gpt_indx),1) /* if Graham's Law */ > 0.053L/powl(graham(gpt_indx),1.25L)) /* satified, increment */ ++tot_roots; /* total number roots, */ else /* otherwise, try more */ if (even(gpt_indx)*Z(graham(gpt_indx),4) /* precise calculation */ > 0.053L/powl(graham(gpt_indx),1.25L)) /* of Z(t)... */ ++tot_roots; /* */ else /* still not satisfied? */ { /* then we have a */ /* Graham block... */ gblock_start = gpt_indx - 1; /* identify left index */ /* of Graham block */ while (even(gpt_indx) /* test succeeding */ *Z(graham(gpt_indx),1) /* Graham points until */ <= /* Graham's Law once */ 0.053L/powl(graham(gpt_indx), /* again satisfied */ 1.25L)) /* */ ++gpt_indx; /* */ gblock_end = gpt_indx; /* identify right index */ /* of Graham block */ gblock_roots = /* determine number of */ grahamblock(gblock_start,gblock_end, /* roots inside */ results_file); /* the Graham block */ ++tot_gblocks; /* increment root and */ tot_gblock_roots += gblock_roots; /* Graham block */ tot_roots += gblock_roots; /* counters */ if (gblock_end-gblock_start /* determine if the new */ > l_gblock_end-l_gblock_start) /* Graham block is */ { /* longest so far */ l_gblock_start = gblock_start; /* */ l_gblock_end = gblock_end; /* */ } /* */ } if (gpt_indx%1000 == 0) /* write progress */ printf("%d roots located to Graham point" /* message to screen */ " %d\n", /* */ tot_roots,gpt_indx); /* */ } fprintf(results_file,"Number of roots between" /* write results to */ " graham point N = %d and" /* output file: */ " graham point N = %d is %d\n", /* */ BEGIN,gpt_indx,tot_roots); /* ...total roots, */ fprintf(results_file,"Number of Graham" /* */ " blocks found is %d\n", /* ...total Graham */ tot_gblocks); /* blocks, */ fprintf(results_file,"Total number of roots" /* */ " in these %d Graham blocks is %d\n", /* ...total roots in */ tot_gblocks,tot_gblock_roots); /* Graham blocks, */ fprintf(results_file,"Longest Graham block" /* */ " found is [%d, %d]\n", /* ...longest Graham */ l_gblock_start,l_gblock_end); /* block */ fclose(results_file); /* close output file */ return(0); } int even(int n) /*************************************************************************** * * * function which returns -1 if argument is odd and +1 if argument is even * * * ***************************************************************************/ { if (n%2 == 0) return(1); else return(-1); } int grahamblock(int n, int m, FILE *outfile) /*************************************************************************** * * * Function which searches for the expected number of roots between Graham * * points n and m. Warnings of possible missed roots along with the * * Graham block are written to outfile for further analysis * * * ***************************************************************************/ { int num; /* root counter */ int p; /* partitions between graham points */ int j; /* index of current graham point */ int k; /* index of current partition point */ long double Z(long double,int); /* Riemann-Siegel Z function */ long double graham(int); /* Graham point evaluation function */ long double gp1; /* Graham interval partition point */ long double gp2; /* Graham interval partition point */ p = 1; num = 0; while (num < (m - n) && p < 1024) /* while no. roots */ { /* less than expected...*/ num = 0; /* refine inter-Graham- */ p = 2*p; /* point partition, go */ j = n; /* to beginning of */ gp2 = graham(j); /* Graham block... */ while (j <= (m-1)) /* for each Graham pt. */ { /* except the last... */ k = 0; /* step through */ while (k <= (p-1)) /* successive partition */ { /* points... */ gp1 = gp2; /* determine */ /* coordinates of */ /* adjacent partition */ /* points... */ /* */ /* */ /* */ gp2 = graham(j) /* */ * (1.0L-((long double) (k+1)) /* */ / ((long double) p)) /* */ + /* */ graham(j+1) /* */ * ((long double) (k+1)) /* */ / ((long double) p); /* */ if (Z(gp1,1) /* roughly evaluate */ > 0.053L/powl(gp1,1.25L) /* Z(t) at adjacent */ && /* partition points... */ Z(gp2,1) /* */ < -0.053L/powl(gp2,1.25L)) /* if the signs differ */ ++num; /* increment the root */ else if (Z(gp1,1) /* count... */ < -0.053L/powl(gp1,1.25L) /* */ && /* */ Z(gp2,1) /* */ > 0.053L/powl(gp2,1.25L)) /* */ ++num; /* */ ++k; /* next partition point */ } ++j; /* next Graham point */ } } p = 1; /* reset partition size */ while (num < (m - n)) /* if no. roots still */ { /* less than expected, */ /* start over using more*/ /* accurate Z(t)... */ num = 0; /* refine inter-Graham- */ p = 2*p; /* point partition, go */ j = n; /* to beginning of */ gp2 = graham(j); /* Graham block... */ while (j <= (m-1)) /* for each Graham pt. */ { /* except the last... */ k = 0; /* step through */ while (k <= (p-1)) /* successive */ { /* partition points... */ gp1 = gp2; /* determine */ /* coordinates of */ /* adjacent partition */ /* points... */ /* */ /* */ /* */ gp2 = graham(j) /* */ * (1.0L-((long double) (k+1)) /* */ / ((long double) p)) /* */ + /* */ graham(j+1) /* */ * ((long double) (k+1)) /* */ /((long double) p); /* */ if (Z(gp1,4) /* accurately evaluate */ > 0.053L/powl(gp1,1.25L) /* Z(t) at adjacent */ && /* partition points... */ Z(gp2,4) /* */ < -0.053L/powl(gp2,1.25L)) /* if the signs differ */ ++num; /* increment the root */ else if (Z(gp1,4) /* count... */ < -0.053L/powl(gp1,1.25L) /* */ && /* */ Z(gp2,4) /* */ > 0.053L/powl(gp2,1.25L)) /* */ ++num; /* */ ++k; /* next partition point */ } ++j; /* next Graham point */ } if (p == 65536) break; /* break if partition */ /* depth getting too */ } /* deep */ if (num < (m - n)) /* if unable to find */ { /* expected number of */ fprintf(outfile, /* roots, write warning */ "Partition depth reached in" /* to output file */ " Graham block [%d,%d],\n" /* */ " possible violation of" /* */ " Rossers rule\n",n,m); /* */ } return(num); } long double Z(long double t, int n) /*************************************************************************** * * * The Z(t) function from the Riemann-Siegel formula. This functions takes * * an additional integer argument which is the number of terms to use in * * the remainder. This integer argument can vary from 0 to 4 which * * corresponds to the first through fifth remainder terms. * * * ***************************************************************************/ { long double ZZ; /* the result */ long double tt; /* theta(t) */ long double p; /* fractional part of sqrt(t/(2.0*pi)) */ long double theta(long double); /* theta function */ long double C(int,long double); /* coefficient of (2*pi/t)^(k*0.5) */ long double R; /* remainder term */ int even(int); /* -1,+1 parity function */ int j; /* summation index for Z(t) function */ int k; /* summation index for remainder term */ int N; /* integer part of sqrt(t/(2.0*pi)) */ const long double pi = 3.1415926535897932385L; /* initializations... */ ZZ = 0.0L; /* */ R = 0.0L; /* */ j = 1; /* */ k = 0; /* */ N = sqrtl(t/(2.0L * pi)); /* */ p = sqrtl(t/(2.0L * pi)) - N; /* */ tt = theta(t); /* */ while (j <= N) /* add up terms of */ { /* main series... */ ZZ = ZZ + 1.0L/sqrtl((long double) j ) /* */ * cosl(fmodl(tt /* */ -t*logl((long double) j), /* */ 2.0L*pi)); /* */ ++j; /* */ } /* */ ZZ = 2.0L * ZZ; /* */ while (k <= n) /* add up terms of */ { /* remainder... */ R = R + C(k,2.0L*p-1.0L) /* */ * powl(2.0L*pi/t, /* */ ((long double) k)*0.5L); /* */ ++k; /* */ } /* */ R = even(N-1) * powl(2.0L * pi / t,0.25L) * R; /* */ return(ZZ + R); } long double theta(long double t) /*************************************************************************** * * * Approximation to theta(t)=Im{log[Pi(it/2-3/4)]}-t/2*log(pi) * * * ***************************************************************************/ { const long double pi = 3.1415926535897932385L; return(t/2.0L*logl(t/2.0L/pi) - t/2.0L - pi/8.0L + 1.0L/48.0L/t + 7.0L/5760.0L/t/t/t); } long double graham(int n) /*************************************************************************** * * * Function which locates n^th Graham point using Newton's Method. * * * ***************************************************************************/ { long double t_n; /* n^th approximation */ long double t_n1; /* (n+1)^st approx. */ const long double pi = 3.1415926535897932385L; t_n=0.0L; t_n1=0.5L*n+20.0L; while(fabsl(t_n-t_n1) > 0.00000001L) { t_n=t_n1; t_n1 = t_n-(t_n*log(t_n/pi/2.0L)/2.0L-t_n/2.0L- pi/8.0L+1.0L/t_n/48.0L+7.0L/5760.0L/(t_n*t_n*t_n)- ((long double) n)*pi) /(log(t_n/pi/2.0L)/2.0L-1.0L/(t_n*t_n)/48.0L- 7.0L/1920.0L/(t_n*t_n*t_n*t_n)); } return(t_n1); } long double C(int n, long double z) /*************************************************************************** * * * Coefficients of remainder terms; n can range from 0 to 4. * * * ***************************************************************************/ { if (n==0) return(.38268343236508977173L * powl(z, 0.0L) +.43724046807752044936L * powl(z, 2.0L) +.13237657548034352332L * powl(z, 4.0L) -.01360502604767418865L * powl(z, 6.0L) -.01356762197010358089L * powl(z, 8.0L) -.00162372532314446528L * powl(z,10.0L) +.00029705353733379691L * powl(z,12.0L) +.00007943300879521470L * powl(z,14.0L) +.00000046556124614505L * powl(z,16.0L) -.00000143272516309551L * powl(z,18.0L) -.00000010354847112313L * powl(z,20.0L) +.00000001235792708386L * powl(z,22.0L) +.00000000178810838580L * powl(z,24.0L) -.00000000003391414390L * powl(z,26.0L) -.00000000001632663390L * powl(z,28.0L) -.00000000000037851093L * powl(z,30.0L) +.00000000000009327423L * powl(z,32.0L) +.00000000000000522184L * powl(z,34.0L) -.00000000000000033507L * powl(z,36.0L) -.00000000000000003412L * powl(z,38.0L) +.00000000000000000058L * powl(z,40.0L) +.00000000000000000015L * powl(z,42.0L)); else if (n==1) return(-.02682510262837534703L * powl(z, 1.0L) +.01378477342635185305L * powl(z, 3.0L) +.03849125048223508223L * powl(z, 5.0L) +.00987106629906207647L * powl(z, 7.0L) -.00331075976085840433L * powl(z, 9.0L) -.00146478085779541508L * powl(z,11.0L) -.00001320794062487696L * powl(z,13.0L) +.00005922748701847141L * powl(z,15.0L) +.00000598024258537345L * powl(z,17.0L) -.00000096413224561698L * powl(z,19.0L) -.00000018334733722714L * powl(z,21.0L) +.00000000446708756272L * powl(z,23.0L) +.00000000270963508218L * powl(z,25.0L) +.00000000007785288654L * powl(z,27.0L) -.00000000002343762601L * powl(z,29.0L) -.00000000000158301728L * powl(z,31.0L) +.00000000000012119942L * powl(z,33.0L) +.00000000000001458378L * powl(z,35.0L) -.00000000000000028786L * powl(z,37.0L) -.00000000000000008663L * powl(z,39.0L) -.00000000000000000084L * powl(z,41.0L) +.00000000000000000036L * powl(z,43.0L) +.00000000000000000001L * powl(z,45.0L)); else if (n==2) return(+.00518854283029316849L * powl(z, 0.0L) +.00030946583880634746L * powl(z, 2.0L) -.01133594107822937338L * powl(z, 4.0L) +.00223304574195814477L * powl(z, 6.0L) +.00519663740886233021L * powl(z, 8.0L) +.00034399144076208337L * powl(z,10.0L) -.00059106484274705828L * powl(z,12.0L) -.00010229972547935857L * powl(z,14.0L) +.00002088839221699276L * powl(z,16.0L) +.00000592766549309654L * powl(z,18.0L) -.00000016423838362436L * powl(z,20.0L) -.00000015161199700941L * powl(z,22.0L) -.00000000590780369821L * powl(z,24.0L) +.00000000209115148595L * powl(z,26.0L) +.00000000017815649583L * powl(z,28.0L) -.00000000001616407246L * powl(z,30.0L) -.00000000000238069625L * powl(z,32.0L) +.00000000000005398265L * powl(z,34.0L) +.00000000000001975014L * powl(z,36.0L) +.00000000000000023333L * powl(z,38.0L) -.00000000000000011188L * powl(z,40.0L) -.00000000000000000416L * powl(z,42.0L) +.00000000000000000044L * powl(z,44.0L) +.00000000000000000003L * powl(z,46.0L)); else if (n==3) return(-.00133971609071945690L * powl(z, 1.0L) +.00374421513637939370L * powl(z, 3.0L) -.00133031789193214681L * powl(z, 5.0L) -.00226546607654717871L * powl(z, 7.0L) +.00095484999985067304L * powl(z, 9.0L) +.00060100384589636039L * powl(z,11.0L) -.00010128858286776622L * powl(z,13.0L) -.00006865733449299826L * powl(z,15.0L) +.00000059853667915386L * powl(z,17.0L) +.00000333165985123995L * powl(z,19.0L) +.00000021919289102435L * powl(z,21.0L) -.00000007890884245681L * powl(z,23.0L) -.00000000941468508130L * powl(z,25.0L) +.00000000095701162109L * powl(z,27.0L) +.00000000018763137453L * powl(z,29.0L) -.00000000000443783768L * powl(z,31.0L) -.00000000000224267385L * powl(z,33.0L) -.00000000000003627687L * powl(z,35.0L) +.00000000000001763981L * powl(z,37.0L) +.00000000000000079608L * powl(z,39.0L) -.00000000000000009420L * powl(z,41.0L) -.00000000000000000713L * powl(z,43.0L) +.00000000000000000033L * powl(z,45.0L) +.00000000000000000004L * powl(z,47.0L)); else return(+.00046483389361763382L * powl(z, 0.0L) -.00100566073653404708L * powl(z, 2.0L) +.00024044856573725793L * powl(z, 4.0L) +.00102830861497023219L * powl(z, 6.0L) -.00076578610717556442L * powl(z, 8.0L) -.00020365286803084818L * powl(z,10.0L) +.00023212290491068728L * powl(z,12.0L) +.00003260214424386520L * powl(z,14.0L) -.00002557906251794953L * powl(z,16.0L) -.00000410746443891574L * powl(z,18.0L) +.00000117811136403713L * powl(z,20.0L) +.00000024456561422485L * powl(z,22.0L) -.00000002391582476734L * powl(z,24.0L) -.00000000750521420704L * powl(z,26.0L) +.00000000013312279416L * powl(z,28.0L) +.00000000013440626754L * powl(z,30.0L) +.00000000000351377004L * powl(z,32.0L) -.00000000000151915445L * powl(z,34.0L) -.00000000000008915418L * powl(z,36.0L) +.00000000000001119589L * powl(z,38.0L) +.00000000000000105160L * powl(z,40.0L) -.00000000000000005179L * powl(z,42.0L) -.00000000000000000807L * powl(z,44.0L) +.00000000000000000011L * powl(z,46.0L) +.00000000000000000004L * powl(z,48.0L)); } /**************************************************************************/ /**************************************************************************/