/* Copyright (C) 2024 Alessandro Languasco */ /********** ****** COMPUTATION OF sum_{p = a mod d} log p/ (p^s(p-1)); Re(s)> 0; in AP ****** FOR Euler's constants in AP (L.-Moree) **********/ \\ Global variables global(CHI_VECTOR, NUM_PRIME_FACTORS_D,P_PRIMES, VALCHI_CONJ_A_MATRIX, VALCHI_P_MATRIX); global(DFACTORS, TWOPII, S5VECTOR, S2VECTOR, GAMMA1VECTOR, LOG_DFACTORS, LOG_P, LOG_P_OVER_P_TO_S); global(SQUARE_FREE_UP_TO_J, LOG_DER_AT_1, SJCHI_VECTOR, NUMPRIMES_UP_TO_P, ONE_OVER_P_MINUS_ONE); global(ONE_OVER_P_D_MINUS_ONE, ORD_CHI_J_MATRIX, CHI_J_MATRIX, ISPRIMEDIVISOR_OF_D, NUMPRIMES_UP_TO_D, PRIMES_UP_TO_D); \\ Routine to compute the value of L'/L(s,chi) when chi is a Dirichlet character (also imprimitive) {L_derlog(L, chi, G, d, s) = local(conductor, Gstar, chistar, corr_chi_imprim, val_chistar_d1, ud1, L_value, L_prime_value, L_derlog_value, d1); L_value = lfun(L, s); \\ value of L at s for the attached primitive character L_prime_value = lfun(L, s, 1); \\ value of L'at s for the attached primitive character L_derlog_value = L_prime_value / L_value; \\ value of logarithmic derivatie at s for the attached primitive character conductor = zncharconductor(G, chi); if (conductor < d, \\ chi is imprimitive \\ computing the correction sum for the logarithmic derivative [Gstar,chistar] = znchartoprimitive(G, chi); \\ chistar is the primitive character such that chi is induced by chistar corr_chi_imprim = 0.; for(i=1, NUM_PRIME_FACTORS_D, d1 = DFACTORS[i]; \\print("d1 = ", d1); if(conductor % d1 <> 0, \\ d1 not divides both d and conductor of chi and conductor <> 1 ud1 = chareval(Gstar, chistar, d1); val_chistar_d1 = exp(TWOPII * ud1); \\ value of chistar(d1) corr_chi_imprim += val_chistar_d1 * LOG_DFACTORS[i] / (d1^s - val_chistar_d1) \\corr_chi_imprim += val_chistar_d1 * log(d1) / (d1^s - val_chistar_d1 ) ; ); ); L_derlog_value += corr_chi_imprim ; ); return(L_derlog_value); } \\ Routine to compute the sum p <= P of chi(p) log(p) / (p^s - chi(p))); to get the truncated log derivative L(s, chi) {chi_sum_fin(P, chi, G, d, s) = local(fin_sum_chi, val_chi_p, up, p); fin_sum_chi = 0.; \\ computation of the truncated sum over chi(p) log(p) / (p^s - chi(p) \\ working on (p,d) == 1 \\i = 1 ; \\forprime (p = 2, d, for (i = 1, NUMPRIMES_UP_TO_D, \\if ( d % p <> 0, \\ p does not divide d if ( ISPRIMEDIVISOR_OF_D[i] == 0, p = PRIMES_UP_TO_D[i]; up = chareval(G, chi , p ); val_chi_p = exp(TWOPII * up); fin_sum_chi += val_chi_p * LOG_P[i] /(p^s - val_chi_p); ); \\i+=1; ); \\i = NUMPRIMES_UP_TO_D + 1 ; \\ working on (p,d) == 1 \\forprime (p = d+1, P, for (i = NUMPRIMES_UP_TO_D + 1, NUMPRIMES_UP_TO_P, p = P_PRIMES[i]; up = chareval(G, chi , p ); val_chi_p = exp(TWOPII * up); fin_sum_chi += val_chi_p * LOG_P[i] /(p^s- val_chi_p); ); return(fin_sum_chi); } {lazygamma_AP(d, s, prec=19) = local(dminusone, A, K, P, err1, G, moebius_vector, Smallprimes_sum, chi, dcoprimes, toterr, frac_part_sigma, int_part_sigma, t, real_s, resfile, num_square_free, L2, elaptimecomp, seconds, minutes, millisec, ua, Sjchi, Sjkchi, err2, factorp, i, accuracy, phid, J, ktimesj, aux, sum_chars, numdcoprimes, p, chi_to_j, j, L2_derlog_value, L2trunc_derlog, fin_sum_chi_j, results); if (d > 100, error("d must be <= 100")); if (d <= 2, error("d must be >= 3")); if (floor(prec)<> prec, prec = floor(prec); print("prec forced to be an integer"); ); if (prec < 5, error("prec must be an integer >= 5")); default(realprecision, prec+10); \\ set working precision; dminusone = d-1; phid = eulerphi(d); resfile = fileopen("lazygamma_AP-results.csv", "a"); \\filewrite(resfile,"d; a; s ; lazygamma_das; accuracy"); print("****** COMPUTATION OF sum_{p = a mod d} log p/ (p^s(p-1)); Re(s)> 0; in AP *******"); print("****** WITH A PRECISION OF AT LEAST ", prec," DECIMAL DIGITS *******"); \\ settings of the main parameters for a precision of about 100 decimal digits gettime(); real_s = real(s); if (real_s <=0, error ("Re(s) must be positive")); t = imag(s); int_part_sigma = floor(real_s); frac_part_sigma = real_s - int_part_sigma; A = 3200; P = A*d; \\ is P in the file if (P > 9600, P = 9600); accuracy = 10^(-prec-10); K = 2; err1 = log(P) / ((K +frac_part_sigma - 1)*(P^(K-1+frac_part_sigma))*(P-1)); ; until (err1 < accuracy/2, K += 1; err1 = log(P) / ((K +frac_part_sigma - 1)*(P^(K-1+frac_part_sigma))*(P-1)); ); \\J = 26; J = 2; err2 = P^2/(P^((J+1)*int_part_sigma)*(P^(J+1)-1.0)); until (err2 < accuracy/2, J += 1; err2 = P^2/(P^((J+1)*int_part_sigma)*(P^(J+1)-1.0)); ); print("-------------------------- "); print ("Accuracy parameters: "); if (prec == 19 , print("default precision is with prec = 19")); print("prec = ", prec, "; A = ", A, "; P = ", P, "; K = ", K, "; J = ", J); print("Accuracy is 10^(-",prec+10,")"); print("err1 = ", err1); print("err2 = ", err2); toterr = err1+err2; print("toterr = ", toterr); print("-------------------------- "); \\ setting the group G = znstar(d, 1); TWOPII = 2*Pi*I; \\ PRECOMPUTATIONS moebius_vector = vector(J); moebius_vector[1] = 1; for (j = 2, J, moebius_vector[j] = moebius(j) * 1.0); num_square_free = hammingweight(moebius_vector); SQUARE_FREE_UP_TO_J = vector(num_square_free); SQUARE_FREE_UP_TO_J[1] = 1; i=2; for (j = 2, J, if (moebius_vector[j] <> 0, SQUARE_FREE_UP_TO_J[i] = j; i+=1 ); ); NUMPRIMES_UP_TO_P = primepi(P); LOG_P = vector(NUMPRIMES_UP_TO_P); LOG_P_OVER_P_TO_S = vector(NUMPRIMES_UP_TO_P); ONE_OVER_P_MINUS_ONE = vector(NUMPRIMES_UP_TO_P); P_PRIMES = vector(NUMPRIMES_UP_TO_P); i=1; forprime ( p=2, P, factorp = p; P_PRIMES[i] = factorp; LOG_P[i] = log(factorp); LOG_P_OVER_P_TO_S[i] = LOG_P[i]/factorp^s; ONE_OVER_P_MINUS_ONE[i] = 1.0/(factorp - 1.0); i += 1 ); \\ generating the list of integers <= d-1 coprime with d dcoprimes = vector(phid); i=1; dcoprimes[1]=1; for(m=2, dminusone, if (gcd(m,d) == 1, i+=1; dcoprimes[i]=m ) ); numdcoprimes = #dcoprimes; \\ generating the non-principal Dirichlet characters mod d \\ and the orders of his powers CHI_VECTOR = vector(phid); CHI_J_MATRIX = matrix(phid, J); i=1; foreach(dcoprimes, m, \\ m = 1 corresponds to the principal character mod d CHI_VECTOR[i]=znconreychar(G,m); CHI_J_MATRIX[i,1] = CHI_VECTOR[i]; \\for ( j = 2, J, foreach(SQUARE_FREE_UP_TO_J,j, CHI_J_MATRIX[i,j] = charpow(G, CHI_J_MATRIX[i,1], j); ); i+=1; ); \\ DFACTORS contains the list of prime factors of d \\ it is needed to handle the L functions associated to imprimitive characters DFACTORS=factor(d)[,1]; NUM_PRIME_FACTORS_D = #DFACTORS; \\print(DFACTORS); \\print(NUM_PRIME_FACTORS_D); LOG_DFACTORS = vector(NUM_PRIME_FACTORS_D); ONE_OVER_P_D_MINUS_ONE = vector(NUM_PRIME_FACTORS_D); for(i=1, NUM_PRIME_FACTORS_D, p = DFACTORS[i]; factorp = p*1.0; ONE_OVER_P_D_MINUS_ONE[i] = 1.0/(factorp - 1.0); LOG_DFACTORS[i] = log(factorp); ); \\ ISPRIMEDIVISOR_OF_D contains 1 or 0 depending if p<=d divides d or not NUMPRIMES_UP_TO_D = primepi(d); ISPRIMEDIVISOR_OF_D = vector(NUMPRIMES_UP_TO_D); PRIMES_UP_TO_D = vector(NUMPRIMES_UP_TO_D); i = 1 ; forprime (p = 2, d, PRIMES_UP_TO_D[i] = p; if ( d % p == 0, ISPRIMEDIVISOR_OF_D[i] = 1 ); i+=1 ); \\ S5VECTOR stores the value of the sum over each non trivial character S5VECTOR = vector(numdcoprimes); \\ initialised with zeros SJCHI_VECTOR = vector(phid); for(m=1, phid, \\ we run over the characters chi = CHI_VECTOR[m]; Sjchi = 0; \\for (j=1, J , foreach(SQUARE_FREE_UP_TO_J,j, chi_to_j = CHI_J_MATRIX[m,j]; L2 = lfuncreate([G, chi_to_j]); Sjkchi = 0; for(k = 1 , K, ktimesj = j*(k + s); fin_sum_chi_j = chi_sum_fin(P, chi_to_j, G, d, ktimesj); L2_derlog_value = L_derlog(L2, chi_to_j, G, d, ktimesj); L2trunc_derlog = L2_derlog_value + fin_sum_chi_j; Sjkchi += L2trunc_derlog; ); Sjchi += moebius_vector[j]* Sjkchi; ); SJCHI_VECTOR[m] = Sjchi; ); VALCHI_CONJ_A_MATRIX = matrix(numdcoprimes,phid); i=0; foreach(dcoprimes, a, i += 1 ; sum_chars = 0; for(m=1, phid, \\ we run over the characters chi = CHI_VECTOR[m]; ua = chareval(G, chi , a ); VALCHI_CONJ_A_MATRIX[i,m] = exp(-TWOPII * ua); \\ conjugate character sum_chars += VALCHI_CONJ_A_MATRIX[i,m]* SJCHI_VECTOR[m]; ); S5VECTOR[i] = -sum_chars/phid; ); S2VECTOR = vector(numdcoprimes); i=0; foreach(dcoprimes, a, i += 1 ; Smallprimes_sum = 0; j = 0; forprimestep(p = 2, P, Mod( a, d), j+=1; Smallprimes_sum += log(p)/( p^s*(p-1) ); ); S2VECTOR[i] = Smallprimes_sum; ); print(" RESULTS "); elaptimecomp=gettime(); print("-------------------------- "); i=0; results = matrix(numdcoprimes,2); foreach(dcoprimes, a, i += 1 ; aux = S5VECTOR[i] + S2VECTOR[i] ; if (t == 0, aux = real(aux)); results[i,1] = aux; results[i,2] = toterr; print("lazygamma_AP(",d,";",a,";",s,") = ", aux); filewrite1(resfile,d); filewrite1(resfile,";"); filewrite1(resfile,a); filewrite1(resfile,";"); filewrite1(resfile,s); filewrite1(resfile,";"); filewrite1(resfile,aux); filewrite1(resfile,";"); default(format,e); default(realprecision,5); filewrite(resfile,toterr); default(format,f); default(realprecision, prec+10); ); print("-------------------------- "); \\fileclose(resfile); print("-------------------------- "); print("****** COMPUTATION TIME ********"); seconds=floor(elaptimecomp/1000)%60; minutes=floor(elaptimecomp/60000); millisec=elaptimecomp - minutes*60000 - seconds*1000; print("Computation time: ", minutes, " min, ", seconds, " sec, ", millisec, " millisec"); print("****** END PROGRAM ********"); \\return(results); } /************ Last login: Sun Apr 21 17:16:48 on ttys009 languasc@languasc-pro programs % gp2c-run -pmy_ -g -W lazygamma_AP.gp Warning:lazygamma_AP.gp:318: variable undeclared e Warning:lazygamma_AP.gp:321: variable undeclared f Reading GPRC: /Users/languasc/.gprc GPRC Done. GP/PARI CALCULATOR Version 2.15.5 (released) arm64 running darwin (aarch64/GMP-6.3.0 kernel) 64-bit version compiled: Mar 2 2024, Apple clang version 15.0.0 (clang-1500.1.0.2.5) threading engine: pthread (readline v8.1 enabled, extended help enabled) Copyright (C) 2000-2022 The PARI Group PARI/GP is free software, covered by the GNU General Public License, and comes WITHOUT ANY WARRANTY WHATSOEVER. Type ? for help, \q to quit. Type ?18 for how to get moral (and possibly technical) support. parisizemax = 2048000000, primelimit = 500000, nbthreads = 8 ? lazygamma_AP(3,1) ****** COMPUTATION OF sum_{p = a mod d} log p/ (p^s(p-1)); Re(s)> 0; in AP ******* ****** WITH A PRECISION OF AT LEAST 19 DECIMAL DIGITS ******* -------------------------- Accuracy parameters: default precision is with prec = 19 prec = 19; A = 3200; P = 9600; K = 8; J = 4 Accuracy is 10^(-29) err1 = 1.8160335279673509769170764563 E-32 err2 = 1.3862135372243990823444496355 E-32 toterr = 3.2022470651917500592615260917 E-32 -------------------------- RESULTS -------------------------- lazygamma_AP(3;1;1) = 0.088697734043610556704192075157 lazygamma_AP(3;2;1) = 0.48356682867672584922258373801 -------------------------- -------------------------- ****** COMPUTATION TIME ******** Computation time: 0 min, 0 sec, 139 millisec ****** END PROGRAM ******** ? lazygamma_AP(5,1) ****** COMPUTATION OF sum_{p = a mod d} log p/ (p^s(p-1)); Re(s)> 0; in AP ******* ****** WITH A PRECISION OF AT LEAST 19 DECIMAL DIGITS ******* -------------------------- Accuracy parameters: default precision is with prec = 19 prec = 19; A = 3200; P = 9600; K = 8; J = 4 Accuracy is 10^(-29) err1 = 0.000000000000000000000000000000018160335279673509769170764563 err2 = 0.000000000000000000000000000000013862135372243990823444496355 toterr = 0.000000000000000000000000000000032022470651917500592615260917 -------------------------- RESULTS -------------------------- lazygamma_AP(5;1;1) = 0.032215765502861114070237406669 lazygamma_AP(5;2;1) = 0.41190952219648495888544674866 lazygamma_AP(5;3;1) = 0.21326981117369976623624394485 lazygamma_AP(5;4;1) = 0.017499616336937163237350619148 -------------------------- -------------------------- ****** COMPUTATION TIME ******** Computation time: 0 min, 0 sec, 285 millisec ****** END PROGRAM ******** ? lazygamma_AP(5,1+I) ****** COMPUTATION OF sum_{p = a mod d} log p/ (p^s(p-1)); Re(s)> 0; in AP ******* ****** WITH A PRECISION OF AT LEAST 19 DECIMAL DIGITS ******* -------------------------- Accuracy parameters: default precision is with prec = 19 prec = 19; A = 3200; P = 9600; K = 8; J = 4 Accuracy is 10^(-29) err1 = 0.000000000000000000000000000000018160335279673509769170764563 err2 = 0.000000000000000000000000000000013862135372243990823444496355 toterr = 0.000000000000000000000000000000032022470651917500592615260917 -------------------------- RESULTS -------------------------- lazygamma_AP(5;1;1 + I) = -0.021444216980944757854997960099 - 0.0093785312310658457543520028333*I lazygamma_AP(5;2;1 + I) = 0.23650242252504415423045698513 - 0.26247577449567019857128800777*I lazygamma_AP(5;3;1 + I) = 0.061148727522832592503673327855 - 0.16700786331834430922134334949*I lazygamma_AP(5;4;1 + I) = -0.012533372292999425603456557975 + 0.0026226846602237623863835390936*I -------------------------- -------------------------- ****** COMPUTATION TIME ******** Computation time: 0 min, 0 sec, 914 millisec ****** END PROGRAM ******** *****************/