/*** Computation of C(q,a) A*q=B truncation limit of primes *** *** chi is a character mod q; q >=3 init_MertensConstant2();mertensconstants(positive integer) *******************************************************************************************************/ \\ Global variables: \\ bernvector: vector of the needed Bernoulli numbers \\ factvector: vector of the needed factorials \\ zetavector: vector of the needed values of the Riemann zeta function \\ B: level A*q of the primes used for the finite products \\ U: is the minimun value obtained in computing values of Dirichlet L-functions with the Euler-McLaurin formula \\ isqprime: is 1 if q is prime; 0 otherwise global(bernvector,factvector, zetavector,characters, B, U, isqprime); global(defaulterror=0); {mertensconstants(q)=local(minutes,millisec,seconds, A,K,M,N,T, alphapqa, S,P,phiq,a, u, i, j, maxim, col, freq, maxim1,lines, rows, C,Cresults,Cvector,prodC,deltaprodC,TruncCresults, TruncprodC, chi, chars, finalerr,finalerrvector,totalerr, err1, err2, err4, elaptimeprecomp, elaptimefinalcomp, Schivector, correctdigits, defaultprecision ); \\ Local variables: \\ minutes,millisec,seconds: used just to compute the elapsed computation time \\ A,K,M,N,T: MAIN PARAMETERS; see the paper for their meaning \\ S,P,phiq,a,i,j, maxim, col, freq, maxim1,lines, rows: indexes or aux variables \\ alphapqa: exponend in the Mertens formula for arithmetic progrssion; used in finite products \\ finalerr,finalerrvector,totalerr, err1, err2, err4: errors, see the paper \\ C,Cresults,Cvector,prodC,deltaprodC, TrunCresults, TruncprodC: values od the constants; of their product and its error from the expected value \\ chi, chars : vector of a character mod q; matrix of the characters mod q; see the function genchar \\ defaultprecision: used to fix the precision used in the computations \\ correctdigits: number of correct digits of our output \\ Schivector : with q fixed, computes for every character chi the summation over m and k; doesn't depend on a print("************ A. LANGUASCO & A. ZACCAGNINI *************"); print("****** COMPUTATION OF THE MERTENS' CONSTANTS MOD q; 3 <= q<= 100 *******"); print("****** WITH A PRECISION OF AT LEAST 100 DECIMAL DIGITS *******"); print("q= ",q); \\ primality test on q phiq=eulerphi(q); isqprime=isprime(q); \\ precision setting defaultprecision=120; \\defaultprecision=510; default(realprecision,defaultprecision); defaulterror=10^(-defaultprecision); \\ Initialization of the matrix of characters and their computation chars = matrix(phiq,q+3,lines,rows, 0); chars = genchar(q); \\ print(chars); /********************************************************************************* Start of the Computation of Schivector: this is the most time-consuming part *********************************************************************************/ print("precomputations (quite long); depends on phi(q)."); \\\ Initialization of U that computes the error in estimating the L(s,chi) function using the Euler Mc-Laurin Formula U = 100000000; \\ settings of the main parameters for a precision of about 100 decimal digits A=3200; \\ IMPORTANT: N has to be a multiple of q; T has to be even and positive \\if(q>=3 && q<=10, A=100; K=10; M=10; N=840; T = 6); \\\ faster on the mac laptop if(q>=3 && q<=10, M = 26; K = 26; N = ((840*20\q)+1)*q; T = 88); \\ for the 3-10 range if(q>=11 && q<=20, M = 26; K = 26; N = ((840*25\q)+1)*q; T = 106); \\ for the 11-20 range if(q>=21 && q<=30, M = 26; K = 26; N = ((840*32\q)+1)*q; T = 116); \\ for the 21-30 range if(q>=31 && q<=40, M = 26; K = 26; N = ((840*36\q)+1)*q; T = 128); \\ for the 31-40 range if(q>=41 && q<=50, M = 26; K = 26; N = ((840*40\q)+1)*q; T = 138); \\ for the 41-50 range if(q>=51 && q<=60, M = 26; K = 26; N = ((840*44\q)+1)*q; T = 146); \\ for the 51-60 range if(q>=61 && q<=70, M = 26; K = 26; N = ((840*66\q)+1)*q; T = 128); \\ for the 61-70 range if(q>=71 && q<=80, M = 26; K = 26; N = ((840*70\q)+1)*q; T = 134); \\ for the 71-80 range if(q>=81 && q<=90, M = 26; K = 26; N = ((840*70\q)+1)*q; T = 144); \\ for the 81-90 range if(q>=91 && q<=100, M = 26; K = 26; N = ((840*73\q)+1)*q; T = 152); \\ for the 91-100 range /***** OLD if(q>=61 && q<=70, M = 26; K = 26; N = ((840*48\q)+1)*q; T = 154); \\ for the 61-70 range if(q>=71 && q<=80, M = 26; K = 26; N = ((840*40\q)+1)*q; T = 206); \\ for the 71-80 range if(q>=81 && q<=90, M = 26; K = 26; N = ((840*45\q)+1)*q; T = 206); \\ for the 81-90 range \\ if(q>=91 && q<=100, M = 26; K = 26; N = ((840*49\q)+1)*q; T = 210); \\ for the 91-100 range if(q>=91 && q<=100, M = 26; K = 26; N = ((840*72\q)+1)*q; T = 152); \\ for the 91-100 range ******/ print("Main parameters: for their meaning see the paper"); print("A= ",A); print("K= ",K); print("M= ",M); print("N= ",N); print("T= ",T); gettime(); \\ Initialization of bernvector: a vector of precomputed Bernoulli numbers; \\ bernvector[i] contains the i-1-th Bernoulli number \\ this depends on the fact that Gp indexes start from 1 maxim=max(T,M*K)+1; bernvector=vector(maxim,col,0); factvector=vector(maxim,col,0); for(j=1,maxim,bernvector[j]=bernfrac(j-1); factvector[j]=j!); \\ Initialization of zetavector: a vector of precomputed values of the Riemann zeta function; \\ zetavector[i] contains zeta(i) starting from i=2 \\ the position corresponding to i=1 isn't used maxim1=max(M,K); zetavector=vector(M*K,col,0); freq=vector(M*K,col,0); for(j=2,maxim1,zetavector[j]=zeta(j); freq[j]=freq[j]+1); for(m=2, M, for (k=2,K,zetavector[m*k]=zeta(m*k); freq[m*k]=freq[m*k]+1)); \\ print("Frequency of the computed values of the Riemann zeta", freq); \\ initialization of the vector of the used character and of the vector \\ of the results of the summation over m and k (depends on chi and q; doesn't depend on a) chi = vector(q+3); Schivector= vector(phiq-1,col,0); \\ Computation of the sum over primes depending on a character; \\ The results are stored in a vector Schivector for (j=2, phiq, \\ for j=1 we have the principal character; we don't sum on it, see genchar \\ for(i=1,q, chi[i]=chars[j,i]); \\ just for testing \\ print("******************** j= ",j); \\ just for testing chi=chars[j,]; \\ print("chi = ",chi); \\ just for testing for (m=1,M, \\ print("m= ",m); \\ just for testing Schivector[j-1] = Schivector[j-1] + (m^(-1) * Smchi(m, q, chi, A, K, T, N)) ); ); \\ B is the level of prime used in the finite products B=A*q; if (B>9600, B=9600); write("~/sessionsRIC.tex","************ A. LANGUASCO & A. ZACCAGNINI *************"); write("~/sessionsRIC.tex", " "); write("~/resultsRIC.tex","************ A. LANGUASCO & A. ZACCAGNINI *************"); write("~/resultsshortRIC.tex","************ A. LANGUASCO & A. ZACCAGNINI *************"); write("~/rawresultsRIC.tex","************ A. LANGUASCO & A. ZACCAGNINI *************"); write("~/sessionsRIC.tex","****** COMPUTATION OF THE MERTENS' CONSTANTS MOD q; 3 <= q<= 100 *******"); write("~/sessionsRIC.tex", " "); write("~/resultsRIC.tex","****** COMPUTATION OF THE MERTENS' CONSTANTS MOD q; 3 <= q<= 100 *******"); write("~/resultsshortRIC.tex","****** COMPUTATION OF THE MERTENS' CONSTANTS MOD q; 3 <= q<= 100 *******"); write("~/rawresultsRIC.tex","****** COMPUTATION OF THE MERTENS' CONSTANTS MOD q; 3 <= q<= 100 *******"); write("~/sessionsRIC.tex","****** WITH A PRECISION OF AT LEAST 100 DECIMAL DIGITS *******"); write("~/sessionsRIC.tex", " "); write("~/resultsRIC.tex","****** WITH A PRECISION OF AT LEAST 100 DECIMAL DIGITS *******"); write("~/resultsshortRIC.tex","****** WITH A PRECISION OF AT LEAST 100 DECIMAL DIGITS *******"); write("~/rawresultsRIC.tex","****** WITH A PRECISION OF AT LEAST 100 DECIMAL DIGITS *******"); write("~/sessionsRIC.tex","q= ",q); write("~/sessionsRIC.tex", " "); write("~/sessionsRIC.tex","Main parameters: for their meaning see the paper"); write("~/sessionsRIC.tex", " "); write("~/sessionsRIC.tex","A= ",A); write("~/sessionsRIC.tex", " "); write("~/sessionsRIC.tex","K= ",K); write("~/sessionsRIC.tex", " "); write("~/sessionsRIC.tex","M= ",M); write("~/sessionsRIC.tex", " "); write("~/sessionsRIC.tex","N= ",N); write("~/sessionsRIC.tex", " "); write("~/sessionsRIC.tex","T= ",T); write("~/sessionsRIC.tex", " "); print("Schivector =", Schivector); \\ just for testing write("~/sessionsRIC.tex", "Schivector =", Schivector); write("~/sessionsRIC.tex", " "); print("********** ERRORS **************"); print("Errors: see the paper for their definitions and estimations"); write("~/sessionsRIC.tex","********** ERRORS **************"); write("~/sessionsRIC.tex", " "); write("~/sessionsRIC.tex","Errors: see the paper for their definitions and estimations"); write("~/sessionsRIC.tex", " "); \\ err1 using the estimate |L_Aq|< (Aq)^(-k*m) * min (Aq/(k*m-1); err1= (2*B)/(K^2*(B^K-1)*(B-1)); print("err1= ", err1*1.0); write("~/sessionsRIC.tex","err1= ", err1*1.0); write("~/sessionsRIC.tex", " "); err2= B/(B^M*(B-1)*M*(M-1)); print("err2= ", err2*1.0); write("~/sessionsRIC.tex","err2= ", err2*1.0); write("~/sessionsRIC.tex", " "); print("U = ", U); write("~/sessionsRIC.tex","U = ", U); write("~/sessionsRIC.tex", " "); err4 = (2 * (M * K + T-1)^(T-2) * abs(bernvector[T+1]) * q^T)/(U * factvector[T] * N^(T-1) *(N-1)) ; print("err4= ", err4); write("~/sessionsRIC.tex","err4= ", err4); write("~/sessionsRIC.tex", " "); totalerr= (err1+err2+err4) * (phiq-1)/phiq; print("totalerr= ",totalerr); write("~/sessionsRIC.tex","totalerr= ",totalerr); write("~/sessionsRIC.tex", " "); elaptimeprecomp=gettime(); \\print(elaptimeprecomp); /********************************************************************************** END computations of Schivector ***********************************************************************************/ gettime(); print("used level of primes: B= min(A*q,9600) = ", B); write("~/sessionsRIC.tex","used level of primes: B= min(A*q,9600) = ", B); write("~/sessionsRIC.tex", " "); \\ Initialization of the variables for the results Cvector = vector(phiq); finalerrvector = vector(phiq); Cresults = matrix(phiq,4); TruncCresults= matrix(phiq,4); i=1; for (a=1,q, if (gcd(a,q)==1, \\ Computation of the truncated product over primes P=1; forprime(p=2, B , if (p%q == a, alphapqa=phiq-1, alphapqa=-1); \\print ("p=", p," alphapqa = ", alphapqa); P = P * (1-1/p)^(alphapqa); ); \\ Computation of the sum over prime and characters S=0; for (j=2, phiq, \\ print("******************** j= ",j); chi=chars[j,]; \\ print("chi = ",chi); S= S + Schivector[j-1] * conj(chi[a]); ); \\ print("S = ", S); print("****** FINAL RESULT ********"); write("~/sessionsRIC.tex","****** FINAL RESULT ********"); write("~/sessionsRIC.tex", " "); C= P^(1/phiq) * exp(-S/phiq - Euler/phiq); Cvector[i]=real(C); print("C(q,a) with q = ", q, " and a = ",a, " is = ", real(C)); write("~/sessionsRIC.tex","C(q,a) with q = ", q, " and a = ",a, " is = ", real(C)); write("~/sessionsRIC.tex", " "); \\ FINAL ERRORS finalerr=abs(C)*abs(totalerr); finalerrvector[i]=finalerr; print("with error <= of ", finalerr); write("~/sessionsRIC.tex","with error <= of ", finalerr); write("~/sessionsRIC.tex", " "); correctdigits=floor(abs(log(finalerr)/log(10)))-2; if (finalerr>=1, write("~/sessionsRIC.tex","err4 is LARGE !!!!, change N and/or T"); write("~/sessionsRIC.tex", " "); print("err4 is LARGE !!!!, change N and/or T"), write("~/sessionsRIC.tex","number of the correct decimal digits = ", correctdigits); write("~/sessionsRIC.tex", " "); print("number of the correct decimal digits = ", correctdigits)); Cresults[i,]=[q,a,real(C),correctdigits]; TruncCresults[i,]=[q,a,truncate(real(C)*10^40)/10^40.,correctdigits]; i=i+1; ); if ((q==3 && a==1), write("~/sessionsRIC.tex", "error from Finch's computation of C(3,1) = ", 1.4034774468-real(C)); write("~/sessionsRIC.tex", " "); print("error from Finch's computation of C(3,1) = ", 1.4034774468-real(C))); if ((q==3 && a==2), write("~/sessionsRIC.tex", "error from Finch's computation of C(3,2) = ", 0.6000732161-real(C)); write("~/sessionsRIC.tex", " "); print("error from Finch's computation of C(3,2) = ", 0.6000732161-real(C))); if ((q==4 && a==1), write("~/sessionsRIC.tex", "error from Finch's computation of C(4,1) = ", 1.2923041571-real(C)); write("~/sessionsRIC.tex", " "); print("error from Finch's computation of C(4,1) = ", 1.2923041571-real(C))); if ((q==4 && a==3), write("~/sessionsRIC.tex", "error from Finch's computation of C(4,3) = ", 0.8689277682-real(C)); write("~/sessionsRIC.tex", " "); print("error from Finch's computation of C(4,3) = ", 0.8689277682-real(C))); ); \\ product verification print("****** PRODUCT VERIFICATION ********"); print("VERIFICATION OF EQ. (24) OF THE PAPER"); write("~/sessionsRIC.tex","****** PRODUCT VERIFICATION ********"); write("~/sessionsRIC.tex", " "); write("~/sessionsRIC.tex","VERIFICATION OF EQ. (24) OF THE PAPER"); write("~/sessionsRIC.tex", " "); prodC=1; \\ print(Cvector, finalerrvector); for (u=1, phiq, prodC = prodC * Cvector[u]); deltaprodC = prodC - exp(-Euler)*q/phiq; \\OUTPUTS on three files: results= complete results in TeX format; \\ resultsshort = results truncated to 40 digits in TeX format; \\ rawresults: complete results in Pari Format; needed for subsequent checks write("~/resultsRIC.tex","These are the constants for q = ", q); write("~/resultsRIC.tex", " "); write("~/resultsRIC.tex","The first column is the modulus q"); write("~/resultsRIC.tex", " "); write("~/resultsRIC.tex","The second column is the class a"); write("~/resultsRIC.tex", " "); write("~/resultsRIC.tex","The third column is the constant C(q,a)"); write("~/resultsRIC.tex", " "); write("~/resultsRIC.tex","The fourth column is number of correct digits in C(q,a)"); write("~/resultsRIC.tex", " "); write("~/resultsRIC.tex","$$"); writetex("~/resultsRIC.tex",Cresults); write("~/resultsRIC.tex","$$"); write("~/resultsRIC.tex", " "); write("~/resultsRIC.tex","The product of these constants is = ", prodC); write("~/resultsRIC.tex", " "); write("~/resultsRIC.tex","The error from its expected value exp(-gamma)q/phi(q) is = ", deltaprodC); write("~/resultsRIC.tex", " "); write("~/rawresultsRIC.tex",Cresults); write("~/resultsshortRIC.tex","These are the constants for q = ", q); write("~/resultsshortRIC.tex", " "); write("~/resultsshortRIC.tex","The first column is the modulus q"); write("~/resultsshortRIC.tex", " "); write("~/resultsshortRIC.tex","The second column is the class a"); write("~/resultsshortRIC.tex", " "); write("~/resultsshortRIC.tex","The third column is the constant C(q,a) (we are just writing 40 decimal digits)"); write("~/resultsshortRIC.tex", " "); write("~/resultsshortRIC.tex","The fourth column is number of correct digits in C(q,a)"); write("~/resultsshortRIC.tex", " "); write("~/resultsshortRIC.tex","$$"); default(realprecision,42); writetex("~/resultsshortRIC.tex",TruncCresults); write("~/resultsshortRIC.tex","$$"); write("~/resultsshortRIC.tex", " "); TruncprodC = truncate(prodC*10^40)/10^40.; write("~/resultsshortRIC.tex","The product of these constants is = ", TruncprodC); write("~/resultsshortRIC.tex", " "); default(realprecision, defaultprecision); print("The product of these constants is = ", prodC); print("The error from its expected value exp(-gamma)q/phi(q) is = ", deltaprodC); write("~/sessionsRIC.tex","The error from its expected value exp(-gamma)q/phi(q) is = ", deltaprodC); write("~/sessionsRIC.tex", " "); write("~/resultsshortRIC.tex","The error from its expected value exp(-gamma)q/phi(q) is = ", deltaprodC); write("~/resultsshortRIC.tex", " "); elaptimefinalcomp=gettime(); print("****** ELAPSED TIME ********"); write("~/sessionsRIC.tex", "****** ELAPSED TIME ********"); write("~/sessionsRIC.tex", " "); write("~/resultsshortRIC.tex", "****** ELAPSED TIME ********"); write("~/resultsshortRIC.tex", " "); write("~/resultsRIC.tex", "****** ELAPSED TIME ********"); write("~/resultsRIC.tex", " "); seconds=floor(elaptimeprecomp/1000)%60; minutes=floor(elaptimeprecomp/60000); millisec=elaptimeprecomp- minutes*60000 - seconds*1000; write("~/resultsshortRIC.tex", "Elapsed time for the precomputation: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); write("~/resultsshortRIC.tex", " "); write("~/resultsRIC.tex", "Elapsed time for the precomputation: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); write("~/resultsRIC.tex", " "); print("Elapsed time for the precomputation: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); write("~/sessionsRIC.tex", "Elapsed time for the precomputation: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); write("~/sessionsRIC.tex", " "); seconds=floor(elaptimefinalcomp/1000)%60; minutes=floor(elaptimefinalcomp/60000); millisec=elaptimefinalcomp- minutes*60000 - seconds*1000; write("~/resultsshortRIC.tex","Elapsed time for the final computation: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); write("~/resultsshortRIC.tex", " "); write("~/resultsRIC.tex","Elapsed time for the final computation: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); write("~/resultsRIC.tex", " "); print("Elapsed time for the final computation: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); write("~/sessionsRIC.tex", "Elapsed time for the final computation: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); write("~/sessionsRIC.tex", " "); seconds=floor((elaptimefinalcomp+elaptimeprecomp)/1000)%60; minutes=floor((elaptimefinalcomp+elaptimeprecomp)/60000); millisec=elaptimefinalcomp+elaptimeprecomp- minutes*60000 - seconds*1000; write("~/resultsshortRIC.tex","Total elapsed time: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); write("~/resultsshortRIC.tex", " "); write("~/resultsRIC.tex","Total elapsed time: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); write("~/resultsRIC.tex", " "); print("Total elapsed time: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); write("~/sessionsRIC.tex","Total elapsed time: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); write("~/sessionsRIC.tex", " "); write("~/timeRIC.tex", q, " & ", minutes, " & ", seconds, " & ", millisec, " \cr"); write("~/timeRIC.tex", " "); print("****** END PROGRAM ********"); write("~/sessionsRIC.tex","****** END PROGRAM ********"); write("~/sessionsRIC.tex", " "); write("~/resultsshortRIC.tex","****** END PROGRAM ********"); write("~/resultsshortRIC.tex", " "); write("~/resultsRIC.tex","****** END PROGRAM ********"); write("~/resultsRIC.tex", " "); } /*** Computation of the n-Bernoulli-chi number *** *** chistar is a primitive character modulo f; f is the conductor of chistar ***/ {chibern(n,chistar,f)= local(S, r); \\ print("start chibern"); \\just for testing S=0; for(r=1,f-1, S = S + chistar[r]*bernpol(n,r/f)); S= S*f^(n-1); return(S); } /*** n-Bernoulli polynomial computed in x ***/ {bernpol(n,x)= local(S, j); \\ print("start berpol");\\just for testing S=0; for(j=0,n, S = S + binomial(n,j)*bernvector[j+1]*x^(n-j)); \\ print("fine bernpol"); return(S); } /*** Computation of Dirichlet L(chi,s); Re(s)>1, using the chi-formula of Euler-McLaurin *** *** chi is a primitive character mod q; *** *** truncation level N = lq (positive integer) of the finite sum for L(s,chi) *** *** truncation level T (even positive integer > 3) of the finite Euler-Mclaurin sum for L(s,chi) ***/ {LEuler(s, chi, q, T, N)= local(S, S1, prods, j, j1, i, LTN, V, conf); \\ print("start L(s,chi,q, T, N)"); \\just for testing S=0; S1=0; prods = s; for (j=1, N-1, j1= j%q; if (j1==0, j1=q); S1 = S1 + chi[j1]/(j^s)); for (i=3, T, prods = prods * (s+i-2); S = S + (-1)^(i-1) * chibern(i,chi,q) * prods * N^(-i+1)/ factvector[i] ); V= N^(-s); S= V * S; \\ \\ print("s=",s); \\ B_0(chi)=0 because chi is a Dirichlet character \\ LTN = - chibern(1,chi,q) * V + chibern(2,chi,q) * s * V * N^(-1)/ 2 - S + S1; \\ \\ print("termine ",chibern(1,chi,q) * V ); \\ print("LTN ",LTN); \\ needed to evaluate the final errors \\ conf=abs(LTN); if ( conf < U, U = conf); \\ \\ print("LTN = ", LTN); \\ print("fine L(s, chi, q, T, N)"); \\ return(LTN) } /*** Computation of the rootnumber(chi); for q<=100 the Gauss sum is used *****/ {rootnumber(chistar, f, e) = local(S, j, j1, W, e1); S=0; e1=0; if (e==0, e1=1, e1=-I); \\ e=0,1 depending on the parity of chi \\ \\ if chistar is quadratic then W(chistar)=1 \\ if (chistar[f+1]==2, W=1; \\ print ("chistar is a quadratic character") ,if (f<=100, \\ \\ print ("chistar is not a quadratic character"); \\ print("chistar = ", chistar); \\ computation of the GAUSS sum \\ remark: chistar[f]=0 \\ for (j=1, f-1, S = S+ chistar[j] * exp(2*Pi*I*j/f)); W = e1 * S/sqrt(f), \\ \\ remark: in this case the sum is an infinite series; I compute just 100000 terms (is it enough ???) \\ for (j=1, 10000, j1= j%f; if (j1==0, j1=f); S = S+ chistar[j1] * exp(-Pi*j^2)); W = e1 * S/conj(S) ); ); \\ \\ print("f = ", f ," e = ", e, " e1 = ", e1); \\ print("Gauss sum = ", S ," module = ", abs(S)); \\ print("rootnumber = ", W ," module = ", abs(W)); return(W) } /*** Computation of the tail of the Euler product for L(s,chi); Re(s)>1, *** **** B1 is the limit of the primes used ***/ {Lprodtail(s, B1, chi, q, T, N, princ) = local(P, P1, P2, p1, q1, u, L, W, e, valchi, Ltail, prim, cond, order, f, chistar, phif, checkchar, induces,a1); \\ print("start Lprodtail"); \\ print("s=", s); L=1; P=1; P1=1; P2=1; u=-s; prim=chi[q+2]; cond=chi[q+3]; order=chi[q+1]; \\ \\ computation of the finite Euler product for L(s,chi) with primes less or equal than B1=B \\ forprime(p=2, B1 , p1= p%q; if (p1==0, p1=q); P=P*(1-chi[p1]*p^u)); \\check if chi is the principal character \\ if ((princ==1 && s<>1), \\ print ("** principal character and s<>1; I can use the Riemann zeta"); \\ \\ computation of the finite product needed to pass from zeta(s) to L(s,chi_0), chi_0 modulo q forprime(p=2, q , q1= q%p; if (q1==0, P1=P1*(1-p^u))); \\ \\ computation of L(s,chi_0) \\ L = zetavector[s] * P1; \\ \\ computation of the tail of L(s,chi_0) with primes > B Ltail= L * P; \\ \\ print ("Ltail= ", Ltail); return(Ltail) ); \\ chi is not the principal character \\ \\ chi is a primitive character ?? \\ if (prim == 0, \\ computation of chistar that induces the imprimitive character chi \\ print("imprimitive character"); f=cond; phif=eulerphi(f); \\ print("f= ",f, " cond= ", cond, " phif= ", phif); chistar=vector(f+3); checkchar = matrix(phif,f+3); checkchar = genchar(f); \\ print(checkchar); for (k=1,phif, induces=1; \\ print("k= ", k); for(a=2, q, if(gcd(a,q)==1, \\ print("a= ", a); a1=a%f; if (a1==0, a1=f); if (abs(chi[a]-checkchar[k,a1])>=defaulterror, induces=0*induces, induces=1*induces) ) ); if (induces ==1, chistar=checkchar[k,]; \\ printtex("chistar = ", chistar); \\ print("induces = ", induces); \\ printtex("chi = ", chi); \\ print("checkchar = ", checkchar[k,]); break; ); ); \\ \\ computation of the finite product for the imprimitive character \\ forprime (p=2, q, if(q%p==0, p1=p%f; if (p1==0, p1=f); P2 = P2 * (1-chistar[p1]*p^u)); ) , \\ print("primitive character"); f=q; chistar=vector(q+3); chistar=chi; \\ print("chi = chistar = ", chistar); ); \\ computation of the parity of chistar; e=0 if chistar(-1)=1; e=1 if chistar(-1)=-1 \\ valchi = chistar[f-1]; \\print("valchi=", valchi); \\print("defaulterror=", defaulterror); if (abs(valchi-1)9600, B=9600); Smtail=0; chik=vector(q+3); \\ print("chi = ", chi); for(k=1,K, for(j=1,q, chik[j]=(chi[j])^k); if ((k%chi[q+1]) == 0, princ=1, princ=0; if(isqprime, chik[q+2]=1; chik[q+3]=q; chik[q+1]= chi[q+1]/gcd(chi[q+1],k), \\ print("primitivity test on chik and computation of its conductor"); \\ ndiv = numdiv(q); \\print("ndiv", ndiv); div = vector(ndiv); div = divisors(q); \\print("div", div); testprim=1; cond=q; \\print("cond", cond); \\print("k = ", k); for (i = 2, ndiv-1, flag=0; \\print("i=",i); forstep (a=div[i]+1, q, div[i], \\print("a=",a); if (gcd(a,q)==1 && abs(chik[a]-1)>= defaulterror, flag=1; break); ); \\ print("flag = ", flag); if (flag == 0, testprim=0; cond = div[i]; break); flag=0; ); chik[q+2]=testprim; chik[q+3]=cond; chik[q+1]= chi[q+1]/gcd(chi[q+1],k); ); ); \\ print("m= ", m," k= ", k, " chik=", chik); \\ print("princ =", princ); \\ if (princ==0, print("chi^k is NOT principal")); \\ if (princ==1, print("chi^k is principal")); Smtail = Smtail + (moebius(k) * log(Lprodtail(k*m, B, chik, q, T, N, princ))/k) ); \\ print("end Smchi"); return(Smtail) } /***************** DIRICHLET CHARACTERS GENERATION ****************/ /********** GENERAL CASE *******************/ /*** Dirichlet characters modulo an integer q such that Zq^* is not cyclic ******/ {genchar(q)= local(phiq, u, l, i, k, vec, numprimefactors, car1,car2, order1,order2, prim1,prim2, cond,cond1, cond2, carresult, factorization,q1,q2,numcar1,numcar2,l1,l2,ndiv,div,flag,testprim, charactersgen); charactersgen=matrix(phiq,q+3,lines,rows, 0); \\ \\ we have a row for every character \\ and a column for every class in Z_q \\ the q+1 column contains the order of the character in \hat{Z_q*} \\ the q+2 column is 1 if the character is primitive; 0 otherwise \\ the q+3 column contains the conductor f of the character (=q if primitive) \\ that induces this one \\\\ q is prime? if (isprime(q)==1, return(gencharmodprime(q))); \\Zq^* is cyclic ? vec = vector(3); vec = znstar(q); if (length(vec[2])==1, return(gencharcyclic(q))); \\ print("Z/qZ^* with q= ",q," is not cyclic"); \\ q is a power of two >=8? (the case q=2,4 are detected with isprime and cyclic) \\ factorization is a matrix containing the complete factorization of q factorization = factor(q); \\ calcolo e controllo del numero di fattori primi distinti di q numprimefactors= length(factorization[,1]); if ((numprimefactors==1 && factorization[1,1]==2), \\ print("q= ",q," is a power of two >=8"); return(gencharpowersoftwo(q,factorization[1,2]))); \\ GENERAL CASE: Zq^* is not cyclic and q is not power of two >=8 phiq=vec[1]; \\ print("phiq= ",phiq); car1=matrix(phiq,q+3,lines,rows, 0); car2=matrix(phiq,q+3,lines,rows, 0); carresult=matrix(phiq,q+3,lines,rows, 0); \\ We build a character mod q starting from characters defined mod prime powers dividing q \\ the case of powers of two has to be considered in a special way \\\\\\\\\\\\\\\\\\\\\\\ order1= vector(phiq); order2= vector(phiq); prim1= vector(phiq); prim2= vector(phiq); cond1= vector(phiq); cond2= vector(phiq); if (factorization[1,1]==2, if (factorization[1,2]==1, q1=2; car1=gencharmodprime(2); \\print(car1); \\car1=characters , q1=(2^factorization[1,2]) ; car1=gencharpowersoftwo(q1,factorization[1,2]); \\print(car1); \\car1=characters ), q1=factorization[1,1]^factorization[1,2]; car1=gencharcyclic(q1); \\car1=characters; \\print(car1); ); numcar1= length(car1[,1]); order1=car1[,q1+1]; prim1=car1[,q1+2]; cond1=car1[,q1+3]; \\ print("car1",car1); for (k=2, numprimefactors, q2=factorization[k,1]^factorization[k,2]; car2=gencharcyclic(q2); \\car2=characters; numcar2= length(car2[,1]); order2=car2[,q2+1]; prim2=car2[,q2+2]; cond2=car2[,q2+3]; for (i=1, numcar1, for (j=1, numcar2, for(l=1,q, l1= l%q1; if (l1==0, l1=q1); l2= l%q2; if (l2==0, l2=q2); carresult[(i-1)*numcar2+j,l] = car1[i,l1] * car2[j,l2] ); carresult[(i-1)*numcar2+j,q+1]=lcm(order1[i],order2[j]); \\ primitivity test on the generated character and computation of the conductor testprim=1; cond=q; ndiv = numdiv(q); div = vector(ndiv); div = divisors(q); for (u=2, ndiv-1, flag=0; forstep (a=div[u]+1, q, div[u], if (gcd(a,q)==1 && (round(real(carresult[(i-1)*numcar2+j,a]))) <> 1, flag=1; break); ); \\ print("flag = ", flag); if (flag == 0, testprim=0; cond = div[u]; break); flag=0; ); carresult[(i-1)*numcar2+j,q+2]=testprim; carresult[(i-1)*numcar2+j,q+3]=cond; ) ); car1=carresult; q1=q1*q2; numcar1= numcar1*numcar2; order1=car1[,q+1]; prim1=car1[,q+2]; cond1=car1[,q+3]; ); charactersgen=carresult; \\ print("the characters modulo q = " ,q," are "); \\ printtex(charactersgen); \\ printp(charactersgen); return(charactersgen) } /*** Dirichlet characters modulo a prime p ***/ /*** this is the easiest case ***/ {gencharmodprime(p)= local(j, k, u, v, g, g1,charactersp); if (isprime(p)<>1, print("p= ",p," is not a prime"); return); \\ print(p," is a prime"); \\ \\ we have a row for every character \\ and a column for every class in Z_p \\ the p+1 column contains the order of the character in \hat{Z_p*} \\ the p+2 column is 1 if the character is primitive; 0 otherwise \\ the p+3 column contains the conductor f of the character (=q if primitive) \\ that induces this one charactersp=matrix(p-1,p+3); g1=znprimroot(p); g=lift(g1); u=2*Pi*I/(p-1); \\ print("u= ", u); \\ print("g= ",g); for (k=1,p-1, charactersp[1,k]=1); charactersp[1,p]=0; charactersp[1,p+1]=1; charactersp[1,p+2]=0; \\ the principal character mod p, p odd, is not primitive if (p==2,charactersp[1,p+2]=1); \\ the principal character mod 2, is primitive charactersp[1,p+3]=p; \\ every character mod a prime has f=p for(j=2,p-1, charactersp[j,1]=1; charactersp[j,p]=0; charactersp[j,p+2]=1; \\ every non principal character mod a prime is primitive charactersp[j,p+3]=p; \\ every character mod a prime has f=p \\ v=exp(u*(j-1)); \\ \\ order of the j-th character (it is identified by the (j-1)th power \\ of the prmitive root g \\ charactersp[j,p+1]= (p-1)/gcd(j-1,p-1); \\ \\ print("j= ", j); \\ print("v= ", v); \\ print("order of the character= ", charactersp[j,p+1]); \\ for (k=1,p-1, \\ print("k= ", k); print("g^k%p = ", g^k%p); print("v^k = ", v^k); charactersp[j,g^k%p]=v^k ); ); \\ print("the characters mod q = " ,p," are "); \\ printtex(charactersp); \\ printp(charactersp); return(charactersp) } /*** Dirichlet characters modulo an integer q such that Zq^* is cyclic ***/ /*** similar to the previous case, we have a primitive root ***/ {gencharcyclic(q)= local(phiq, u, v, l, g, g1, g2, k, div, ndiv, a, testprim, cond, flag, characterscycl); if (isprime(q)==1, return(gencharmodprime(q))); \\ print("Z/qZ^* with q= ",q," is cyclic"); \\ \\ we have a row for every character \\ and a column for every class in Z_q \\ the q+1 column contains the order of the character in \hat{Z_q*} \\ the q+2 column is 1 if the character is primitive; 0 otherwise \\ the q+3 column contains the conductor f of the character (=q if primitive) \\ that induces this one g1=znprimroot(q); g=lift(g1); phiq=eulerphi(q); \\ print("phiq= ",phiq); characterscycl=matrix(phiq,q+3); u=2*Pi*I/phiq; \\ print("u= ", u); \\ print("g= ",g); ndiv = numdiv(q); div = vector(ndiv); div = divisors(q); for(l=1,phiq, characterscycl[l,1]=1; characterscycl[l,q]=0; v=exp(u*(l-1)); \\ \\ order of the l-th character (it is identified by the (l-1)th power \\ of the primitive root g) \\ remark: if l=1, pari gives gcd(0,phiq)=phiq \\ characterscycl[l,q+1]= phiq/gcd(l-1,phiq); \\ \\ print("l= ", l); \\ print("v= ", v); \\ print("order of the character= ", characterscycl[l,q+1]); \\ for (k=1,q-1, \\ print("k= ", k); print("g^k%q = ", g^k%q); print("v^k = ", v^k); g2=g^k; if (gcd(g2,q)==1, characterscycl[l,g2%q]=v^k, characterscycl[l,g2%q]=0); ); \\ \\ primitivity test on the l-th character and computation of the conductor \\ testprim=1; cond=q; for (i=2, ndiv-1, flag=0; forstep (a=div[i]+1, q, div[i], if (gcd(a,q)==1 && (round(real(characterscycl[l,a]))) <> 1, flag=1; break); ); \\ print("flag = ", flag); if (flag == 0, testprim=0; cond = div[i]; break); flag=0; ); characterscycl[l,q+2]=testprim; characterscycl[l,q+3]=cond; ); \\ print("the characters mod q = " ,q," are "); \\ printtex(characterscycl); \\ printp(characterscycl); return(characterscycl) } /*** Dirichlet characters modulo an integer q =2^a with a>=3 ***/ {gencharpowersoftwo(q,power)= local(phiq, u,u2, v, l, g, g2, k, div, ndiv, a, testprim, cond, flag, veceulerphi, tabindex, tabchar, mult, elem, aux,vecelem,j1,j2, exponent, index,car,P, characterstwo); \\ we have a row for every character \\ and a column for every class in Z_q \\ the q+1 column contains the order of the character in \hat{Z_q*} \\ the q+2 column is 1 if the character is primitive; 0 otherwise \\ the q+3 column contains the conductor f of the character (=q if primitive) \\ that induces this one \\ we generate the characters mod a power of two as described in Davenport's book phiq=eulerphi(q); \\ print("phiq= ",phiq); characterstwo=matrix(phiq,q+3,lines,rows, 0); /*** "primitive roots" for powers of two >=8 *****/ g = vector(2); g2 = vector(2); u = vector(2); u2 = vector(2); v = vector(2); veceulerphi = vector(2); index = vector(2); tabindex = matrix(phiq,2); tabchar = matrix(phiq,2); exponent = power; \\ this is an integer g[1]=-1; u[1]=Pi*I; veceulerphi[1]=2; g[2]=5; veceulerphi[2]=2^(exponent-2); u[2]=2*Pi*I/veceulerphi[2]; \\ print("exponent = ", exponent); mult=vector(2); mult[1]=1; mult[2]=veceulerphi[2]; \\ print("mult = ", mult); for (j1=1, veceulerphi[1], for (j2=1, veceulerphi[2], l=(j1-1)*mult[2]+j2; tabindex[l,1]=j1-1; \\ print("riga = ", l); print("j1= ", j1); tabindex[l,2]=j2-1; \\ print("j2= ", j2); )); /*** This part is to generate the table of the indexes for the characters ***/ tabchar=tabindex; \\ print("tabchar =", tabchar); \\ print("tabindex =", tabindex); /*** vecelem collects the generated elements ***/ vecelem = vector(phiq); /*** Used to test the primitivity *****/ ndiv = numdiv(q); div = vector(ndiv); div = divisors(q); /*** Generation of the characters ***/ /*** k counts the characters ***/ /*** l counts the elements ***/ /*** the elements are generated using "primitive roots" ***/ for(k=1,phiq, characterstwo[k,1]=1; \\ characterstwo[k,q]=0; (character is initialized with zeros) car=tabchar[k,]; \\ print("car = ", car); v[1]=exp(u[1]*car[1]); v[2]=exp(u[2]*car[2]); \\ \\ order of the k-th character \\ remark: pari gives gcd(0,phiq)=phiq \\ characterstwo[k,q+1]= lcm(veceulerphi[1]/gcd(car[1],veceulerphi[1]), veceulerphi[2]/gcd(car[2],veceulerphi[2])); \\ \\ print("l= ", l); \\ print("v= ", v); \\ print("order of the character= ", characterstwo[k,q+1]); for(l=1,phiq, index=tabindex[l,]; \\ print("index = ",index); P=2^exponent; elem = Mod(g[1]^index[1] * g[2]^index[2], P); \\ print("elem = ", elem); vecelem[l]=lift(elem); elem=vecelem[l]; aux=1; for (i=1, 2, aux=aux*v[i]^index[i]); characterstwo[k,elem]= aux; \\ print("index= ", index); print("g^index%q = ", elem); \\ print("v^index = ", aux); ); \\ primitivity test on the k-th character and computation of the conductor testprim=1; cond=q; for (i=2, ndiv-1, flag=0; forstep (a=div[i]+1, q, div[i], if (gcd(a,q)==1 && (round(real(characterstwo[k,a]))) <> 1, flag=1; break); ); \\ print("flag = ", flag); if (flag == 0, testprim=0; cond = div[i]; break); flag=0; ); characterstwo[k,q+2]=testprim; characterstwo[k,q+3]=cond; ); \\ print("vecelem = ", vecelem); \\ print("the characters mod q = " ,q," are "); \\ printtex(characterstwo); \\ printp(characterstwo); return(characterstwo) } /*** Function to evaluate err4 for optimizing T ***/ {valutT(q,N,K,M,L)=local(prev,P,T,j,finalerror); default(realprecision,100); \\ set working precision; q=4; N=840*100; K=25; M=20; L=200; bernvector=vector(L+1); factvector=vector(L+1); for(j=1,L+1,bernvector[j]=bernfrac(j-1); factvector[j]=j!); finalerror=0; for(s=1,L, P=1; prev=10^100; print("s= ", s); forstep(T=2,L,2, for(j=0,T-2,P=P*(s+j)); print("T= ", T); finalerror= (2 * (M * K + T-1)^(T-2) * abs(bernvector[T+1]) * q^T)/(factvector[T] * N^(T-1) *(N-1)); \\ print(finalerror);print(prev); if(finalerror>prev,print("Optimal T is ",T," for s= ",s," with error: ",prev); return); prev=finalerror; ); ); print("Optimal error for the last T: ",finalerror); } /*** Function to evaluate err4 for optimizing T ***/ {valT(q)=local(N,K,M,T,finalerror,err1,err2,A); default(realprecision,110); \\ set working precision; \\ IMPORTANT: N has to be a multiple of q; T has to be even and positive A=3200; for(q=41, 50, if(q>=3 && q<=10, M = 26; K = 26; N = ((840*20\q)+1)*q; T = 88); \\ for the 3-10 range if(q>=11 && q<=20, M = 26; K = 26; N = ((840*25\q)+1)*q; T = 106); \\ for the 11-20 range if(q>=21 && q<=30, M = 26; K = 26; N = ((840*32\q)+1)*q; T = 116); \\ for the 21-30 range if(q>=31 && q<=40, M = 26; K = 26; N = ((840*36\q)+1)*q; T = 128); \\ for the 31-40 range if(q>=41 && q<=50, M = 26; K = 26; N = ((840*40\q)+1)*q; T = 138); \\ for the 41-50 range if(q>=51 && q<=60, M = 26; K = 26; N = ((840*44\q)+1)*q; T = 146); \\ for the 51-60 range if(q>=61 && q<=70, M = 26; K = 26; N = ((840*66\q)+1)*q; T = 128); \\ for the 61-70 range if(q>=71 && q<=80, M = 26; K = 26; N = ((840*70\q)+1)*q; T = 134); \\ for the 71-80 range if(q>=81 && q<=90, M = 26; K = 26; N = ((840*70\q)+1)*q; T = 144); \\ for the 81-90 range if(q>=91 && q<=100, M = 26; K = 26; N = ((840*73\q)+1)*q; T = 152); \\ for the 91-100 range /***** OLD if(q>=61 && q<=70, M = 26; K = 26; N = ((840*48\q)+1)*q; T = 154); \\ for the 61-70 range if(q>=71 && q<=80, M = 26; K = 26; N = ((840*40\q)+1)*q; T = 206); \\ for the 71-80 range if(q>=81 && q<=90, M = 26; K = 26; N = ((840*45\q)+1)*q; T = 206); \\ for the 81-90 range \\ if(q>=91 && q<=100, M = 26; K = 26; N = ((840*49\q)+1)*q; T = 210); \\ for the 91-100 range if(q>=91 && q<=100, M = 26; K = 26; N = ((840*72\q)+1)*q; T = 152); \\ for the 91-100 range ********/ B=A*q; if (B>9600, B=9600); \\ err1 using the estimate |L_Aq|< (Aq)^(-k*m) * min (Aq/(k*m-1); err1= (2*B)/(K^2*(B^K-1)*(B-1)); print("err1= ", err1*1.0); err2= B/(B^M*(B-1)*M*(M-1)); print("err2= ", err2*1.0); \\print ("M = ", M, " K = ", K, " N = ", N, " T = ", T); finalerror= (2 * (M * K + T-1)^(T-2) * abs(bernfrac(T)) * q^T)/(T! * N^(T-1) *(N-1)); print("Error: ",finalerror*1.0); print("q = ", q ," num digits: ",floor(abs(log(finalerror)/log(10)))); ); } /*** Function to see the number of times we compute zeta(km) ***/ {times()=local(K,M,maxim1,col,err2,A,numbers); M = 26; K = 26; maxim1=max(M,K); numbers=vector(M*K,col,0); for(j=2,maxim1,numbers[j]=numbers[j]+1); for(m=2, M, for (k=2,K,numbers[m*k]=numbers[m*k]+1)); print(numbers); return() }