/* Copyright (C) 2021 Alessandro Languasco Computation of the following functions using the algorithms in the paper "Efficient computation some special functions" All these programs require in input a file named primroot.res containing: on the first line: a prime number q; on the second line: a primitive root g of q. 1) zetah(s, k1, k2), s>1 (for FFT applications; results saved on file); computes zeta(s,(g^k mod q)/q) for k1<=k<=k2; 2) zetahprime(s, k1, k2), s>1 (for FFT applications; results saved on file) computes the Hurwitz zeta function zetaprime(s,(g^k mod q)/q) for k1<=k<=k2; 3) zetahtotal(s, k1, k2), s>1 (for FFT applications; results saved on file) computes the first derivative of the Hurwitz zeta function zeta(s,(g^k mod q)/q) and zetaprime(s,(g^k mod q)/q) for k1<=k<=k2;; 4) hurwitz_DIF_even (compute both zetah and zetahprime with reflection fornulae; even case; propositions 5 and 6 of the paper; results saved on file) computes zeta(s,(g^k mod q)/q) + zeta(s,1-(g^k mod q)/q) for k1<=k<=k2 and zetaprime(s,(g^k mod q)/q) + zetaprime(s,1-(g^k mod q)/q) for k1<=k<=k2 5) hurwitz_DIF_odd (compute both zetah and zetahprime with reflection fornulae; odd case; propositions 5 and 6 of the paper; results saved on file) computes zeta(s,(g^k mod q)/q) - zeta(s,1-(g^k mod q)/q) for k1<=k<=k2 and zetaprime(s,(g^k mod q)/q) - zetaprime(s,1-(g^k mod q)/q) for k1<=k<=k2. REMARKS: 1) the values of the Riemann zeta function, of zetaprime, of Euler's Beta function, and of the differences of the digamma function are precomputed at runtime each time the main functions are called; another important gain in the running time would follow from having such values stored once for all (for example, they could be computed and stored when the first time this package is called) 2) The precision parameter is fixed to 128 bits. 3) toward the bottom of this file you will also find some examples about how to use these functions and a comparison of their running times with the analogue ones obtained with the internal PARI/GP functions */ /***************** Precomputations for the external C program ****************/ \\ \p needed to decide the precision default (format, f) for the floating point \\ Hurwitz zeta function computed in a/q, q prime, 1 <= a <= q-1 {zetah(s,x,y)= local(vec,boundindex, aoverq, S, oneminusaoverq, minusaoverq, q, gammaratiosvector, zetas, maxr1, aux, g, a, u , coeff, r1, m, maxindex, ceils, ceilsplusone, zetaRvector, minutes, millisec, abslogaoverq, abslogoneminusaoverq, aux1, magnitude, decprecision, seconds, elaptimeprecomp, zetahfile, name, namefile, start, correction, coeffvector); default(format , f); \\ print numbers with floating point notation (not E notation) vec = readvec("./primroot.res"); \\getting q and g from primroot.res q=vec[1]; g=vec[2]; print("Hurwitz-zeta-values for q = ",q, " and saving on file"); if (s <= 1, error("s >1 !!!!")); if (s-floor(s) == 0, s=floor(s)); if (x<0, x=0); \\ needed to avoid problem with the first interval if (x>=y, error("first exponent should be less than the second")); if (y>q-2, y=q-2); \\ needed to avoid problem with the last interval name="./precomp_zetah"; namefile = Strprintf("%s-%012d-%012d.res", name, x, y); zetahfile = fileopen(namefile, "w"); if (x==0,filewrite(zetahfile,q)); \\ print(q);print(g); magnitude = ceil ( s* abs(log(q)/log(10)) ) ; \\ maximal order of magnitude; attained at 1/q magnitude = ceil(1.1*magnitude); m=128; \\ binary precision required for our algorithm \\ +5 to be safer; decprecision = magnitude + ceil( log(2^m)/ log(10) ) + 5; default( realprecision, decprecision ); print("Internal precision: # decimal digits: ", decprecision); coeff = (m+1)*log(2); zetas = zeta(s); \\ valueatonehalf =(2^s-1)*zetas; \\ q is prime -> a/q<>1/2 ceils = ceil(s); ceilsplusone = ceils+1; aux1 = lngamma(ceilsplusone+s)-lngamma(s)-lngamma(ceilsplusone+1)+abs(log(zeta(ceilsplusone+s)-1)); boundindex = aux1 / log(2) ; boundindex = ceil(boundindex); maxindex = m + 10 + boundindex; zetaRvector = vector(maxindex,i, 0) ; zetaRvector[1] = zetas - 1; for ( i = 2, maxindex, zetaRvector[i] = zeta(s+i-1)-1); gammaratiosvector = vector(maxindex,i, 0); gammaratiosvector[1] = s; for ( i = 2, maxindex, gammaratiosvector[i] = gammaratiosvector[i-1]*(s+i-1)/i ); coeffvector = vector(maxindex,i, 0); for ( i = 1, maxindex-1, coeffvector[i] = zetaRvector[i+1] * gammaratiosvector[i] ); gettime(); maxr1 = 0; a=g^x%q; aux = coeff + abs( log ( gammaratiosvector[ceilsplusone] * zetaRvector[ceilsplusone+1] ) ); for(k=x, y, aoverq = (a/q)*1.0; \\ already sorted for the final summation in external functionx S=0; \\ \\if (aoverq == 0.5, print("1/2"); S = valueatonehalf, \\ q is prime => a/q <> 1/2 \\ a/q <>1/2 minusaoverq = -aoverq; oneminusaoverq = 1 - aoverq; abslogaoverq = abs(log(aoverq)); \\ da controllare abslogoneminusaoverq = abs(log(oneminusaoverq)); \\ da controllare if (aoverq > 0.5, \\ a/q >1/2 r1 = ceil ( (aux + abslogaoverq)/ abslogoneminusaoverq); \\ number of summands if(r1 < 2, r1= 2); if(r1 < ceils, r1 = ceils); start = oneminusaoverq; correction = 1/(aoverq)^s + zetaRvector[1] , \\ a/q <1/2 r1 = ceil ( (aux + abslogoneminusaoverq)/ abslogaoverq); \\ number of summands if(r1 < 2, r1= 2); if(r1 < ceils, r1 = ceils); start = minusaoverq; correction = 1/(aoverq)^s + 1/(1+aoverq)^s + zetaRvector[1] ); if (r1>maxr1, maxr1 = r1); u = start; \\ starting point for computing using a repeated product strategy for (k=1, r1, \\S+= u * zetaRvector[k+1] * gammaratiosvector[k]; S+= u * coeffvector[k]; u*= start; ); S = correction + S; \\ correction values\\ zeta(s) is the zero term filewrite(zetahfile, S); a=(a*g)%q; \\a=g^k%q; ); elaptimeprecomp=gettime(); fileflush(zetahfile); fileclose(zetahfile); /* print computation time */ seconds=floor(elaptimeprecomp/1000)%60; minutes=floor(elaptimeprecomp/60000); millisec=elaptimeprecomp- minutes*60000 - seconds*1000; print("Hurwitz-zeta-values computation time (I/O included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); } \\ First derivative of the Hurwitz zeta function computed in a/q, q prime, 1 <= a <= q-1 {zetahprime(s,x,y)= local(vec,boundindex, aoverq, S, oneminusaoverq, minusaoverq, q, psivector, zetas, maxr1, aux, g, a, u , coeff, r1, m, maxindex, ceils, ceilsplusone, zetaRvector, minutes, millisec, abslogaoverq, abslogoneminusaoverq, aux1, magnitude, decprecision, seconds, elaptimeprecomp, zetahprimefile, name, namefile, start, correction, coeffprimevector, betacoeff, zetaRprimevector, oneplusaoverq, logoneplusaoverq, logaoverq, gammaratiosvector); default(format , f); \\ print numbers with floating point notation (not E notation) vec = readvec("./primroot.res"); \\getting q and g from primroot.res q=vec[1]; g=vec[2]; print("Hurwitz-zeta-prime-values for q = ",q, " and saving on file"); if (s <= 1, error("s >1 !!!!")); if (s-floor(s) == 0, s=floor(s)); if (x<0, x=0); \\ needed to avoid problem with the first interval if (x>=y, error("first exponent should be less than the second")); if (y>q-2, y=q-2); \\ needed to avoid problem with the last interval name="./precomp_zetah_prime"; namefile = Strprintf("%s-%012d-%012d.res", name, x, y); zetahprimefile = fileopen(namefile, "w"); if (x==0,filewrite(zetahprimefile,q)); magnitude = ceil ( s* abs(log(q)/log(10)) ) ; \\ maximal order of magnitude; attained at 1/q magnitude = ceil(1.1*magnitude); m=128; \\ binary precision required for our algorithm \\ +5 to be safer; decprecision = magnitude + ceil( log(2^m)/ log(10) ) + 5; default( realprecision, decprecision ); print("Internal precision: # decimal digits: ", decprecision); ceils = ceil(s); ceilsplusone = ceils+1; coeff = (m+1+2*ceils)*log(2); zetas = zeta(s); \\valueatonehalf = 2^s*log(2)*zetas +(2^s-1)*zetahurwitz(s,1,1); \\ q is prime -> a/q<>1/2 betacoeff = abs(lngamma(ceilsplusone+s)-lngamma(s)-lngamma(ceilsplusone+1)); aux1 = betacoeff + abs(log(zeta(ceilsplusone+s)-1)); boundindex = aux1 / log(2) ; boundindex = ceil(boundindex); maxindex = m + 10 + boundindex; zetaRvector = vector(maxindex,i, 0) ; zetaRvector[1] = zetas - 1; for ( i = 2, maxindex, zetaRvector[i] = zeta(s+i-1)-1); zetaRprimevector = vector(maxindex,i, 0) ; zetaRprimevector[1] = zetahurwitz(s,1,1); \\ zeta'(s) per i=1 for ( i = 2, maxindex, zetaRprimevector[i] = zetahurwitz(s+i-1,1,1) ); \\ zeta'(s+i-1) gammaratiosvector = vector(maxindex,i, 0) ; gammaratiosvector[1] = s; for ( i = 2, maxindex, gammaratiosvector[i] = gammaratiosvector[i-1]*(s+i-1)/i ); psivector = vector(maxindex,i, 0) ; \\initialisation of the needed coefficients psivector[1] = 1/s; \\psi(s+1) - psi(s) for ( i = 2, maxindex, psivector[i] = psivector[i-1] +1/(s+i-1) ); \\ psi(s+i) - psi(s) coeffprimevector = vector(maxindex,i, 0); for ( i = 1, maxindex-1, coeffprimevector[i] = (zetaRvector[i+1]*psivector[i]+zetaRprimevector[i+1] ) * gammaratiosvector[i] ); gettime(); maxr1 = 0; a=g^x%q; aux = coeff + betacoeff; for(k=x, y, aoverq = (a/q)*1.0; \\ already sorted for the final summation in external function S=0; \\ \\if (aoverq == 0.5, print("1/2"); S = valueatonehalf, \\ q is prime => a/q <> 1/2 \\ a/q <>1/2 minusaoverq = -aoverq; oneminusaoverq = 1 - aoverq; oneplusaoverq = 1 + aoverq; logoneplusaoverq = log(oneplusaoverq); logaoverq = log(aoverq); abslogaoverq = abs(logaoverq); abslogoneminusaoverq = abs(log(oneminusaoverq)); if (aoverq > 0.5, \\ a/q >1/2 r1 = ceil ( (aux + abslogaoverq)/ abslogoneminusaoverq); \\ number of summands if(r1 < 2, r1= 2); if(r1 < ceils, r1 = ceils); start = oneminusaoverq; correction = -logaoverq/(aoverq)^s + zetaRprimevector[1] , \\ a/q <1/2 r1 = ceil ( (aux + abslogoneminusaoverq)/ abslogaoverq); \\ number of summands if(r1 < 2, r1= 2); if(r1 < ceils, r1 = ceils); start = minusaoverq; correction = -logaoverq/(aoverq)^s - logoneplusaoverq/(oneplusaoverq)^s + zetaRprimevector[1] ); if (r1>maxr1, maxr1 = r1); u = start; \\ starting point for computing using a repeated product strategy for (k=1, r1, \\S+= u * (zetaRvector[k+1]*psivector[k]+zetaRprimevector[k+1] ) * gammaratiosvector[k]; S+= u *coeffprimevector[k]; u*= start; ); S = correction + S; \\ correction values\\ zetaprime(s) is the zero term filewrite(zetahprimefile, S); a=(a*g)%q; \\a=g^k%q; ); elaptimeprecomp=gettime(); fileflush(zetahprimefile); fileclose(zetahprimefile); /* print computation time */ seconds=floor(elaptimeprecomp/1000)%60; minutes=floor(elaptimeprecomp/60000); millisec=elaptimeprecomp- minutes*60000 - seconds*1000; \\print(elaptimeprecomp); print("Hurwitz-zeta-prime-values computation time (I/O included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); } \\ computes the Hurwitz zeta function and its first derivative at a/q, q prime, 1 <= a <= q-1 {zetahtotal(s,x,y)= local(vec,boundindex, aoverq, S, oneminusaoverq, minusaoverq, q, psivector, zetas, maxr1, aux, g, a, u , coeff, r1, m, maxindex, ceils, ceilsplusone, zetaRvector, minutes, millisec, abslogaoverq, abslogoneminusaoverq, aux1, magnitude, decprecision, seconds, elaptimeprecomp, zetahprimefile, name, namefile, start, correction, coeffprimevector, betacoeff, zetaRprimevector, oneplusaoverq, logoneplusaoverq, logaoverq, gammaratiosvector, coeffvector, powaqsrec, powoneplusaqsrec); default(format , f); \\ print numbers with floating point notation (not E notation) vec = readvec("./primroot.res"); \\getting q and g from primroot.res q=vec[1]; g=vec[2]; print("Hurwitz-zeta-prime-values for q = ",q, " and saving on file"); if (s <= 1, error("s >1 !!!!")); if (s-floor(s) == 0, s=floor(s)); if (x<0, x=0); \\ needed to avoid problem with the first interval if (x>=y, error("first exponent should be less than the second")); if (y>q-2, y=q-2); \\ needed to avoid problem with the last interval name="./precomp_zetah_prime"; namefile = Strprintf("%s-%012d-%012d.res", name, x, y); zetahprimefile = fileopen(namefile, "w"); if (x==0,filewrite(zetahprimefile,q)); name="./precomp_zetah"; namefile = Strprintf("%s-%012d-%012d.res", name, x, y); zetahfile = fileopen(namefile, "w"); if (x==0,filewrite(zetahfile,q)); magnitude = ceil ( s* abs(log(q)/log(10)) ) ; \\ maximal order of magnitude; attained at 1/q magnitude = ceil(1.1*magnitude); m=128; \\ binary precision required for our algorithm \\ +5 to be safer; decprecision = magnitude + ceil( log(2^m)/ log(10) ) + 5; default( realprecision, decprecision ); print("Internal precision: # decimal digits: ", decprecision); ceils = ceil(s); ceilsplusone = ceils+1; coeff = (m+1+2*ceils)*log(2); zetas = zeta(s); \\valueatonehalf =(2^s-1)*zetas; \\ q is prime -> a/q<>1/2 \\valueprimeatonehalf =2^s*log(2)*zetas +(2^s-1)*zetahurwitz(s,1,1); \\ q is prime -> a/q<>1/2 betacoeff = abs(lngamma(ceilsplusone+s)-lngamma(s)-lngamma(ceilsplusone+1)); aux1 = betacoeff + abs(log(zeta(ceilsplusone+s)-1)); boundindex = aux1 / log(2) ; boundindex = ceil(boundindex); maxindex = m + 10 + boundindex; zetaRvector = vector(maxindex,i, 0) ; zetaRvector[1] = zetas - 1; for ( i = 2, maxindex, zetaRvector[i] = zeta(s+i-1)-1); zetaRprimevector = vector(maxindex,i, 0) ; zetaRprimevector[1] = zetahurwitz(s,1,1); \\ zeta'(s) per i=1 for ( i = 2, maxindex, zetaRprimevector[i] = zetahurwitz(s+i-1,1,1) ); \\ zeta'(s+i-1) gammaratiosvector = vector(maxindex,i, 0) ; gammaratiosvector[1] = s; for ( i = 2, maxindex, gammaratiosvector[i] = gammaratiosvector[i-1]*(s+i-1)/i ); psivector = vector(maxindex,i, 0) ; \\initialisation of the needed coefficients psivector[1] = 1/s; \\psi(s+1) - psi(s) for ( i = 2, maxindex, psivector[i] = psivector[i-1] +1/(s+i-1) ); \\ psi(s+i) - psi(s) coeffvector = vector(maxindex,i, 0); for ( i = 1, maxindex-1, coeffvector[i] = zetaRvector[i+1] * gammaratiosvector[i] ); coeffprimevector = vector(maxindex,i, 0); for ( i = 1, maxindex-1, coeffprimevector[i] = (zetaRvector[i+1]*psivector[i]+zetaRprimevector[i+1] ) * gammaratiosvector[i] ); gettime(); maxr1 = 0; a=g^x%q; aux = coeff + betacoeff; for(k=x, y, aoverq = (a/q)*1.0; \\ already sorted for the final summation in external function Sprime=0; S=0; \\ \\if (aoverq == 0.5, print("1/2"); S = valueatonehalf, \\ q is prime => a/q <> 1/2 \\ a/q <>1/2 minusaoverq = -aoverq; oneminusaoverq = 1 - aoverq; oneplusaoverq = 1 + aoverq; logoneplusaoverq = log(oneplusaoverq); logaoverq = log(aoverq); abslogaoverq = abs(logaoverq); abslogoneminusaoverq = abs(log(oneminusaoverq)); powaqsrec= 1/(aoverq)^s; powoneplusaqsrec = 1/(oneplusaoverq)^s; if (aoverq > 0.5, \\ a/q >1/2 r1 = ceil ( (aux + abslogaoverq)/ abslogoneminusaoverq); \\ number of summands if(r1 < 2, r1= 2); if(r1 < ceils, r1 = ceils); start = oneminusaoverq; correctionprime = -logaoverq*powaqsrec + zetaRprimevector[1]; correction = powaqsrec + zetaRvector[1] , \\ a/q <1/2 r1 = ceil ( (aux + abslogoneminusaoverq)/ abslogaoverq); \\ number of summands if(r1 < 2, r1= 2); if(r1 < ceils, r1 = ceils); start = minusaoverq; correctionprime = -logaoverq*powaqsrec - logoneplusaoverq*powoneplusaqsrec + zetaRprimevector[1]; correction = powaqsrec + powoneplusaqsrec + zetaRvector[1] ); if (r1>maxr1, maxr1 = r1); u = start; \\ starting point for computing using a repeated product strategy for (k=1, r1, Sprime+= u *coeffprimevector[k]; S+= u * coeffvector[k]; u*= start; ); Sprime = correctionprime + Sprime; \\ correction values\\ zetaprime(s) is the zero term S = correction + S; \\ correction values\\ zeta(s) is the zero term filewrite(zetahprimefile, Sprime); filewrite(zetahfile, S); a=(a*g)%q; \\a=g^k%q; ); elaptimeprecomp=gettime(); fileflush(zetahprimefile); fileclose(zetahprimefile); fileflush(zetahfile); fileclose(zetahfile); /* print computation time */ seconds=floor(elaptimeprecomp/1000)%60; minutes=floor(elaptimeprecomp/60000); millisec=elaptimeprecomp- minutes*60000 - seconds*1000; \\print(elaptimeprecomp); print("Hurwitz-values computation time (I/O included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); } /* DECIMATION in Frequency - odd case - both zetah and zetahprime in a/q */ {hurwitz_DIF_odd(s,x,y)= local(vec,boundindex, aoverq, S_odd, S_prime_odd, oneminusaoverq, minusaoverq, q, psivector, zetas, maxr1, aux, g, a, u_odd , coeff, r1, m, maxindex, ceils, ceilsplusone, fattore_odd, segno_odd, zetaRvector, minutes, millisec, abslogaoverq, abslogoneminusaoverq, aux1, magnitude, decprecision, twominusaoverq, logoneminusaoverq, logtwominusaoverq, correction_prime_odd, seconds, elaptimeprecomp, zetahfile_odd, zetahprimefile_odd, name_odd, namefile_odd, name_prime_odd, namefile_prime_odd, start_odd, correction_odd, pow1, pow2, pow3, pow4, coeffvector, coeffprimevector, betacoeff, zetaRprimevector, oneplusaoverq, logoneplusaoverq, logaoverq, gammaratiosvector); default(format , f); \\ print numbers with floating point notation (not E notation) vec = readvec("./primroot.res"); \\getting q and g from primroot.res q=vec[1]; g=vec[2]; print("ODD DIF CASE Hurwitz-zeta-zetaprime-values for q = ",q, " and saving on file"); if (s <= 1, error("s >1 !!!!")); if (s-floor(s) == 0, s=floor(s)); if (x<0, x=0); \\ needed to avoid problem with the first interval if (x>=y, error("first exponent should be less than the second")); if (y>q-2, y=q-2); \\ needed to avoid problem with the last interval name_odd="./precomp_zetah_DIFodd"; namefile_odd = Strprintf("%s-%012d-%012d.res", name_odd, x, y); name_prime_odd="./precomp_zetah_primeDIFodd"; namefile_prime_odd = Strprintf("%s-%012d-%012d.res", name_prime_odd, x, y); zetahfile_odd = fileopen(namefile_odd, "w"); zetahprimefile_odd = fileopen(namefile_prime_odd, "w"); if (x==0,filewrite(zetahfile_odd,q);filewrite(zetahprimefile_odd,q)); magnitude = ceil ( s* abs(log(q)/log(10)) ) ; \\ maximal order of magnitude; attained at 1/q magnitude = ceil(1.1*magnitude); m=128; \\ binary precision required for our algorithm \\ +5 to be safer; decprecision = magnitude + ceil( log(2^m)/ log(10) ) + 5; default( realprecision, decprecision ); print("Internal precision: # decimal digits: ", decprecision); ceils = ceil(s); ceilsplusone = ceils+1; coeff = (m+1+2*ceils)*log(2); zetas = zeta(s); \\ valueatonehalf =(2^s-1)*zetas; \\ q is prime -> a/q<>1/2 \\valueprimeatonehalf =2^s*log(2)*zetas +(2^s-1)*zetahurwitz(s,1,1); \\ q is prime -> a/q<>1/2 betacoeff = abs(lngamma(ceilsplusone+s)-lngamma(s)-lngamma(ceilsplusone+1)); aux1 = betacoeff + abs(log(zeta(ceilsplusone+s)-1)); boundindex = aux1 / log(2) ; boundindex = ceil(boundindex); maxindex = m + 10 + boundindex; maxindex = ceil(maxindex/2); zetaRvector = vector(maxindex,i, 0) ; zetaRvector[1] = zetas - 1; for ( i = 2, maxindex, zetaRvector[i] = zeta(s+2*(i-1)-1)-1); zetaRprimevector = vector(maxindex,i, 0) ; zetaRprimevector[1] = zetahurwitz(s,1,1); \\ zeta'(s) for i=1 for ( i = 2, maxindex, zetaRprimevector[i] = zetahurwitz(s+2*(i-1)-1,1,1) ); \\ zeta'(s+2*i-1) gammaratiosvector = vector(maxindex,i, 0) ; gammaratiosvector[1] = s; for ( i = 2, maxindex, gammaratiosvector[i] = gammaratiosvector[i-1]*(s+2*i-2)/(2*i-1)*(s+2*i-3)/(2*i-2) ); psivector = vector(maxindex,i, 0) ; psivector[1] = 1/s; \\psi(s+1) - psi(s) for ( i = 2, maxindex, psivector[i] = psivector[i-1] + 1/(s+2*i-3) + 1/(s+2*i-2) ); \\ psi(s+2*i-1) - psi(s) coeffvector = vector(maxindex,i, 0); for ( i = 1, maxindex-1, coeffvector[i] = zetaRvector[i+1] * gammaratiosvector[i] ); coeffprimevector = vector(maxindex,i, 0); for ( i = 1, maxindex-1, coeffprimevector[i] = (zetaRvector[i+1]*psivector[i]+zetaRprimevector[i+1] ) * gammaratiosvector[i] ); gettime(); maxr1 = 0; a=g^x%q; aux = coeff + betacoeff; for(k=x, y, aoverq = (a/q)*1.0; \\ already sorted for the final summation in external function S_odd=0; S_prime_odd=0; \\ \\if (aoverq == 0.5, print("1/2"); S = valueatonehalf, \\ q is prime => a/q <> 1/2 \\ a/q <>1/2 minusaoverq = -aoverq; oneminusaoverq = 1 + minusaoverq; twominusaoverq = 2 + minusaoverq; oneplusaoverq = 1 + aoverq; logoneplusaoverq = log(oneplusaoverq); logoneminusaoverq = log(oneminusaoverq); logtwominusaoverq = log(twominusaoverq); logaoverq = log(aoverq); abslogaoverq = abs(logaoverq); abslogoneminusaoverq = abs(logoneminusaoverq); pow1 = 1/(aoverq)^s; pow2 = 1/(oneminusaoverq)^s; pow3 = 1/(twominusaoverq)^s; pow4 = 1/(oneplusaoverq)^s; if (aoverq > 0.5, \\ a/q >1/2 r1 = ceil ( (aux + abslogaoverq)/ abslogoneminusaoverq); \\ number of summands if(r1 < 2, r1= 2); if(r1 < ceils, r1 = ceils); start_odd = oneminusaoverq; fattore_odd = oneminusaoverq^2; segno_odd = 1.; correction_odd = pow1 - pow2 - pow3; correction_prime_odd = -logaoverq * pow1 + logoneminusaoverq * pow2 + logtwominusaoverq * pow3 , \\ a/q <1/2 r1 = ceil ( (aux + abslogoneminusaoverq)/ abslogaoverq); \\ number of summands if(r1 < 2, r1= 2); if(r1 < ceils, r1 = ceils); start_odd = aoverq; fattore_odd = aoverq^2; segno_odd = -1.; correction_odd = pow1 + pow4 - pow2; correction_prime_odd = -logaoverq * pow1 + logoneminusaoverq * pow2 - logoneplusaoverq * pow4 ); if (r1>maxr1, maxr1 = r1); r1= ceil(r1/2); \\ starting point for computing using a repeated product strategy u_odd= start_odd; for (k=1, r1, \\S_odd += u_odd * zetaRvector[k+1] * gammaratiosvector[k]; \\S_prime_odd += u_odd * (zetaRvector[k+1]*psivector[k]+zetaRprimevector[k+1] ) * gammaratiosvector[k]; S_odd += u_odd * coeffvector[k]; S_prime_odd += u_odd *coeffprimevector[k]; u_odd *= fattore_odd; ); S_odd = correction_odd + segno_odd*2*S_odd; S_prime_odd = correction_prime_odd + segno_odd*2*S_prime_odd; \\ correction values\\ zeta(s) is the zero term filewrite(zetahfile_odd, S_odd); filewrite(zetahprimefile_odd, S_prime_odd); a=(a*g)%q; \\a=g^k%q; ); \\print(maxr1); \\print(maxindex); elaptimeprecomp=gettime(); fileflush(zetahfile_odd); fileclose(zetahfile_odd); fileflush(zetahprimefile_odd); fileclose(zetahprimefile_odd); /* print computation time */ seconds=floor(elaptimeprecomp/1000)%60; minutes=floor(elaptimeprecomp/60000); millisec=elaptimeprecomp- minutes*60000 - seconds*1000; \\print(elaptimeprecomp); print("ODD DIF CASE Hurwitz-zeta-zetaprime-values computation time (I/O included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); } /* DECIMATION in Frequency - even case - both zetah and zetahprime in a/q */ {hurwitz_DIF_even(s,x,y)= local(vec,boundindex, aoverq, S_even, S_prime_even, oneminusaoverq, minusaoverq, q, psivector, zetas, maxr1, aux, g, a, u_even , coeff, r1, m, maxindex, ceils, ceilsplusone, fattore_even, segno_even, zetaRvector, minutes, millisec, abslogaoverq, abslogoneminusaoverq, aux1, magnitude, decprecision, twominusaoverq, logoneminusaoverq, logtwominusaoverq, correction_prime_even, seconds, elaptimeprecomp, zetahfile_even, zetahprimefile_even, name_even, namefile_even, name_prime_even, namefile_prime_even, start_even, correction_even, pow1, pow2, pow3, pow4, coeffvector, coeffprimevector, betacoeff, zetaRprimevector, oneplusaoverq, logoneplusaoverq, logaoverq, gammaratiosvector); default(format , f); \\ print numbers with floating point notation (not E notation) vec = readvec("./primroot.res"); \\getting q and g from primroot.res q=vec[1]; g=vec[2]; print("EVEN DIF CASE Hurwitz-zeta-zetaprime-values for q = ",q, " and saving on file"); if (s <= 1, error("s >1 !!!!")); if (s-floor(s) == 0, s=floor(s)); if (x<0, x=0); \\ needed to avoid problem with the first interval if (x>=y, error("first exponent should be less than the second")); if (y>q-2, y=q-2); \\ needed to avoid problem with the last interval name_even="./precomp_zetah_DIFeven"; namefile_even = Strprintf("%s-%012d-%012d.res", name_even, x, y); name_prime_even="./precomp_zetah_primeDIFeven"; namefile_prime_even = Strprintf("%s-%012d-%012d.res", name_prime_even, x, y); zetahfile_even = fileopen(namefile_even, "w"); zetahprimefile_even = fileopen(namefile_prime_even, "w"); if (x==0,filewrite(zetahfile_even,q);filewrite(zetahprimefile_even,q)); \\ print(q);print(g); magnitude = ceil ( s* abs(log(q)/log(10)) ) ; \\ maximal order of magnitude; attained at 1/q magnitude = ceil(1.1*magnitude); m=128; \\ binary precision required for our algorithm \\ +5 to be safer; decprecision = magnitude + ceil( log(2^m)/ log(10) ) + 5; default( realprecision, decprecision ); print("Internal precision: # decimal digits: ", decprecision); ceils = ceil(s); ceilsplusone = ceils+1; coeff = (m+1+2*ceils)*log(2); zetas = zeta(s); \\ valueatonehalf =(2^s-1)*zetas; \\ q is prime; a/q<>1/2 \\valueprimeatonehalf =2^s*log(2)*zetas +(2^s-1)*zetahurwitz(s,1,1); \\ q is prime -> a/q<>1/2 betacoeff = abs(lngamma(ceilsplusone+s)-lngamma(s)-lngamma(ceilsplusone+1)); aux1 = betacoeff + abs(log(zeta(ceilsplusone+s)-1)); boundindex = aux1 / log(2) ; boundindex = ceil(boundindex); maxindex = m + 10 + boundindex; maxindex = ceil(maxindex/2); zetaRvector = vector(maxindex,i, 0) ; zetaRvector[1] = zetas - 1; for ( i = 2, maxindex, zetaRvector[i] = zeta(s+2*(i-1))-1); zetaRprimevector = vector(maxindex,i, 0) ; zetaRprimevector[1] = zetahurwitz(s,1,1); \\ zeta'(s) for i=1 for ( i = 2, maxindex, zetaRprimevector[i] = zetahurwitz(s+2*(i-1),1,1) ); \\ zeta'(s+2i-1) gammaratiosvector = vector(maxindex,i, 0) ; gammaratiosvector[1] = s*(s+1)/2; for ( i = 2, maxindex, gammaratiosvector[i] = gammaratiosvector[i-1]*(s+2*i-1)/(2*i)*(s+2*i-2)/(2*i-1) ); psivector = vector(maxindex,i, 0) ; psivector[1] = 1/s + 1/(s+1); \\psi(s+2) - psi(s) for ( i = 2, maxindex, psivector[i] = psivector[i-1] + 1/(s+2*i-1) + 1/(s+2*i-2) ); \\ psi(s+2*i) - psi(s) coeffvector = vector(maxindex,i, 0); for ( i = 1, maxindex-1, coeffvector[i] = zetaRvector[i+1] * gammaratiosvector[i] ); coeffprimevector = vector(maxindex,i, 0); for ( i = 1, maxindex-1, coeffprimevector[i] = (zetaRvector[i+1]*psivector[i]+zetaRprimevector[i+1] ) * gammaratiosvector[i] ); gettime(); maxr1 = 0; a=g^x%q; aux = coeff + betacoeff; for(k=x, y, aoverq = (a/q)*1.0; \\ already sorted for the final summation in external function S_even=0; S_prime_even=0; \\ \\if (aoverq == 0.5, print("1/2"); S = valueatonehalf, \\ q is prime => a/q <> 1/2 \\ a/q <>1/2 minusaoverq = -aoverq; oneminusaoverq = 1 + minusaoverq; twominusaoverq = 2 + minusaoverq; oneplusaoverq = 1 + aoverq; logoneplusaoverq = log(oneplusaoverq); logoneminusaoverq = log(oneminusaoverq); logtwominusaoverq = log(twominusaoverq); logaoverq = log(aoverq); abslogaoverq = abs(logaoverq); abslogoneminusaoverq = abs(logoneminusaoverq); pow1 = 1/(aoverq)^s; pow2 = 1/(oneminusaoverq)^s; pow3 = 1/(twominusaoverq)^s; pow4 = 1/(oneplusaoverq)^s; if (aoverq > 0.5, \\ a/q >1/2 r1 = ceil ( (aux + abslogaoverq)/ abslogoneminusaoverq); \\ number of summands if(r1 < 2, r1= 2); if(r1 < ceils, r1 = ceils); start_even = oneminusaoverq^2; fattore_even = oneminusaoverq^2; segno_even = 1.; correction_even = pow1 + pow2 + pow3 + 2*zetaRvector[1]; correction_prime_even = -logaoverq * pow1 - logoneminusaoverq * pow2 - logtwominusaoverq * pow3 + 2*zetaRprimevector[1] , \\ a/q <1/2 r1 = ceil ( (aux + abslogoneminusaoverq)/ abslogaoverq); \\ number of summands if(r1 < 2, r1= 2); if(r1 < ceils, r1 = ceils); start_even = aoverq^2; fattore_even = aoverq^2; segno_even = 1.; correction_even = pow1 + pow4 + pow2 + 2*zetaRvector[1]; correction_prime_even = -logaoverq * pow1 - logoneminusaoverq * pow2 - logoneplusaoverq * pow4 + 2*zetaRprimevector[1] ); if (r1>maxr1, maxr1 = r1); r1= ceil(r1/2); \\ starting point for computing using a repeated product strategy u_even= start_even; for (k=1, r1, \\S_even += u_even * zetaRvector[k+1] * gammaratiosvector[k]; \\S_prime_even += u_even * (zetaRvector[k+1]*psivector[k]+zetaRprimevector[k+1] ) * gammaratiosvector[k]; S_even += u_even * coeffvector[k]; S_prime_even += u_even * coeffprimevector[k]; u_even *= fattore_even; ); S_even = correction_even + segno_even*2*S_even; S_prime_even = correction_prime_even + segno_even*2*S_prime_even; \\ correction values\\ zeta(s) is the zero term filewrite(zetahfile_even, S_even); filewrite(zetahprimefile_even, S_prime_even); a=(a*g)%q; \\a=g^k%q; ); elaptimeprecomp=gettime(); fileflush(zetahfile_even); fileclose(zetahfile_even); fileflush(zetahprimefile_even); fileclose(zetahprimefile_even); /* print computation time */ seconds=floor(elaptimeprecomp/1000)%60; minutes=floor(elaptimeprecomp/60000); millisec=elaptimeprecomp- minutes*60000 - seconds*1000; print("EVEN DIF CASE Hurwitz-zeta-zetaprime-values computation time (I/O included): ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); } /************************************ Performances tests performed with the "-nosave" version to remove the time needed to write on a file from the total computation time. ***********/