#include #include /*************************************************************************** * * * Program to find the adjustments h_m used in Turing's method for bounding * * S(g_n)= N(g_n) - theta(g_n) / pi - 1. * * * * For g_END, the program finds values of h_(END-SAMPLES+1), ..., h_END * * such that * * * * 1.) g_(END-SAMPLES+1)+h_(END-SAMPLES+1), ..., g_END+h_END is * * strictly increasing, and * * * * 2.) (-1)^n*Z(g_n+h_n) > 0 for n = END-SAMPLES+1, ..., END. * * * * Output to the RESULTFILE consists of n, g_n, h_n, g_n+h_n, Z(g_n+h_n), * * for n = END-SAMPLES+1, ..., END. * * * * on RedHat Linux compile with * * * * gcc -lm -o turingtest turingtest2.c * * * * and run with * * * * ./turingtest * * * ***************************************************************************/ #define END 12193892 /* final graham point */ #define SAMPLES 20 /* number of samples to compute */ #define RESULTFILE "turing.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 */ long double g[SAMPLES]; /* list of Graham points */ long double h[SAMPLES]; /* list of adjustment terms */ long double h_total; /* total of adjustment terms */ long double step; /* step size for adjustment terms */ int g_initial; /* index of first Graham point in range */ int even(int); /* -1,1 parity function */ int num_samples; /* number of samples for test */ int j; /* loop index */ int k; /* loop index */ int n; /* loop index */ int m; /* loop index */ num_samples = SAMPLES; /* */ j = 0; /* */ k = 0; /* */ m = 0; /* initializations */ g_initial = END - num_samples + 1; /* */ step = 0.1L; /* */ h_total = 0.0L; /* */ results_file = fopen(RESULTFILE, "w"); /* open output file */ while (j < num_samples) /* evaluate Graham pts */ { /* in range; */ g[j] = graham(g_initial+j); /* initialize */ h[j] = 0.0L; /* adjustment vector */ ++j; /* */ } /* */ while (k < num_samples) /* for each Graham */ { /* point ... */ n = 0; /* set no. of steps to 0*/ while (even(g_initial+k)*Z(g[k]-n*step,4) /* keep stepping one */ < 0.053L/powl(g[k]-n*step,1.25L) /* step to both the */ && /* left and right of */ even(g_initial+k)*Z(g[k]+n*step,4) /* g_n until */ < 0.053L/powl(g[k]+n*step,1.25L)) /* (-1)^n*Z(g_n+h_n)>0 */ { /* is satisfied, where */ ++n; /* h_n = +/- n*step */ } /* */ if (even(g_initial+k)*Z(g[k]-n*step,4) /* store the minimal */ > 0.053L/powl(g[k]-n*step,1.25L)) /* adjustment found in */ h[k] = -n*step; /* in the h[] vector */ else /* */ h[k] = n*step; /* */ h_total += h[k]; /* increase adj. total */ if ((k > 0) && (g[k]+h[k]<=g[k-1]+h[k-1])) /* if the sequence */ { /* g_k+h_k so far is */ h_total = 0.0L; /* not strictly */ step = step/2.0L; /* increasing, reduce */ k = 0; /* the step size and */ } /* start over, otherwise*/ else ++k; /* go to next Graham pt.*/ } while (m < num_samples) { /* write out the */ fprintf(results_file, /* results */ "%d & " /* */ "%11.3Lf & " /* */ "%5.3Lf & " /* */ "%11.3Lf & " /* */ "%5.3Lf \\\\ \n", /* */ g_initial+m, /* */ g[m], /* */ h[m], /* */ g[m]+h[m], /* */ Z(g[m]+h[m],4)); /* */ ++m; /* */ } /* */ fprintf(results_file,"h_total = %5.3Lf\n", /* */ h_total); /* */ fclose(results_file); /* close output file */ printf("%20.10Lf,%20.10Lf\n",graham(1181229), graham(1181235)); 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); } 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 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; ZZ = 0.0L; R = 0.0L; j = 1; k = 0; N = sqrtl(t/(2.0L * pi)); p = sqrtl(t/(2.0L * pi)) - N; while (j <= N) { ZZ = ZZ + 1.0L/sqrtl((long double) j ) * cosl(fmodl(theta(t)-t*logl((long double) j),2.0L*pi)); ++j; } ZZ = 2.0L * ZZ; while (k <= n) { 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 approximation */ const long double pi = 3.1415926535897932385L; t_n=0.0L; t_n1=0.5655L*n+855.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)); } /**************************************************************************/