A note on the computation of the Euler-Kronecker constants for cyclotomic fields

The goal of this note is to introduce an alternative method to compute the Euler-Kronecker constants for cyclotomic fields and to compare it with other two different ways of computing the same quantity. The new algorithm requires the values of the generalised gamma functions $\Gamma_1$, also known as ${_2}\Gamma$, at some rational arguments in $(0,1)$. Using such method we were able to get $EK_{964477901}= -0.182374\dotsc$, thus giving an independent confirmation of Theorem 4 of Ford-Luca-Moree, $EK^+_{964477901} = 10.402223\dotsc$, and to compute the values of $EK_q$ and $EK^+_q$ for every odd prime $q\le 100000$. We also computed the value of $EK_q$ for some larger prime number $q$ but with no success in finding another negative value. Moreover, as a by-product, we will also provide more data on the generalised Euler constants in arithmetic progressions. The programs used to performed the computations here described and the numerical results obtained are available at the following web address: \url{http://www.math.unipd.it/~languasc/EK-comput.html}.


1.
The goal of this note is to introduce an alternative method to compute the Euler-Kronecker constants for cyclotomic fields and to compare it with other two different ways of computing the same quantity. Moreover, as a by-product, we will also provide more data on the generalised Euler constants in arithmetic progressions.
The definition of the Euler-Kronecker constant for number fields is given in section 1.3 of Ford-Luca-Moree [7], see eq. (1.14) there, but we are here just interested in the special case of cyclotomic fields. In this situation we can use eq. (2.6) of [7]: if q is an odd prime then we define the Euler-Kronecker constant for the cyclotomic field Q(ζ q ) as where ζ q is a primitive q-root of unity, γ is the Euler constant, χ are the non-trivial Dirichlet characters mod q and χ 0 is the trivial Dirichlet character mod q. In [7] the quantity EK q is denoted as γ q but this conflicts with notations used in literature. As we will see later, other quantities related to EK q are the generalised Euler constants in arithmetic progressions, sometimes also called Stieltjes constants in arithmetic progressions, denoted as γ k (a, q), k ∈ N, q ≥ 1, 1 ≤ a ≤ q, which are defined by γ k (a, q) := lim N→+∞ 0<m≤N m≡a mod q (log m) k m − (log N) k+1 q(k + 1) by, e.g., eq. (3)-(4) of Bohman-Fröberg [1]. Remark that γ 0 = γ = 0.577215664901 . . . , the Euler-Mascheroni constant. The quantities in (2) and, as we will see in sections 2-3 below, the one in (1), are hence connected with the values of ψ n , n ≥ 1, which is the logarithmic derivative of Γ n , a generalised Gamma function, see Deninger [4], Dilcher [6] and Katayama [10], whose definition for n = 1 is given in section 3.2. In some sense we can say that the ψ n -functions, n ≥ 1, are the analogue of the usual digamma function. In the following we will denote as ψ the standard digamma function Γ /Γ; we also remark that it can be represented as the function ψ 0 defined in (3).
The theoretical part about the proofs of the identities that we will use in sections 2-3 is classical; a key tool is the functional equation for the Dirichlet L-functions. Such proofs can be found in Cohen's books [2]- [3], for instance. Other useful references are also the papers of Dilcher [5,6] and Deninger [4].
where T(x) is defined in (5) (pay attention to the change of sign in (5) with respect to eq. (3.2) of [7]). Summarising, we finally get Formula (7), which was the one used in the paper by Ford-Luca-Moree [7], let us see that we can compute EK q via (1) using the values of ψ(a/q) and T(a/q), which is connected to ψ 1 (a/q) via (5), together with the values of the non-trivial Dirichlet characters mod q.
From a computational point of view it is clear that in (7) we first have to evaluate T(a/q) and ψ(a/q) for every 1 ≤ a ≤ q − 1. For the ψ(a/q)-values we can rely on the PARI/Gp function psi while, for the T(a/q)-values we can use the summing function sumnum. We'll see more on these computations in section 4.

A : D '
3.1. χ χ 0 is a primitive odd Dirichlet character. Recall that q is an odd prime, let χ χ 0 be a primitive odd Dirichlet character mod q and let τ( χ) := q a=1 χ(a) e(a/q), e(x) := exp(2πix), be the Gauss sum associated with χ. The functional equation for L(s, χ), see, e.g., the proof of Theorem 3.5 of Gun-Murty-Rath [9], gives and hence L (s, χ) which, evaluated at s = 0, gives By the Lerch identity about values of the Hurwitz zeta-function, see, e.g., either eq. (3.1) of Gun-Murty-Rath [9] or Proposition 10.3.5 of Cohen [3], and the orthogonality of Dirichlet characters, we get since, see Corollary 10.3.2 of Cohen [3], we have L(0, χ) = −( q−1 a=1 a χ(a))/q. Summarising we obtain From a computational point of view, in (8) we need to compute the log Γ(a/q) -values instead of the ψ(a/q) ones as in (7); to do so we can rely on internal PARI/Gp functions. In the next paragraph we will see that we have to reuse such values for the even Dirichlet characters case.
We finally remark that the computation in this section reveals that the Euler-Kronecker constant EK + q for the maximal real subfield of Q(ζ q ) is directly connected with the S-function values at a/q since, according to eq. (10) of Moree [11] and (11), we have and hence it seems that in this case the relevant information is encoded in the S-function. By (8) it is then trivial to get that First of all we notice that PARI/Gp, v. 2.11.1, has the ability to generate the Dirichlet L-functions (and many other L-functions) and hence the computation of EK q can be performed using (1) with few instructions of the gp scripting language. This computation has a linear cost in the number of calls of the lfun function of PARI/Gp and, at least for 271 ≤ q ≤ 20011, is, on our Macbook laptop, slower than both the approaches we are about to describe below.
Comparing (8) and (11) with (7), we see that in both cases we can rely on internal PARI/Gp functions to compute either the log Γ(a/q) -values or the ψ(a/q)-values, 1 ≤ a ≤ q − 1, and finally we have to evaluate the T and S series respectively involved. Remarking that all these functions have a pole in 0 and that we will take q very large, it is also relevant to know their order of magnitude for Hence for x → 0 + , we can expect that log Γ(x) and S(x) will be exponentially smaller than ψ(x) and T(x). Another difference is that, for the odd Dirichlet characters, eq. (8) before does not involve the estimation of an infinite series. So it seems reasonable to compare the following two approaches: a) use the T-series formulae and the ψ(a/q)-values for computing both the odd and the even Dirichlet characters cases like in [7]; b) use the S-series formulae for the even case and the finite sum of the a χ(a)-values for the odd one; remark that in both cases we have to evaluate a sum of the log Γ(a/q) -values. This way we can double check the computation performed in [7] not only because we are developing a different implementation of the same formulae, but also because we can approach the problem in an alternative way. In the computation we will use the PARI/Gp scripting language to exploit its ability in accurately evaluate the infinite sums involved in the definition of the T and S series previously described via the predefined summing function sumnum.
To check the correctness of the practical computations it is possible to use the following formulae; recalling γ = 0.577215664901 . . . and ζ (0) = −2.006356455908 . . . , we have that Formula (12) follows from Gauß' multiplication formula, see, e.g, section 12.15 of Whittaker-Watson [13], formula (13) is an immediate consequence of Theorem 2.5 of Deninger [4] and formulae (14)-(15) follow respectively from equations (7.9)-(7.10) of Dilcher [5]. We also remark that approaches a)-b) trivially require a quadratic number of products to perform the computations in (7)-(8) and (11), but this can be improved by using the Discrete Fourier Transform (DFT) and the following argument. Focusing on (7)-(8) and (11), we remark that, since q is prime, it is enough to get g, a primitive root of q, and χ 1 , the Dirichlet character mod q given by χ 1 (g) = e 2πi/(q−1) , to see that the set of the non-trivial characters mod q is { χ j 1 : j = 1, 2, . . . , q − 2}. Hence, if, for every k ∈ {0, . . . , q − 2}, we denote g k ≡ a k ∈ {1, . . . , q − 1}, every summation in (7)-(8) and (11) is either of the type q−1 k=1 e 2πi j k/(q−1) f (a k /q) or q−1 k=1 e −2πi j k/(q−1) f (a k /q), where j ∈ {1, . . . , q − 2} and f is a suitable function. As a consequence, such quantities are the DFT, or its inverse transformation, of the sequence { f (a k /q) : k = 0, . . . , q − 2}. This was used in [7] to speed-up the computation of these quantities via the use of DFT-dedicated software libraries. Unfortunately in the scripting language of PARI/Gp the DFT-functions work only if q = 2 + 1, for some ∈ N. So we had to trivially perform these summations and hence, in practice, this part is the most time consuming one in both the approaches a) and b) since it has a quadratic cost in q; this is the reason why for q > 2011 the direct approach using the lfun function of PARI/Gp becomes faster than the others (trivially performing the sums over a = 1, . . . , q − 1).
Being aware of such limitations, we performed the computation of EK q with these three approaches for every q prime, q ≤ 300, on a Macbook Air laptop ("Early 2015", 8Gb RAM, 1.6 Ghz Intel Core i5, two cores) using a precision of 30 decimal digits, see Table 1 of section 6. The results coincide up the desired precision and are coherent with the data of Table 1 on page 1472 of [7]. Moreover it seems that the version which uses the T-function is a bit faster than the one with the S-function probably because the internal sequence involved in its summation has a simpler form with respect to the one involved in the sum defining S, although the general term of such series have roughly the same decay order. In fact the computation of the values of Table 1 needed 1 minute and 7 seconds using the T-function, 1 minute and 18 seconds using the S-function and 1 minute and 23 seconds via the direct approach. We also computed the values of EK q for q = 1009, 2003, 3001, 4001, 5003, 6007, 7001, 8009, 9001, 10007, 20011, 30011, as you can see in Table 2 of section 6. These numbers were chosen to heuristically evaluate how the computational cost depends on the size of q. In this case, in the fifth column of Table 2 we also reported the running time of the direct approach, i.e. using (1), the third and fourth columns are respectively the running times of the approaches a) and b), while the sixth one indicates the used precision. For these values of q it became clear that the computation time spent in performing the sums having the Dirichlet characters values as coefficients was the longest one. This means that inserting a DFT-algorithm in the approaches a) and b) is fundamental to further improve their performances.
Hence for larger values of q we used the gp2c compiler tool to obtain suitable C programs to perform the precomputations of the needed T, ψ, S and log(Γ)-values. Then we passed such values to other C programs which used the fftw [8] library to perform the final stage. In this final stage the performances were extremely good in the sense that such a part was thousands-times faster than the same one trivially performed. This way we computed the values of EK q for q = 40009, 42611, 50021, 60013, 70001, 80021, 90001, 100003, 305741, 1000003,  4178771, 6766811, 10000019, 28227761, with double, long double and quadruple precisions, see Table 3. Some of these q-values were chosen for their dimension, while others because their measures using the "greedy sequence of prime offsets", http://oeis.org/A135311, are larger than 1.2 so that they are good candidates to have a negative Euler-Kronecker constant, see sect. 1.4 of [7] or sect. 2-3 of [11]. Such computations were performed with and Dell OptiPlex-3050, equipped with an Intel i5-7500 processor, 3.40GHz, four cores, 8 GB of RAM and running Ubuntu 18.04. We remark that the quadruple precision computation performances are affected from a lack of hardware support of the FLOAT128 type of the C programming language.
After having evaluated the running times of the previous examples, we decided to provide the scattered plots, see Figures 1-2, of the normalised values of EK q and EK + q (both in long double precision) for every odd prime q ≤ 10 5 thus doubling the known range of such data, see [7]. Such computation required the use of the cluster of the Mathematical Department of the University of Padova; for a description of the used cluster see http://www.math.unipd. it/~languasc/EKcomput/Description-Cluster-Math-Unipd.pdf. The minimal value of EK q /log q, 3 ≤ q ≤ 10 5 , q prime, is 0.23449 . . . and it is attained at q = 42611, as expected.
For even larger values of q the precomputations, if performed on a single desktop computer, would require too much time; hence we parallelised them on the cluster previously mentioned. This way we were able to obtain an independent confirmation of Theorem 4 of Further improvements on our programs were then performed to lower the RAM and disk occupation and to use a dedicated FFTW interface, called guru64, for being able to work on transforms whose length is ≥ (2 32 −1)/2 (the int bound for the C programming language). This let us to evaluate the Euler-Kronecker constants for larger "good" candidates q (in the sense that their measures using the greedy sequence of prime offsets were larger than 1.2). This way, after having used the cluster to get the values of T, ψ, S and log(Γ), in about 90 minutes of computation time on the same machine mentioned before we got that EK 2918643191 = 0.302789 . . . and EK + 2918643191 = 12.573983 . . . using the long double precision. In this case it seems that double precision version performed using the T-function is much less stable than the one with the S-function probably because of the fact that T(x) and ψ(x) are, for x → 0 + , respectively much larger that S(x) and log(Γ(x)).
The PARI/Gp scripts and the C programs used and the computational results obtained are available at the following web address: http://www.math.unipd.it/~languasc/EK-comput. html.

O E γ k (a, q)
Recall that q is an odd prime. In the case we are using the T-series, we have to precompute their values at a/q and the ones of ψ at the same arguments. Hence, as a by-product we can also obtain the values of the generalised Euler constants γ 0 (a, q) and γ 1 (a, q) as you can see in paragraphs 5.1-5.2. In practice this is obtained by activating an optional flag in the main gp script. The case about γ k (a, q), k ≥ 2, is described in paragraph 5.3. γ 0 (a, q). For γ 0 (a, q) with 1 ≤ a ≤ q − 1, q odd prime, by (2) we have

Generalised Euler constants
Moreover, since ψ(1) = −γ and T(1) = 0, we also have Using the formulae in the previous two paragraphs we computed γ 0 (a, q) and γ 1 (a, q) with q prime, 3 ≤ q ≤ 100, 1 ≤ a ≤ q, in about 9 seconds of computation time with a precision of 30 digits. Such results are listed at the bottom of the gp-script file that can be downloaded here: http://www.math.unipd.it/~languasc/EK-comput.html.

5.
3. The general case γ k (a, q), k ≥ 2. The general case γ k (a, q), k ∈ N, k ≥ 2, q ≥ 1, 1 ≤ a ≤ q, do not follow from the data already computed for the Euler-Kronecker constants since we need information about the values of ψ n (x), for every 2 ≤ n ≤ k. Such a direct computation of both ψ n (a/q) and γ n can be easily performed via eq. (2)-(3) using the PARI/Gp summing function sumnum paying attention to submit a sufficiently fast convergent sum. For example, to compute γ n , n ∈ N, we used the formulae and which both easily follow from (4). We get, in less than 15 seconds of time and with a precision of at least 40 digits, the results in Table 4 of section 6; to be sure about the correctness of such results we computed them twice using the formulae (16)-(17) and then we compared their results. Such values are in agreement with the data on page 282 of Bohman-Fröberg [1] for n = 0, . . . , 20. For larger n's the formulae in (16)-(17) seem to be not good enough to get precise results via the sumnum function with this precision level. To compute ψ n (a/q) and, as a consequence, γ k (a, q), we can proceed in a similar way as we did for T(a/q) and γ 1 (a, q), see the program file here http://www.math.unipd.it/~languasc/ EK-comput.html. At the bottom of such program file you can find a large list (too long to be included here) of computed values of γ k (a, q) for 1 ≤ k ≤ 20, 1 ≤ q ≤ 9, 1 ≤ a ≤ q, with a precision of 20 digits. In about 1 minute and 33 seconds of computation time we replicated Dilcher's computations since the values we got are in agreement with the data on pages S21-S24 of [5].
Acknowledgements. I wish to thank Karim Belabas and Bill Allombert for a couple of key suggestions about libpari and gp2c and Luca Righi for his help in developing the quadruple precision versions of the fft-programs, in designing the parallelised precomputations and in organising the computational runs on the cluster of the Math. Dept. of the University of Padova.