Fast multi-precision computation of some Euler products

TL;DR

This paper introduces a method for rapidly computing certain Euler products over subsets of the multiplicative group modulo q, enabling precise calculations relevant to number theory.

Contribution

It presents a novel algorithm for double exponential time computation of Euler products over specific subsets, with an implementation in Sage and extensions to polynomial ratios.

Findings

Efficient computation of Euler products for given subsets

Implementation of Sage script for practical calculations

Application to precise evaluation of number-theoretic Euler products

Abstract

For every modulus $q\ge3$, we define a family of subsets $\mathcal{A}$ of the multiplicative group $(\mathbb{Z}/{q}\mathbb{Z})^\times$ for which the Euler product $\prod_{p\text{mod}q\in\mathcal{A}}(1-p^{-s})$ can be computed in double exponential time, where $s>1$ is some given real number. We provide a Sage script to do so, and extend our result to compute Euler products $\prod_{p\in\mathcal{A}}F(1/p)/G(1/p)$ where $F$ and $G$ are polynomials with real coefficients, when this product converges absolutely. This enables us to give precise values of several Euler products intervening in Number Theory.

Equations

$$\alpha_{0}^{(3)}=\frac{1}{2^{1/2}3^{1/4}}\prod_{p\equiv 2[3]}\biggl{(}1-\frac{1}{p^{2}}\biggr{)}^{-1/2}$$

$$\beta_{0}=\frac{3^{1/4}\sqrt{\pi}}{2^{5/4}}\frac{\log(2+\sqrt{3})^{1/4}}{\Gamma(1/4)}\prod_{p\equiv 5,7,11[12]}\biggl{(}1-\frac{1}{p^{2}}\biggr{)}^{-1/2}.$$

$$N(x)={\alpha }_0^{(3)}\frac{x(1+o(1))}{\sqrt{\log x}}.$$

$$N^{\prime }(x)={\beta }_0\frac{x(1+o(1))}{(\log x{)}^{3/4}}.$$

$${\alpha }_0^{(3)}=0.$$

$$33454216923657080763100276496582468971791125286643\cdots$$

$${\beta }_0=0.$$

$$96454588678412888654487524105108994874678139792727\cdots$$

$$\prod _{p\in \mathcal{A}\text{mod}q}(1-p^{-s}{)}^{-1}$$

$$\zeta (s;q,\mathcal{A})=\prod _{p\text{mod}q\in \mathcal{A}}(1-p^{-s}{)}^{-1}$$

$$\zeta (s;q,\mathcal{A})=\prod _{\begin{array}{c}p\text{mod}q\in \mathcal{A},\\ p\ge P\end{array}}(1-p^{-s}{)}^{-1}.$$

$$L_P(s,\chi )=\prod _{p\ge P}(1-\chi (p)/p^s{)}^{-1}.$$

$${\gamma }_s(G_0,t)=\log \prod _{\chi \in G_0^{\perp }}{L}_P(ts,\chi ).$$

$${\Gamma }_s(t)=({\gamma }_s(G_0,t){)}_{G_0\in \mathcal{G}}.$$

$$V_{s}(t)=\bigl{(}\log\zeta_{P}(s;q,\mathcal{A})\bigr{)}_{\mathcal{A}\in G^{\sharp}}.$$

$$\Vert W\Vert =\underset{i}{\max }|W_i|$$

$$M_{1}^{-1}\big{|}_{i=\mathcal{A},j=K}=\begin{cases}\mu(|\langle{\mathcal{A}}\rangle{}/K|)/|G/K|&\text{when $K\subset\langle{\mathcal{A}}\rangle{}$},\\ 0&\text{otherwise}\end{cases}$$

$$N_{d}\big{|}_{i=B_{0},j=B_{1}}=\sum_{\begin{subarray}{c}K\subset B_{0},\\ |KB_{1}/K|=d\end{subarray}}\mu(|B_{0}/K|).$$

$$\biggl{\|}V_{s}(1)-\sum_{0\leq v\leq r-1}(-1)^{v}\sum_{\begin{subarray}{c}d_{1}\cdots d_{v}\leq 2^{r}\end{subarray}}\frac{N_{d_{1}}}{d_{1}}\ldots\frac{N_{d_{v}}}{d_{v}}M_{1}^{-1}\Gamma_{s}(d_{1}\ldots d_{v})\biggr{\|}\\ \leq\frac{1}{2}\biggl{(}1+\frac{r-1}{|G^{\sharp}|}\biggr{)}\biggl{(}\frac{|G^{\sharp}|d(\exp G)}{2}\biggr{)}^{r-1}\frac{1+P/(s2^{r}-1)}{P^{s2^{r}}}$$

$$(\log {\zeta }_P(s;q,\mathcal{A}|r){)}_{\mathcal{A}\in G^{♯}}=\sum _{0\le v\le r-1}(-1{)}^v\sum _{\begin{array}{c}{d}_1\cdots d_v\le 2^r\end{array}}\frac{N_{d_1}}{d_1}\dots \frac{N_{d_v}}{d_v}{M}_1^{-1}{\Gamma }_s(d_1\dots d_v)$$

$$\prod _{\begin{array}{c}p\ge P,\\ p\in \mathcal{A}\end{array}}\frac{F(1/p)}{G(1/p)}=\prod _{2\le j\le J}{\zeta }_P(j;q,\mathcal{A}|r{)}^{b_G(j)-b_F(j)}\times I,$$

$$|\log I|\leq\max(\deg F,\deg G)\biggl{(}\biggl{(}\frac{|G^{\sharp}|d(\exp G)}{2}\biggr{)}^{r-1}\frac{r\beta^{2}}{P^{2^{r+1}}}(1+2^{-r}P)+\frac{4\beta^{J+1}}{JP^{J}}\biggr{)}.$$

$$\beta=\max\Bigl{(}1,\sum_{1\leq k\leq\deg F}|a_{k}|,\sum_{1\leq k\leq\deg G}|b_{k}|\Bigr{)}$$

$$\prod_{p\equiv 1[8]}\biggl{(}1-\frac{4}{p}\biggr{)}\biggl{(}\frac{p+1}{p-1}\biggr{)}^{2}=\begin{aligned} 0.&95694\,53478\,51601\,18343\,69670\,57273\,89182\,87531\\ &74977\,2913914789\,05432\,60424\,60170\,16444\,88885\\ &94814\,40512\,03907\,95084\cdots\end{aligned}$$

$$I$$

$$=\begin{array}{rl}0.& 66974096993707122053892243157176440668837015743648\\ & 24185732985228452467999564571472731506210214359373\cdots \end{array}$$

$$\prod _{p\equiv 1[8]}\frac{p(p-8)}{(p-4{)}^2}=\begin{array}{rl}0.& 8830710047439466714178342990031085346768\\ & 8883488097347071929515939521194699065659\\ & 68857993832860379164\cdots \end{array}$$

$$\lambda$$

$$=\begin{array}{rl}0.& 79220823816754166877545556657902410112893225098622\\ & 11172279734525695141549441249066029538839802752927\cdots \end{array}$$

$$A_{\chi}=\prod_{p\geq 2}\biggl{(}1+\frac{(\chi(p)-1)p}{(p^{2}-\chi(p))(p-1)}\biggr{)}$$

$$\prod _{\begin{array}{c}p\in \mathcal{P},\\ f(p)=-1\end{array}}(1-p^{-s}{)}^2=\frac{{\prod }_{p\in \mathcal{P}}(1-p^{-s})}{{\prod }_{p\in \mathcal{P}}(1-f(p)p^{-s})}\prod _{\begin{array}{c}p\in \mathcal{P},\\ f(p)=-1\end{array}}(1-p^{-2s}).$$

$$\prod _{\begin{array}{c}p\in \mathcal{P},\\ f(p)=-1\end{array}}\frac{(1-p^{-s}{)}^2}{1-p^{-2s}}$$

Full text

1. Introduction

At the beginning of our query lie two constants that appear in the paper [2] by É. Fouvry, C. Levesque and M. Waldschmidt. On following this paper, they are

[TABLE]

and

[TABLE]

Both occur in number theory as densities. The number of integers $n$ of the shape $n=x^{2}-xy+y^{2}$, where $x$ and $y$ are integers (these are the so-called Loeschian numbers, see sequence A003136 of [10]) is given by

[TABLE]

This accounts for our interest in the first constant. The second one occurs because the number of Loeschian numbers that are also sums of two squares (see sequence A301430 of [10]) is given by

[TABLE]

The question we address here is devising a fast manner to compute the intervening Euler products. From sequence A301429 of [10], we know that $\alpha_{0}^{(3)}=0.638909\ldots$ but we would like (much!) more digits. Similarly it is known that $\beta_{0}=0.30231614235\ldots$.

Theorem 1.1**.**

We have

[TABLE]

and

[TABLE]

Our method is more general and allows one to compute Euler products of the shape

[TABLE]

for any $s$ with $\Re s>1$ and some subsets $\mathcal{A}$ of $(\mathbb{Z}/q\mathbb{Z})^{\times}$. We use a set of identities that lead to fast convergent formulae. The use of a similar formula for scientific computations can be found in [14, equation (15)] by D. Shanks. This author’s approach has been put in a general context by P. Moree & D. Osburn in [7, equation (3.2)]. On looking closely, we see that an accurate value of $\alpha_{0}^{(3)}$ already follows from this paper. The formulae we prove have a wider reach, though they fail to exhaust the problem. The reader may want to read subsection 2.2 now to understand the initial idea. In the simplest form, we produce a formula that links for instance $\zeta(s;12,1)=\prod_{\begin{subarray}{c}p\equiv 1[12]\end{subarray}}{(1-p^{-s})^{-1}}$ to $\zeta(2s;12,1)$. We then reuse this formula to change $2s$ in $4s$, and so on, and we finally use $\zeta(2^{r}s;12,1)=1+\mathcal{O}(1/2^{s2^{r}})$. This is analogous to D. Shanks scheme in [14]. In the general case however, we link values at $s$ with values at $ds$ for some $d>1$, but these values are not the one of the same function, but of some companion functions. This means that we have to work simultaneously with several players. Let us first define these companions, which are all the products we propose to compute.

When $K$ is a cyclic subgroup of $(\mathbb{Z}/q\mathbb{Z})^{\times}$, we denote by $A(K)$ the set of elements $x$ from $(\mathbb{Z}/q\mathbb{Z})^{\times}$ such that the subgroup $\langle{x}\rangle{}$ generated by $x$ is equal to $K$. We note that the sets $A(K)$, when $K$ ranges though the set of cyclic subgroups of $(\mathbb{Z}/q\mathbb{Z})^{\times}$, determine a partition of $(\mathbb{Z}/q\mathbb{Z})^{\times}$. A subset $\mathcal{A}$ of $(\mathbb{Z}/q\mathbb{Z})^{\times}$ is said to be a lattice-invariant class if it is of the form $A(K)$ for some cyclic subgroup $K$ of $(\mathbb{Z}/q\mathbb{Z})^{\times}$, i.e. if all its elements generate the same subgroup (see Definition 3.1 below). Here is a consequence of our approach.

Theorem 1.2**.**

Let $q$ be some modulus and $\mathcal{A}$ be a lattice-invariant class of $(\mathbb{Z}/q\mathbb{Z})^{\times}$. For every $s>1$, the product

[TABLE]

can be computed in double-exponential time.

This theorem applies in particular to $\mathcal{A}=\{1\}$ and to $\mathcal{A}=\{-1\}$ and this is enough to compute $\beta_{0}$ and $\alpha_{0}^{(3)}$. The last section contains numerical examples. The material of this paper has been used to write the script

LatticeInvariantEulerProducts-02.sage

which we shorten below in LIEP.sage and which can be found on the second author website. We give some details about this script when developing the proof below.

We produce in Proposition 7.3 an explicit expression for the number $|G^{\sharp}|=|\mathscr{G}|$ of lattice-invariant classes. Though our formula is only a sum of non-negative summands that are multiplicative expressions, its order of magnitude is not obvious when $q$ has numerous prime factors. We have for instance not been able to establish that $|G^{\sharp}|\ll_{\epsilon}q^{\epsilon}$ (for every positive $\epsilon$) though this was our initial guess.

Notation

When $\mathcal{A}$ is a subset of $(\mathbb{Z}/q\mathbb{Z})^{\times}$, we define $\langle{\mathcal{A}}\rangle{}$ to be the (multiplicative) subgroup generated by $\mathcal{A}$, and when $\mathcal{A}=\{a\}$, we may shorten $\langle{\{a\}}\rangle{}$ in $\langle{a}\rangle{}$.

When additionally $P\geq 2$ is some real parameter, we define

[TABLE]

This is in accordance with the notation of Theorem 1.2. We define further

[TABLE]

Precise statement of the main result

Let $q>1$ be a modulus. Let $G_{0}$ be a subgroup of $G=(\mathbb{Z}/q\mathbb{Z})^{\times}$ and let $G_{0}^{\perp}$ be the subgroup of characters that take the value 1 on $G_{0}$. Let $s>1$ be a real number and $P\geq 2$ be a parameter. We shall compute directly the contribution of the primes $<P$. We define, for any positive integer $t$:

[TABLE]

The parameter $P$ has disappeared from our notation and the reader may stick with $P=2$. When $s$ is a real number, the number $\prod_{\chi\in G_{0}^{\perp}}L_{P}(ts,\chi)$ is indeed a positive real number because, when $\chi$ belongs to $G_{0}^{\perp}$, so does $\overline{\chi}$.

We denote the set of lattice-invariant classes by $G^{\sharp}$ and the set of cyclic subgroups of $G$ by $\mathscr{G}$. Both sets are in an obvious one-to-one correspondence. We consider the vector

[TABLE]

The rows of $\Gamma_{s}(t)$ are indexed by cyclic subgroups of $G$. It is computed by the function GetGamma of the script LIEP.sage from the values of the Hurwitz zeta function. See the implementation notes below. We next define

[TABLE]

The rows of $V_{s}(t)$ are indexed by classes. We control the size of our vectors with the norm

[TABLE]

when $W$ is the vector of coordinates $W_{i}$. We define the square matrix $M_{1}^{-1}$ by

[TABLE]

where $\mathcal{A}$ ranges $G^{\sharp}$ while $K$ ranges $\mathscr{G}$. It is unusual to define a matrix by its inverse. In the natural course of the proof, a matrix $M_{1}$ will occur, whose inverse is the one above; it is computed in Proposition 4.1. The reader will readily check that there are no circularity in our definitions. Let us recall that the exponent of $G$ is the maximal order of an element in $G$ and is denoted by $\exp G$. To each divisor $d>1$ of $\exp G$, we associate the square matrix $N_{d}$ whose columns and rows are indexed by cyclic subgroups of $G$ and whose entries are given by

[TABLE]

The sum is over subgroups $K$. The condition $|KB_{1}/K|=d$ can be replaced by the condition $|B_{1}/K\cap B_{1}|=d$. Here is our main theorem.

Theorem 1.3**.**

For any integer $r\geq 2$, we have

[TABLE]

where $d_{1},\ldots,d_{r}$ are all divisors of $\exp G$ excluding 1.

When $v=0$, we use $d_{1}\ldots d_{v}=1$ and $N_{d_{1}}\cdots N_{d_{v}}=\operatorname{Id}$. We provide in Section 5 the numerical datas modulo 7 that will enable the reader to follow the proof step by step in this case. This example may also be used to check our routines.

Extending the computations

Now that we know how to compute some Euler products $\zeta_{P}(s;q,\mathcal{A})$ in a fast manner, we can extend these computations to more general Euler products, though still on the same sets of primes. To do so, we add a definition:

[TABLE]

Theorem 1.4**.**

Let $F,G\in\mathbb{R}[X]$ be two coprime polynomials satisfying $F(0)=G(0)=1$ such that $(F(X)-G(X))/X^{2}\in\mathbb{R}[X]$. Let $\beta\geq 2$ be an upper bound for the maximum modulus of the inverses of the roots of $F$ and of $G$. Let $P\geq 2\beta$ be a parameter. Then, for any parameters $J\geq 3$ and $r\geq 2$, we have

[TABLE]

where the integers $b_{G}(j)$ and $b_{F}(j)$ are defined in Lemma 6.1 and

[TABLE]

Remark 1.5*.*

Inequality (53) gives a more precise bound for $|\log I|$ which we will use in the actual script.

Remark 1.6*.*

Lemma 6.3 ensures that we may select

[TABLE]

when $F(X)=1+a_{1}X+\ldots+a_{\delta}X^{\delta}$ and $G=1+b_{1}X+\ldots+b_{\delta^{\prime}}X^{\delta^{\prime}}$.

Remark 1.7*.*

The function GetEulerProds(q, F, G, nbdecimals) gives all these Euler products. The polynomials $F$ and $G$ are to be given as polynomial expressions with the variable $x$.

D. Shanks in [15] (resp. [16], resp. [17]) has already been able to compute an Euler product over primes congruent to 1 modulo 8 (resp. to 1 modulo 4, resp. 1 modulo 8), by using an identity (Lemma of section 2 for [15], equation (5) in [16] and the Lemma of section 3 in [17]) that is a precursor of our Lemma 6.1.

In these three examples, the author has only been able to compute the first five digits, and this is due to three facts: the lack of interval arithmetic package at that time, the relative weakness of the computers and the absence of a proper study concerning the error term. We thus complement these results by giving the first hundred decimals.

Corollary 1.8** (Shank’s Constant).**

We have

[TABLE]

And thus Shank’s constant satisfies

[TABLE]

As explained in [16], the number of primes $\leq X$ of the form $m^{4}+1$ is conjectured to be asymptotic to $I\cdot X^{1/4}/\log X$. The name “Shank’s Constant” comes from Chapter 2, page 90 of [1]. When using the script that we introduce below, this value is obtained with the call

GetEulerProds(8, $1-2*x-7*x^{2}-4*x^{3}$, $1-2*x+x^{2}$, 150, 400).

Corollary 1.9** (Lal’s Constant).**

We have

[TABLE]

And thus Lal’s constant satisfies

[TABLE]

As explained in [17], the number of primes $\leq X$ of the form $(m+1)^{2}+1$ and such that $(m-1)^{2}+1$ is also a prime is conjectured to be asymptotic to $\lambda\cdot X^{1/2}/(\log X)^{2}$. The name “Lal’s Constant” comes from the papers [5] and [17]. When using the script that we introduce below, the first value is obtained with the call

GetEulerProds(8, $1-8*x$, $1-8*x+16*x^{2}$, 100, 400).

We close this section by mentioning another series of challenging constants. In [8], P. Moree computes inter alia the series of constants $A_{\chi}$ defined six lines after Lemma 3, page 452, by

[TABLE]

where $\chi$ is a Dirichlet character. Our theory applies only when $\chi$ is real valued,

Thanks

The authors thank M. Waldschmidt for having drawn their attention of this question, P. Moree and É. Fouvry for helpful discussions on how to improve this paper and X. Gourdon for free exchanges concerning some earlier computations.

2. A general mechanism

We start by presenting the mechanism of Shanks in [13] is a general setting.

Lemma 2.1**.**

Let $\mathcal{P}$ be a set of prime numbers and let $f$ be a function from $\mathcal{P}$ to $\{\pm 1\}$. For every $s$ with $\Re s>1$, we have

[TABLE]

Proof.

The proof is straightforward. We simply write

[TABLE]

as required. ∎

Shanks’s method is efficient to deal with product of primes belonging to a coset modulo a quadratic character. We generalize it as follows.

Lemma 2.2**.**

Let $q>1$ be a modulus. We set $G_{0}$ be a subgroup of $G=(\mathbb{Z}/q\mathbb{Z})^{\times}$ and $G_{0}^{\perp}$ be the subgroup of characters that take the value 1 on $G_{0}$. For any integer $b$, we define $\langle{b}\rangle{}$ to the the subgroup generated by $b$ modulo $q$. We have

[TABLE]

and, for any element $a\notin G_{0}$ of order 2, we have

[TABLE]

where $\hat{G}$ is the set of characters of $G$.

Case $G_{0}=\{1\}$ of the first identity is classical in Dedekind zeta function theory, and can be found in [12, Proposition 13] in a rephrased form. Case $a\neq 1$ will not be required for the general theory. It may however lead quickly to efficient formulae.

Proof.

We note that $\prod_{\chi\in G_{0}^{\perp}}(1-\chi(p)z)^{\chi(a)}=\prod_{\psi\in\hat{H}}(1-\psi(p)z)^{f(\psi)}$ when $\langle{p}\rangle{}=H$ and where

[TABLE]

The condition $\chi\in G_{0}^{\perp}$ can also be written as $\chi|G_{0}=1$, hence we can assume that $\psi|(H\cap G_{0})=1$. We write

[TABLE]

where

[TABLE]

When $a$ lies outside $H{}G_{0}$, this sum vanishes; otherwise it equals $|G/(H{}G_{0})|\psi^{\prime}(a)$. The characters of $H{}G_{0}$ that are trivial on $G_{0}$ are canonically identified with the characters of the cyclic group $(H{}G_{0})/G_{0}$. We thus have

[TABLE]

and this proves our first formula.

When $a^{2}\equiv 1[q]$ and $a\notin G_{0}$, and since $(H{}G_{0})/G_{0}$ is cyclic, of (even) order $h$ say, the characters are given by $\chi(p^{x})=e(cx/h)$ since $p$ is a generator and where $c$ ranges $\{0,\cdots,h-1\}$. We thus have, when $a\in H$,

[TABLE]

The reader will readily complete the proof by setting $K=H{}G_{0}$. ∎

2.1. A special case

Let us select for $G_{0}$ the kernel of a given quadratic character $\chi_{1}$. The subgroup $K$ can take only two values, $G_{0}$ or $G$. We thus get

[TABLE]

which gets converted into

[TABLE]

Lemma 2.1 can also be used to obtain the same result.

2.2. More details modulo 12

Here is the character table modulo 12:

[TABLE]

First Identity

This table enables us to write:

[TABLE]

And thus

[TABLE]

which gives rise to the formula

[TABLE]

This identity reduces the computation of $\zeta(s;12,1)$ to the one of $\zeta(2s;12,1)$ and we can iterate this formula. Note that we can take the required fourth root as only real numbers are involved, when the terms are properly grouped.

Second Identity

Similarly, we find that

[TABLE]

whence

[TABLE]

which we finally write in the form

[TABLE]

This identity again reduces the computation of $\zeta(s;12,11)$ to the one of $\zeta(2s;12,11)$ and we can iterate this formula. Again, we can take the required fourth rooths as only real numbers are involved, when the terms are properly grouped.

Third Identity

We also find that

[TABLE]

whence

[TABLE]

We are exactly in the same position as with the second identity. We again finally write in the form

[TABLE]

This identity again reduces the computation of $\zeta(s;12,5)$ to the one of $\zeta(2s;12,5)$ and we can iterate this formula. Again, we can take the required fourth rooths as only real numbers are involved, when the terms are properly grouped.

Fourth Identity

We can easily produce a similar formula linking $\zeta(s;12,7)$ to $\zeta(2s;12,7)$ or use the fact that the product $\zeta(s;12,1)\zeta(s;12,5)\zeta(s;12,7)\zeta(s;12,11)$ equals $L(s,\chi_{0,12})$, and thus is known, to infer such a formula from the ones above.

3. Products obtained in general

We want to compute Euler products of the shape $\zeta(s;q,\mathcal{A})$ for $s>1$ and some subset $\mathcal{A}$ of $(\mathbb{Z}/q\mathbb{Z})^{\times}$. Computing $L(s,\chi)$ is easier as it can be reduced to sums over integers is some arithmetic progressions. Equation (17) reduces in a special case the computations of $\zeta(s;q,\mathcal{A})$ to the one of $\zeta(s;q,\mathcal{A})$, and we can continue the process. We soon reach $\prod_{p\in\mathcal{A}\mathbin{\,\textrm{mod}\,}q}(1-1/p^{2^{N}s})$ with a large enough $N$ which can be approximated by $1+\mathcal{O}(2^{-2^{N}s})$. The object of this section is to devise a setting to understand which sums we relate together.

Definition 3.1**.**

Two elements $g_{1}$ and $g_{2}$ of the abelian group $G$ are said to be lattice-invariant if and only if they generates the same group.

The map between the set of cyclic subgroups of $G$ and the set of lattice-invariant-classes which, to a subgroup, associates the subset of its generators, is one-to-one.

The function GetLatticeInvariantClasses of the script LIEP.sage gives the two lists: the one of the cyclic subgroups and the one of their generators, ordered similarly and in increasing size of the subgroup.

Any two elements of $(\mathbb{Z}/q\mathbb{Z})^{\times}$ equivalent according to it cannot be distinguished by using the formulae of Lemma 2.2. Conversely, the question is to know whether we are indeed able to distinguish each class. To each class $\mathcal{A}$, we attach the enumerable collection of symbols $(x_{\mathcal{A}}^{r})_{r\geq 1}$. We shall replace each of them according to the rule

[TABLE]

We consider the module of finite formal combinations

[TABLE]

with coefficients $\alpha_{{\mathcal{A}},r}\in\mathbb{Z}$ and indeterminates $x_{\mathcal{A}}^{r}$. The superscript $r$ is not a power. We consider the following special elements. Let $G_{0}\subset K\subset G$ be two subgroups such that $K/G_{0}$ is cyclic. We define

[TABLE]

With that, we find that

[TABLE]

4. Iterating the formula

The first identity of Lemma 2.2 gives us as many identities as there are subgroups $G_{0}$; we know by Definition 3.1 that the number of lattice-invariant-classes equals the one of cyclic subgroups. It turns out that it is enough to restrict our attention to cyclic subgroups $G_{0}$. Let $\mathscr{G}$ be the subset of such subgroups, which we order by inclusion. On recalling definition (7), we may rewrite (31) in the form

[TABLE]

where (this is the case $K=G_{0}$)

[TABLE]

and, where, when $d>1$ (i.e. $G_{0}\subsetneq K$), we have

[TABLE]

Equation (32) gives us a relation between $M_{1}V_{s}(t)$ and $M_{d}V_{s}(dt)$ for several $d$’s that are strictly larger than 1. Our roadmap is to invert the matrix $M_{1}$ and to iterate this formula. We compute explicitly $M_{1}^{-1}$ by using some generalised Moebius inversion, which we first put in place.

The Moebius function associated to $\mathscr{G}$

We follow closely the exposition of Rota in [11]. On the algebra of functions $f$ on couples $(K,L)$ of points of $\mathscr{G}$ such that $K\subset L$ (the so-called incidence algebra, see [11, Section 3]), we define the convolution product

[TABLE]

We consider the $\mathscr{G}$-zeta function which is defined by

[TABLE]

This function is shown to be invertible in the above algebra and its inverse is called the $\mathscr{G}$-Moebius function, denoted by $\mu_{\mathscr{G}}$. By definition, we have the two Moebius inversion formulas:

[TABLE]

and

[TABLE]

We end this reminder with a formula giving the value of $\mu_{\mathscr{G}}(K,H)$.

Computing $\mu_{\mathscr{G}}(K,H)$

Let $C_{p}(K,H)$ be the number of chains of length $p$ going from $K$ to $H$, i.e. the number of $p+1$-uples $K=A_{0}\subsetneq A_{1}\subsetneq A_{2}\subsetneq\ldots\subsetneq A_{p}=H$. Then (cf [11, Proposition 6])

[TABLE]

Since the subgroups of a cyclic group are all cyclic, we only have to consider the chains in $H/K$. There is one and only one subgroup for each divisor of $|H/K|$, and any two such subgroups $L_{1}$ and $L_{2}$ are included according to whether |L_{1}|\big{|}|L_{2}| or not. This transfers the problem on a problem on integers. Let $c_{\ell}(n)$ be the number of $\ell+1$-divisibility chains between 1 and $n$. We have $c_{0}(n)=1\!\!\!1_{n=1}$ while $c_{1}(n)=1\!\!\!1_{n\geq 2}$ and $c_{p+1}(n)=(c_{\ell}\star c_{1})(n)$. This proves that $c_{\ell}(n)=d_{\ell}^{*}(n)$, the number of $p$-tuples $(d_{1},d_{2},\ldots,d_{\ell})$ of divisors of $n$ that are such that $d_{i}\neq 1$ and $d_{1}d_{2}\cdots d_{\ell}=n$. We have

[TABLE]

and thus the generating series of $\sum_{p\geq 0}(-1)^{p}d^{*}_{\ell}(n)$ is

[TABLE]

We have proved that

[TABLE]

Inverting the matrix $M_{1}$

Proposition 4.1**.**

The matrix $M_{1}$ is invertible and the coefficients of its inverse are given by

[TABLE]

Proof.

We find that

[TABLE]

We replace ${\mathcal{A}}$ by the subgroup $B=\langle{{\mathcal{A}}}\rangle{}$ it generates. Inverting $f(K)=|G/K|\sum_{B\subset K}v_{B}$ is done with the Moebius function of $\mathscr{G}$. To do so, simply consider the more general function

[TABLE]

where $v^{*}(H,B)=v_{B}$. This gets inverted in

[TABLE]

which yield, by specializing $H=\{1\}$

[TABLE]

We could also have applied [11, Proposition 2 (**)]. This gives us

[TABLE]

Our proposition is proved. ∎

The function GetM1Inverse of the script LIEP.sage computes $M_{1}^{-1}$.

The recursion formula

We start from (32) and deduce that

[TABLE]

We readily find that $N_{d}=dM_{1}^{-1}M_{d}$ is given by (11).

Proof.

Indeed we have

[TABLE]

This is exactly what we have written in (11). ∎

By considering the exact sequence

[TABLE]

one shows that $|KB_{1}/K|=|B_{1}|/|K\cap B_{1}|$. As a consequence, we see that only the $d$ that divides the exponent of $G$ appear. The function GetNds of the script LIEP.sage computes $(N_{d})_{d}$.

(41)

$V_{s}(t)=-\sum_{\begin{subarray}{c}d|\exp G,\\ d\neq 1\end{subarray}}\frac{N_{d}}{d}V_{s}(dt)+M_{1}^{-1}\Gamma(t).$

Unfolding the recursion

Let $z\geq 1$ and $r\geq 1$ be two parameters. We have

[TABLE]

where $d_{1},\ldots,d_{r}$ are all divisors of $\exp G$ excluding 1. We can incorporate the last summand in the one before by considering as the value for $s=0$.

Proof.

Let us prove this formula by recursion. Case $r=1$ is just (41). Let us see precisely what happens for $r=2$. We start from

[TABLE]

which we rewrite as

[TABLE]

We use again this equation on $V_{s}(d_{1}t)$ when $d_{1}\leq z$, and $z/d_{1}$ rather than $z$, getting

[TABLE]

To go from $r$ to $r+1$, we select the divisors $d_{r}$ that are such that $d_{1}d_{2}\cdots d_{r}\leq z$ and employ (41) on $V_{s}(d_{1}\cdots d_{r}t)$. ∎

Lemma 4.2**.**

The coefficients of a product $N_{d_{1}}N_{d_{2}}\cdots N_{d_{v}}$ are at most (in absolute value) equal to $|G^{\sharp}|^{v-1}$, where $G^{\sharp}$ is the set of lattice-invariant classes (which is also the number of cyclic subgroups of $G$).

End of the proof of Theorem 1.3

The formula (47) with $t=1$ contains most of our proof. We only have to control the error term, which is our next task.

The number of possible $d$’s is at most the number of divisors of $\exp G$ minus 1, so at most $d(\exp G)$. The coefficients of a typical product $N_{d_{1}}\cdots N_{d_{v}}$ are of size at most $|G^{\sharp}|^{v-1}$, we divide each coefficient by $d_{1}\cdots d_{v}$ which is at least $z$, and we have at most $d(\exp G)^{v}$ $v$-tuples $(d_{1},\ldots,d_{v})$. As a consequence, each coordinate, says $y$, of the vector

[TABLE]

satisfies

[TABLE]

We deal similarly with the coordinates of the vector

[TABLE]

except that the denominator $d_{1}\cdots d_{r}$ is not especially larger than $z$; we however select $z=2^{r}$ to ensure this condition. This means that only $d_{1}=d_{2}=\ldots=d_{r}=2$ is admissible. So, on combining both, we see that

[TABLE]

To complete the proof, we simply need a bound for $\max_{D\geq 2^{r}}\|V_{s}(D)\|$ and such a bound is provided by the next lemma.

Lemma 4.3**.**

Let $\mathcal{A}$ be a subset of the $G=(\mathbb{Z}/q\mathbb{Z})^{\times}$. Let $f>1$ be a real parameter. We have

[TABLE]

Proof.

We use

[TABLE]

hence, by using a comparison to an integral, we find that

[TABLE]

∎

5. A detailed example modulo 7

Wet set $G=(\mathbb{Z}/7\mathbb{Z})^{\times}$. We find that

[TABLE]

(indexed in this order) and that

[TABLE]

also indexed in that order. There are 6 Dirichlet characters whose values are given by (with $\zeta_{6}=\exp(2i\pi/6)$)

[TABLE]

We obtain this list with the command

[TABLE]

and the remark $\zeta_{6}-1=\zeta_{6}^{2}$. The 8th component of GetStructure(7) gives the index of the characters that are trivial on the above subgroups, its value is thus

[TABLE]

The vector $\Gamma_{s}(t)$ is given by (it is defined by (7))

[TABLE]

while

[TABLE]

Now that the players and the surrounding environment has been described, let us turn towards the main step of our proof: the recursion (32). We first check that

[TABLE]

Whence the relation

[TABLE]

with

[TABLE]

The call GetM1Inverse(7,GetStructure(7))^(-1) produces the matrix $M_{1}$. The matrices $N_{d}=dM_{1}^{-1}M_{d}$ are obtained by GetNds(7,GetStructure(7)). They are

[TABLE]

In order to check our script, we mention that the call

GetM1Inverse(7,GetStructure(7))^(-1)*GetNd(2,7,GetStructure(7))/2

gives $M_{2}$ for instance (and one can replace the parameter 2 that occurs twice with 3 or 6 to get $M_{3}$ and $M_{6}$). We have reached

[TABLE]

Our objective is $V_{s}(1)$ and we know how to compute $\Gamma_{s}(t)$ while, when $d$ is large, $V_{s}(dt)$ vanishes approximately; it is thus enough to iterate the above formula. We end the numerical example here.

6. Rational Euler Products

Let us recall the Witt decomposition. The readers will find in [6, Lemma 1] a result of the same flavour. We have simply modified the proof and setting as to accomodate polynomials having real numbers for coefficients.

Lemma 6.1**.**

Let $F(t)=1+a_{1}t+\ldots+a_{\delta}t^{\delta}\in\mathbb{R}[t]$ be a polynomial of degree $\delta$. Let $\alpha_{1},\ldots,\alpha_{\delta}$ be the inverses of its roots. Put $s_{F}(k)=\alpha_{1}^{k}+\ldots+\alpha_{\delta}^{k}$. The $s_{F}(k)$ are integers and satisfy the Newton-Girard recursion

[TABLE]

where we have defined $a_{\delta+1}=a_{\delta+2}=\ldots=0$. Put

[TABLE]

Let $\beta\geq 1$ be such that $\beta\geq\max_{j}|1/|\alpha_{j}|$. When $t$ belongs to any segment $\subset(-\beta,\beta)$, we have

[TABLE]

where the convergence is uniform in the given segment.

And how does the mathematician E. Witt enter the scene? In the paper [18] on Lie algebras, Witt produced in equation $(11)$ therein a decomposition that is the prototype of the above expansion.

Proof.

Since we follow the proof of [6, Lemma 1], we shall be rather sketchy. We write $F(t)=\prod_{i}(1-\alpha_{i}t)$. We thus have

[TABLE]

This series is absolutely convergent in any disc $|t|\leq b<1/\beta$ where $\beta=\max_{j}(1/|\alpha_{j}|)$. We may also decompose $tF^{\prime}(t)/F(t)$ in Lambert series as

[TABLE]

as some series shuffling in any disc of radius $b<\min(1,1/\beta)$ shows. The lemma follows readily by integrating the above relation. ∎

Lemma 6.2**.**

We use the hypotheses and notation of Lemma 6.1. Let $\beta\geq 2$ be larger than the inverse of the modulus of all the roots of $F(t)$. We have

[TABLE]

Proof.

We clearly have $|s_{F}(j)|\leq\deg F\cdot\beta^{j},$ so that

[TABLE]

∎

There are numerous easy upper estimates for the inverse of the modulus of all the roots of $F(t)$ in terms of its coefficients. Here is a simplistic one.

Lemma 6.3**.**

Let $F(X)=1+a_{1}X+\ldots+a_{\delta}X^{\delta}$ be a polynomial of degree $\delta$. Let $\rho$ be one of its roots. Show that, either $|\rho|\geq 1$ or $1/|\rho|\leq|a_{1}|+|a_{2}|+\ldots+|a_{\delta}|$.

Proof.

The readers may first notice that

[TABLE]

The conclusion is easy. ∎

Proof of Theorem 1.4.

The proof requires several steps. The very first one is a direct consequence of (49), which leads to the identity

[TABLE]

The absence of the $j=1$ term is due to our assumption that $(F(X)-G(X))/X^{2}\in\mathbb{Z}[X]$. Up to this point (50) is only established as a formal identity. Our second step is to establish (50) for all $t\in\mathbb{C}$ with $|t|<1/\beta$ and to control the rate of convergence. By Lemma 6.2, we know that $|b_{F}(j)-b_{G}(j)|\leq 2\max(\deg F,\deg G)\beta^{j}/j$. Therefore, for any bound $J$, we have

[TABLE]

as soon as $|t|<1/\beta$. We thus have

[TABLE]

where $|\log I_{1}|\leq 2\max(\deg F,\deg G)|t\beta|^{J+1}/[(1-|t\beta|)(J+1)]$.

Now that we have the expansion (52) for each prime $p$, we may combine them. We readily get

[TABLE]

where $I_{2}$ satisfies

[TABLE]

since $P\geq 2$ and $J\geq 3$. As announced earlier, we may rearrange the product over the primes $p$ and get

[TABLE]

The last step is to replace $\zeta_{P}(j;q,\mathcal{A})$ by the approximation, say $\zeta_{P}(j;q,\mathcal{A}|r)$ given by (13). We find that

[TABLE]

where $I_{3}$ satisfies

[TABLE]

with

[TABLE]

Therefore (and since $r\geq 2$)

[TABLE]

and this ends the proof. ∎

7. Counting the number of Lattice-Invariant Classes

It is of interest to count how many lattice-invariant classes there are, i.e. to determine the cardinality of $G^{\sharp}$ which is equally the number of cyclic subgroups, i.e. the cardinality of $\mathscr{G}$. We proceed in several steps.

Lemma 7.1**.**

Let $d\geq 1$ and $q\geq 1$ be two integers. The number $\rho(q;d)$ of solutions to the equation $x^{d}\equiv 1[q]$ is a multiplicative function of the variable $q$. When $p$ is a prime, we find that

[TABLE]

The function $d\mapsto\rho(q;d)$ is also multiplicative.

Proof.

The multiplicative character of $\rho(q;d)$ stems from the Chinese Remainder Theorem. In $\mathbb{Z}/p^{\alpha}\mathbb{Z}$ and $p\neq 2$, the equation $x^{d}\equiv 1[p^{\alpha}]$ has $(d,p^{\alpha-1}(p-1))$ as stated in [9, Corollary 2.42]; it is an easy consequence of the fact that $(\mathbb{Z}/p^{\alpha}\mathbb{Z})^{\times}$ is cyclic in this case.

When $p=2$, the equation $x^{d}\equiv 1[2]$ has exactly one solution, namely $x=1$. When $p=2$ and $\alpha\geq 2$, the multiplicative group $(\mathbb{Z}/p^{\alpha}\mathbb{Z})^{\times}$ is isomorphic to the direct product $(\mathbb{Z}/2\mathbb{Z},+)\times(\mathbb{Z}/2^{\alpha-2}\mathbb{Z},+)$. We find as a consequence of [4, Proposition 4.2.2] that $\rho(2^{\alpha},d)=2(d,2^{\alpha-2})$.

The multiplicativity of the function $d\mapsto\rho(q;d)$ folows from the explicit expression of $\rho(q;d)$: it is a product (over prime factors of $q$) of multiplicative functions of the variable $d$. ∎

Lemma 7.2**.**

The number $\rho^{*}(q;d)$ of elements of order $d$ in $(\mathbb{Z}/q\mathbb{Z})^{\times}$ is given by

[TABLE]

where $\rho(q;d)$ is defined and determined in Lemma 7.1.

Proof.

This is a consequence of the Moebius inversion formula as, by classifying the solution of $x^{d}\equiv 1[q]$ by their order, we find that $\rho(q;d)=\sum_{\ell|d}\rho^{*}(q;\ell)$. ∎

Proposition 7.3**.**

When $q\geq 3$, the number $|\mathscr{G}|$ of cyclic subgroups of $(\mathbb{Z}/q\mathbb{Z})^{\times}$ is given by

[TABLE]

where $\rho(q;d)$ is defined and determined in Lemma 7.1.

We have checked this expression with Sage via the function CardClassList of our script. The values have been checked against a direct count: we have the list of lattice-invariant classes, hence their number.

Proof.

Each cyclic subgroup of order $d$ has $\varphi(d)$ generators. Hence the number of cyclic subgroups of order $d$ is equal to $\rho^{*}(q;d)/\varphi(d)$, whence, by Lemma 7.2,

[TABLE]

To evaluate the inner sum, write $\varphi(q)=h_{1}h_{2}h_{3}$, where $h_{1}$ is the product of the $p^{v_{p}(\varphi(q))}$ with $p|\ell$ and $p|\varphi(q)/\ell$, then $h_{2}$ is the product of the $p^{v_{p}(\varphi(q))}$ with $p|\ell$ but $p\nmid\varphi(q)/\ell$ and $(h_{3},\ell)=1$ is what remains after division by $h_{1}h_{2}$. We readily find that

[TABLE]

This vanishes when $2|h_{3}$, so we can restrict our attention to even $\ell$’s. In which case we get

[TABLE]

We reverse to the variable $d$ rather than $\ell$ to write our lemma. We have also used the condition $q\geq 3$ to ensure that $2|\varphi(q)$. ∎

8. Notes on the implementation

The parameter $r$ is not very large, typically between 2 and 8. Since in (12), several products $d=d_{1}\cdots d_{v}$ are equal, we store the computed values of $\Gamma_{s}(dt)$ in the dictionary ComputedGammas in the function GetVs of the script LIEP.sage. We proceed similarly with the dictionary ComputedProductNdsM1Inverse for the products $N_{d_{1}}\cdots N_{d_{v}}M_{1}^{-1}$ in . Since the list $[d_{1},\cdots,d_{v}]$ cannot be a key for such a dictionary, we simply replace it by the tuple $(d_{1},\cdots,d_{v})$.

Concerning the general structure, the function GetStructure computes all the algebraical quantities that we need: the list of cyclic subgroups, the one of lattice-invariant classes, the exponent of our group, its character group, the set of invertible classes and, for each cyclic subgroup, the set of characters that are trivial on it.

Once the script is loaded via load(’LIEP.sage’), a typical call will be

GetVs(12, 2, 100, 300)

to compute modulo 12 the possible constants with $s=2$, asking for 100 decimal digits and using $P=300$. The output is self explanatory. The number of decimal digits asked for is roughly handled and one may lose precision in between, but this is indicated at the end (we observed no such phenomenon, but it may still happen!). A more precise treatment would first check the output and if the precision attained would not be enough, increase automatically this parameter. We prefer to let the users do that by themselves. The digits presented when WithLaTeX $=1$ are always accurate. Note that we expect the final result to be of size roughly unity, so we ask for is not the relative precision but the number of decimals. Hence, in the function GetGamma, we replace by an approximation of 0 the values that we know are insignificantly small. This is a true time-saver.

There are two subsequent optional parameters Verbose and WithLaTeX. The first one may take the values 0, 1 and 2; when equal to 0, the function will simply do its job and return the list of the invariant classes and the one of the computed lower and upper values. When equal to 1, its default value, some information on the computation is given. At level 2, more informations is given, but that should not concern the casual user. When the parameter Verbose is at least 1 and WithLaTeX is 1, the values of the constants will be further presented in a format suitable for inclusion in a LaTeX-file. For instance, the call

GetVs(12, 2, 100, 100, 1, 1)

is the one used to prepare this document.

To compute the Euler products as explained in Theorem 1.4, we have the function GetEulerProds(q, F, G, nbdecimals, bigP = 100, Verbose = 1, WithLaTeX = 0). Note that the parameter bigP may be increased during the run of the program to ensure that $P\geq 2\beta$ (a condition that is most of the time satisfied). We reused the same structure as the function GetVs, without calling it: this is to also keep all the precomputed datas. Since the coefficients $|b_{F}(j)-b_{G}(j)|$ may increase like $\beta^{j}$, we increase the working precision by $J\log\beta/\log 2$.

Checking

The values given here have been checked in several manners. The co-authors of this paper have computed several of the next values via independent scripts. We also provide the function GetVsChecker(q, s, borne = 10000) which computes approximate values of the same Euler products by simply truncating the Euler product representation. We checked with positive result the stability of our results with respect of the variation of the parameter $P$. This proved to be a very discriminating test.

Furthermore, approximate values for Shank’s and Lal’s constants are known (Finch in [1] gives 10 digits) and we agree on those. Finally, the web site [3] by X. Gourdon and P. Sebah is nowadays difficult to decypher but a postscript version is available on the same page. They give in section 4.4 the first fifty digits of the constant they call $A$ and which is

[TABLE]

Our result match the one of [3].

9. Some results

In this part, we exhibit some results for $s=2$ and small $q$’s. We decided to produce 100 decimal digits each time. Each computation took at most five seconds and we selected uniformly $P=100$.

Modulo 3

[TABLE]

[TABLE]

Modulo 4

[TABLE]

[TABLE]

Modulo 5

[TABLE]

[TABLE]

[TABLE]

Modulo 7

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Modulo 8

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Modulo 9

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Modulo 11

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Modulo 12

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Modulo 13

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Modulo 15

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Modulo 16

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Some notes on timing

We tried several large computations to get an idea of the limitations of our script, with the uniform choice $P=300$ and asking for 100 decimal digits. Since we did not run each computations hundred times to get an average timing, this table has to be taken with a pinch of salt. We present relative timing, knowing that the computation with $q=3$, $q=4$ or $q=4$ took about a tenth of a second.

[TABLE]

[TABLE]

This table shows that the value of $\varphi(q)$ is the main determinant of the time needed. The column with the tag “$\#d_{i}^{\prime}s$” contains the number of tuples $(d_{1},\cdots,d_{v})$ in the main formula.

Here is now a shorter table when asking 1000 decimal digits still with $P=300$. The time needed is still very decent.

[TABLE]

When asking for 5000 decimal digits and only $q=3$, it took about 16 minutes (with $P=500$) to get an answer, which essentially sets the horizon of the present method.