HIGH PRECISION COMPUTATION OF

HARDY-LITTLEWOOD CONSTANTS

Henri Cohen

Abstract. We explain how to compute slowly convergent Euler products or prime sums very efficiently, in particular the Hardy-Littlewood constants associated to prime producing polynomials.

In many contexts, it is necessary to compute slowly convergent Euler products or prime sums. Although it has been known for a long time that this can be done in non-trivial ways, it seems necessary to survey in detail the possible methods which can be used for doing this. In particular, we will explain how to compute to high precision the Hardy-Littlewood constants associated with prime producing polynomials.

We must first make precise what we mean by high precision computation of a constant. In the context of this paper, this means that we can compute the constant to any desired reasonable degree of accuracy (10, 100 or even 1000 decimal digits) in a reasonable amount of time. In particular, methods which apply to the computation of only 5 or 6 decimal digits are not acceptable here.

1. Computation of ζ(k)

As we shall see below, we will need to compute to reasonably high accuracy the ordinary Riemann ζ function at positive integers. We first explain how to do this. There are many possible methods for this task. Let us assume that we want to compute ζ(k) to an accuracy of 10^{−D} (since 1 < ζ(k) ≤ ζ(2) = π^{2}/6 = 1.644 … there is no real difference between absolute and relative accuracy here).

1.1. Using Euler-MacLaurin

Let N ≥ 1 be an integer. The Euler-McLaurin summation formula gives us immediately the expansion:

ζ(k) = ∑_{n=1}^{N} \frac{1}{n^{k}} + ∑_{i=0}^{p} \frac{B_{i}(k + i − 2)!}{(k − 1)!N^{k+i−1}} + R_{p}(k, N) ,

where B_{0} = 1, B_{1} = −1/2, B_{2} = 1/6, B_{3} = 0, B_{4} = −1/30, B_{5} = 0, B_{6} = 1/42, B_{7} = 0, B_{8} = −1/30, etc... are the Bernoulli numbers and |R_{p}(k, N)| is less than

Typeset by AMS-TEX

Key words and phrases. Hardy-Littlewood constant, prime sum, Euler product, prime numbers.

1991 AMS Classification 11Y40

the first non-zero term which is omitted, i.e. for p even

|R_{p}(N)| ≤ \frac{|B_{p+2}|(k + p)!}{(k − 1)!N^{k+p+1}} .

Since for p even we have

|B_{p+2}| = \frac{2ζ(p + 2)(p + 2)!}{(2π)^{p+2}} < \frac{π^{2}}{3} \frac{(p + 2)!}{(2π)^{p+2}} ,

we obtain an explicit upper bound for the error term |R_{p}(N)|. We can then easily find good values of N and p which give us the desired accuracy.

For example, if we want 50 decimal digits, we can choose...

The above will be enough for our purposes. We refer to [Coh-Oli] for more details, including computations of ζ(s) on the whole complex plane.

1.2. Using explicit formulas

We may of course also use the explicit formulas for ζ(k) when k is even. Recall that if k ≥ 2 is even, then

ζ(k) = \frac{|B_{k}|2^{k−1}π^{k}}{k!} .

This formula is extremely useful, but suffers from two handicaps. The first obvious one is that it applies only to k even. The second more subtle handicap is that it is not reasonable to apply it when k is large, since in that case B_{k} is a rational number with huge coefficients, and π^{k} also needs to be computed to high accuracy, although the final value of ζ(k) itself is extremely close to 1. In that case, it is much better to compute ζ(k) simply by computing a small partial sum.

It is not so well known that the first handicap of the explicit formula for ζ(k) can be partially lifted by the use of formulas initially due to S. Ramanujan. For example, he gives the following formula:

ζ(3) = \frac{7π^{3}}{180} − 2 ∑_{n≥1} \frac{1}{n^{3}(e^{2πn} − 1)} .

The general formulas (for k > 1) are as follows.

Theorem 1.1. If k ≡ 3 (mod 4), we have

ζ(k) = \frac{2^{k−1}π^{k}}{(k + 1)!} ∑_{0≤n≤(k+1)/2}(−1)^{n−1}\binom{k + 1}{2n}B_{k+1−2n}B_{2n} − 2 ∑_{n≥1} \frac{1}{n^{k}(e^{2πn} − 1)} .

If k ≡ 1 (mod 4), k ≥ 5, we have

\begin{aligned}ζ(k) & = \frac{(2π)^{k}}{(k + 1)!(k − 1)} ∑_{0≤n≤(k−1)/4}(−1)^{n}(k + 1 − 4n)\binom{k + 1}{2n}B_{k+1−2n}B_{2n} \\ & − 2 ∑_{n≥1} \frac{e^{2πn}(1 + (4πn)/(k − 1)) − 1}{n^{k}(e^{2πn} − 1)^{2}} . \end{aligned}

These formulas are proved by computing the period functions associated to the Eisenstein series of weight k + 1 on PSL_{2}(ℤ).

Note that e^{2π} > 535, hence the convergence of these series, although linear, is rather fast (almost 3 decimal digits per term).

2. Computation of prime sums and simple Euler products

Here we give some examples of sums which can easily be computed from the values of the Riemann ζ function at positive integers.

2.1. Computation of ∑_{p} p^{−m}

As a first example, we explain how to compute

S_{m} = ∑_{p} \frac{1}{p^{m}}

for m ≥ 2 integral, where here and in the rest of this paper the sum over p (likewise products over p) is implicitly over the positive prime numbers.

The link between prime sums and the Riemann ζ function is evidently given by the Euler product for ζ, which can be written

log(ζ(m)) = − ∑_{p} log(1 − p^{−m}) = ∑_{p} ∑_{k≥1} \frac{1}{kp^{km}} = ∑_{k≥1} \frac{S_{km}}{k} .

By Möbius inversion, we obtain

S_{m} = ∑_{k≥1} \frac{μ(k)}{k} log(ζ(km)) .

Note that log(ζ(km)) = O(2^{−km}) so the above series for S_{m} already converges rapidly enough to enable us to compute it to a few decimals. It is however much better to compute the partial sum for the series S_{m} for primes up to A, and use the ζ expansion only for the remainder term. More precisely, set

ζ_{p>A}(s) = ζ(s) ∏_{p≤A} (1 − \frac{1}{p^{s}} ) .

Then the above formula for S_{m} is easily generalized to

S_{m} = ∑_{p≤A} \frac{1}{p^{m}} + ∑_{k≥1} \frac{μ(k)}{k} log(ζ_{p>A}(km)) .

Since log(ζ_{p>A}(km)) = O(A^{−km}), this converges much faster than the above series for S_{m} if A is not taken too small. The optimal value of A depends on the number of desired decimal digits of the result and is easily computed (values of 20 to 100 are reasonable if many digits are desired), and in a few seconds, we can thus compute several hundred decimal digits of S_{m}. For example, we have

S_{2} = 0.45224742004106549850654336483224793417323134323989 …

S_{3} = 0.17476263929944353642311331466570670097541212192614 …

Note that the computation time is less than 0.02 seconds.

2.2. Computation of ∑_{p} f(p)

More generally, assume that we have a function f(p) which can be expanded into an absolutely convergent series

f(p) = ∑_{m≥1} \frac{u(m)}{p^{m}} ,

and we want to compute S_{f} = ∑_{p} f(p) (a necessary condition for convergence is of course u(1) = 0). We could write S_{f} = ∑_{m≥1} u(m)S_{m} and compute S_{m} as explained in the preceding section, but this would be very wasteful since we would need to recompute many times the same quantities. It is much better to first obtain a specific formula as follows:

S_{f} = ∑_{m≥1} u(m)S_{m} = ∑_{m≥1} u(m) ∑_{k≥1} \frac{μ(k)}{k} log(ζ(km))

= ∑_{N≥1} log(ζ(N)) ∑_{k|N} \frac{μ(k)}{k}u(N/k)

= ∑_{N≥1} \frac{u^{∗}(N)}{N} log(ζ(N)) ,

where

u^{∗}(N) = ∑_{k|N} μ(k)(N/k)u(N/k)

is the multiplicative convolution of the arithmetic functions μ(n) and n u(n). As usual, this is written

S_{f} = ∑_{p≤A} f(p) + ∑_{N≥1} \frac{u^{∗}(N)}{N} log(ζ_{p>A}(N)) .

In this way, we compute

∑_{p} \frac{1}{p(p − 1)} = 0.77315666904979512786436745985594239561874133608318 …

(u(m) = 1 for m ≥ 2, u(1) = 0, hence u^{∗}(N) = ϕ(N) − μ(N)),

∑_{p} \frac{1}{(p − 1)^{2}} = 1.3750649947486352879172531305224396991795999601753 …

(u(m) = m − 1, hence u^{∗}(N) = N^{2} ∏_{p|N}(1 − 1/p^{2}) − ϕ(N)) and similar sums.

By taking derivatives or limits, in an analogous manner we easily compute

∑_{p} \frac{log p}{p^{2}} = 0.4930911093687644621978262050564912580555880596367 …

lim_{s→1^{+}} ∑_{p} \frac{1}{p^{s}} −log(ζ(s)) = −0.3157184520538900768510852514737065719905926876787 …

lim_{x→∞} ∑_{p≤x} \frac{1}{p}−log log x = 0.2614972128476427837554268386086958590515666482612 …

2.3. Computation of ∑_{p} 1/(p^{m} log p)

We can also use the same method to compute the series ∑_{p} 1/(p^{m} log p) which converge for m ≥ 1. If we set

I(m) = ∫_{m}^{∞} log ζ(t) dt ,

as above, Möbius inversion shows that

∑_{p} \frac{1}{p^{m} log p} = ∑_{k≥1} \frac{μ(k)}{k^{2}}I(km) .

As usual, it is preferable to use ζ_{p>A}(s) instead of ζ(s). We note that

∫_{m}^{∞} − log(1 − p^{−t}) dt = ∑_{l≥1} \frac{1}{l} ∫_{m}^{∞} p^{−tl} dt = \frac{1}{log p} ∑_{l≥1} \frac{1}{l^{2} p^{ml}} = \frac{1}{log p} Li_{2}(1/p^{m})

where

Li_{2}(x) = ∑_{l≥1} \frac{x^{l}}{l^{2}}

is the dilogarithm function. Thus

∑_{p} \frac{1}{p^{m} log p} = ∑_{p≤A} \frac{1}{p^{m} log p} + ∑_{k≥1} \frac{μ(k)}{k^{2}} (I(km) − ∑_{p≤A} \frac{Li_{2}(p^{−km})}{log p}) .

Thus, the whole problem reduces to the computation of I(m) for integral values of m.

For this, we use once again the Euler-MacLaurin summation formula in the following form. Let N ≥ 1 be an integer and f be a “nice” function (we do not make this precise yet). Then

∫_{m}^{∞} f(t) dt = \frac{1}{N} ∑_{n>mN} f (\frac{n}{N} ) + \frac{f(m)}{2N} + ∑_{i=2}^{p} \frac{B_{i}}{N^{i}i!}f^{(i−1)}(m) + R_{p}(N) ,

where R_{p}(N) is smaller in absolute value than the first neglected term.

If we apply this to f(t) = log(ζ_{p>A}(t)) we can thus reduce the computation of I(m) to the computation of some derivatives of f(t). To compute these derivatives, we approximate f^{(n)}(t) by the formula

f^{(n)}(t) = \frac{∑_{0≤k≤n}(−1)^{k}\binom{n}{k}f(x + (n/2 − k) ∗ ε)}{ε^{n}} + O(ε^{2})

where ε is chosen small. In practice, if we want D decimal digits in the result, we choose ε = 10^{−D/2} and compute the numerator with (n + 1)D/2 decimal digits.

The above method is not valid for m = 1, since f(t) and its derivatives are not defined at t = 1. Instead, we apply it to the function

f_{2}(t) = ∑_{n≥1}(−1)^{n−1}n^{−t} = (1 − 2^{1−t})ζ(t)

which is defined for t > 0 (and can be computed accurately using a modification of Euler-Maclaurin). It is easy to check that

I(1) = \frac{π^{2}}{6 log 2} + ∫_{1}^{∞} f_{2}(t) dt ,

which enables us to compute I(1).

We find for example

I(1) = 1.79756995862873940793025078212153165864605183075709 …

I(2) = 0.536526945921177100961719019548794400667004728755026 …

∑_{p} \frac{1}{p log p} = 1.63661632335126086856965800392186367118159707613129 …

∑_{p} \frac{1}{p^{2} log p} = 0.5077821878591993187743751037947055704669736717043207 …

It is quite remarkable that one is able to compute ∑_{p} 1/(p log p) to any desired degree of accuracy, since not only is this sum over the (apparently irregular) set of primes as the preceding sums that we have considered, but its convergence is extremely slow (in 1/ log(X) where X is the upper bound for the primes). Thus, it would be inconceivable to compute it using the naive method since even with a table of primes up to 10^{20} (already an impossible practical limit), we would obtain less than 2 decimal digits.

2.4. Computation of simple Euler products

The word “simple” refers here to Euler products of the form ∏_{p} g(p) where g(p) can be expanded into an absolutely convergent series

g(p) = 1 + ∑_{m≥1} \frac{v(m)}{p^{m}} ,

where once again a necessary condition for convergence is v(1) = 0. We can immediately reduce to the case treated in the preceding section by taking logarithms. Thus, we set

f(p) = log(g(p)) = ∑_{m≥1} \frac{u(m)}{p^{m}} ,

and the sequence u(m) can easily be computed by induction from the sequence v(m). More precisely, setting F (x) = f(1/x), G(x) = g(1/x) and taking the derivative with respect to the formal variable x, we have

f′(x) = g′(x)/g(x) = ∑_{n≥0}(n + 1)u(n + 1)x^{n}

hence setting v(0) = 1,

( ∑_{m≥0} v(m)x^{m}) (∑_{n≥0}(n + 1)u(n + 1)x^{n}) = ∑_{N≥0}(N + 1)v(N + 1)x^{N} ,

which gives by identification the recursion

v(N + 1) = \frac{1}{N + 1} ∑_{0≤n≤N}(n + 1)u(n + 1)v(N − n) .

We can then use the formulas explained in the preceding section to compute S_{f}, and deduce our Euler product as e^{S_{f}}.

We give two examples. The first one is the computation of the twin prime and Goldbach constant

C_{g} = ∏_{p>2} \frac{p(p − 2)}{(p − 1)^{2}} .

For the sake of completeness, we recall that conjecturally the number of twin primes (q, q + 2) with q ≤ X is asymptotic to

2C_{g} \frac{X}{log^{2} X} ,

and that the number of decompositions of an even number N as a sum of two primes is asymptotic to

2C_{g} \frac{N}{log^{2} N} ∏_{p|N, p>2} \frac{p − 1}{p − 2} .

To compute C_{g}, we proceed as explained above. We have

log (\frac{p(p − 2)}{(p − 1)^{2}}) = − ∑_{m≥1} \frac{2^{m} − 2}{mp^{m}} ,

hence with the above notations, we have u(m) = (2^{m} − 2)/m. Thus,

u^{∗}(N) = − ∑_{k|N} μ(k)(2^{N/k} − 2) ,

and in this way we compute

C_{g} = 0.6601618158468695739278121100145557784326233602847 …

As a second example, we can compute the Artin constant

C_{a} = ∏_{p} (1 − \frac{1}{p(p − 1)}) .

Conjecturally, the number of primes less than or equal to x for which 2 is a primitive root modulo p should have density C_{a} in the set of primes.

We have

u^{∗}(N) = − ∑_{k|N} μ(k) ((\frac{1 + \sqrt{5}}{2})^{N/k} + (\frac{1 − \sqrt{5}}{2})^{N/k} − 1) ,

and in this way we compute

C_{a} = 0.373955813619202288054728054346416415111629248606 …

3. Computation of L(k, χ) for quadratic characters χ

Let χ be a Dirichlet character modulo f. Recall that the L-function associated to χ is

L(s, χ) = ∑_{n≥1} \frac{χ(n)}{n^{s}} = ∏_{p} (1 − \frac{χ(p)}{p^{s}})^{−1} .

In the next section, we will need to compute L(k, χ) for k positive integral. We explain here how this is done. We assume f to be sufficiently large for the naive methods to fail, but not too large that the method explained below becomes too slow. In practice f < 10^{14} is OK.

3.1. A general formula for Dirichlet series with functional equation The computation of L(s, χ) is based on the following theorem, whose proof is

not really difficult. We state in a much more general setting than needed here.

Theorem 3.1. For i = 1 and i = 2, let L_{i}(s) = ∑_{n≥1} a_{i}(n)n^{−s} be a Dirichlet series such that the a_{i}(n) have polynomial growth. For i = 1 and i = 2, let γ_{i}(s) be functions of the form A^{s} ∏_{i} Γ(a_{i}s + b_{i}) with A > 0, a_{i} positive rational numbers and b_{i} complex numbers, such that the following conditions hold.

(1) The functions Λ_{i}(s) = γ_{i}(s)L_{i}(s) extend analytically to the whole complex plane into meromorphic functions which are bounded on vertical strips.

(2) There exist rational functions ϕ_{i}(s) such that Λ_{i}(s) − ϕ_{i}(s) is an entire function, in other words Λ_{i}(s) has only a finite number of poles.

(3) There exists a functional equation

Λ_{1}(k − s) = w · Λ_{2}(s)

for some constant w ∈ ℂ^{∗}, some real number k, valid for all s for which both sides are defined.

Let W_{i}(t) be the inverse Mellin transform of γ_{i}(s), in other words such that

γ_{i}(s) = ∫_{0}^{∞} t^{s}W_{i}(t) \frac{dt}{t} ,

so that for any δ sufficiently large,

W_{i}(t) = \frac{1}{2iπ} ∫_{δ−i∞}^{δ+i∞} t^{−s}γ_{i}(s) ds .

Define the functions F_{i}(s, x) by

F_{i}(s, x) = ∫_{x}^{∞} t^{s}W_{i}(t) \frac{dt}{t} .

and the functions p_{i}(s, x) by

p_{i}(s, x) = ∑_{a} ∫_{x}^{∞} t^{−s} Res_{a}(t^{z}ϕ_{i}(z)) \frac{dt}{t}

where the sum is over all poles a of ϕ_{i}(z) (or equivalently of Λ_{i}(z)).

Then for all t_{0} > 0, we have

Λ_{1}(s) = ∑_{n≥1} \frac{a_{1}(n)}{n^{s}}F_{1}(s, nt_{0}) + w ∑_{n≥1} \frac{a_{2}(n)}{n^{k−s}}F_{2}(k − s, n/t_{0}) + p_{1}(s, 1/t_{0}) .

and symmetrically

Λ_{2}(s) = ∑_{n≥1} \frac{a_{2}(n)}{n^{s}}F_{2}(s, n/t_{0}) + w^{−1} ∑_{n≥1} \frac{a_{1}(n)}{n^{k−s}}F_{1}(k − s, nt_{0}) + p_{2}(s, t_{0}) .

If in addition t_{0} = 1 and ϕ_{i}(s) tends to 0 when |s| tends to infinity, we have p_{i}(s, 1) = ϕ_{i}(s).

Although the statement of this theorem seems technical, its proof is easy. It can be roughly summarized by saying that once we has proved a functional equation, we can compute rapidly the Dirichlet series involved.

Let us consider the special case of Dirichlet characters. If χ is a Dirichlet character modulo d of conductor f | d, then χ = χ_{0}χ_{f} where χ_{0} is the trivial character modulo d and χ_{f} is the primitive character modulo f equivalent to χ. Clearly,

L(s, χ) = ∏_{p∤d} (1 − \frac{χ_{f}(p)}{p^{s}})^{−1} = L(s, χ_{f}) ∏_{p|d, p∤f} (1 − \frac{χ_{f}(p)}{p^{s}}) ,

so to compute L(s, χ) it is sufficient to compute L(s, χ_{f}), and hence we may assume that χ is a primitive character.

Proposition 3.2. Let χ be a primitive Dirichlet character modulo f > 1. Let

Γ(s, x) = ∫_{x}^{∞} e^{−t}t^{s} \frac{dt}{t}

be the incomplete Γ function, and let

G(χ) = ∑_{r mod f} χ(r)e^{2iπr/f}

be the Gauss sum associated to f. Then, if χ is an even character

Γ (\frac{s}{2}) L(s, χ) = ∑_{n≥1} \frac{χ(n)}{n^{s}}Γ (\frac{s}{2}, \frac{πn^{2}}{f} )+G(χ)π^{s−1/2}f^{−s} ∑_{n≥1} \frac{\overline{χ(n)}}{n^{1−s}} Γ (\frac{1 − s}{2}, \frac{πn^{2}}{f} ) ,

while if χ is an odd character

Γ (\frac{s + 1}{2}) L(s, χ) = ∑_{n≥1} \frac{χ(n)}{n^{s}}Γ (\frac{s + 1}{2}, \frac{πn^{2}}{f} )+\frac{G(χ)}{i}π^{s−1/2}f^{−s} ∑_{n≥1} \frac{\overline{χ(n)}}{n^{1−s}} Γ (\frac{2 − s}{2}, \frac{πn^{2}}{f} ) .

In the special case where χ(n) = (\frac{D}{n} ) is a quadratic character modulo D and

D is a fundamental discriminant (so χ is a primitive character), then \overline{χ} = χ, χ is even or odd according to whether D > 0 or D < 0, f = |D|, and G(χ) = \sqrt{D}, i.e. G(χ) = \sqrt{D} if D > 0 and G(χ) = i\sqrt{|D|} if D < 0.

The above series for L(s, χ) converges exponentially, but we need to take approximately O(f^{1/2}) terms before the terms start to become small. Thus, the fundamental limit to the above method is the size of f. As mentioned, f < 10^{14} can be done, although not without difficulty when f is close to this upper limit, but it is not reasonable to use the above method when f is much larger.

3.2. Computation of the incomplete Gamma function

In the general case, we will need to compute the functions F_{i}(s, x). We consider here the simplest case corresponding to Dirichlet characters, i.e. the case where F_{i}(s, x) is an incomplete Γ function.

Recall that

Γ(s, x) = ∫_{x}^{∞} e^{−t}t^{s} \frac{dt}{t} .

To compute Γ(s, x), we proceed as follows. If x is small (say x ≤ 10), we write

Γ(s, x) = Γ(s) − ∫_{0}^{x} e^{−t}t^{s−1} dt = Γ(s) − ∑_{n≥0} \frac{(−1)^{n}x^{n+s}}{n!(n + s)} .

This series converges for every x, but cannot be used for large values of x because of catastrophic cancellation, the same that would occur if we tried to compute e^{−t} for t large by using the power series.

If x is large, we can use instead a classical continued fraction expansion

Γ(s, x) = \frac{e^{−x}x^{s}}{x + 1 − s − \frac{1 · (1 − s)}{x + 3 − s − \frac{2 · (2 − s)}{x + 5 − s − ⋯}}}

Finally, in our case we will only need the values of the L series for positive integral values of s, hence we will need the values of the incomplete gamma function only for integral or half integral values of s, and a simple recursion then reduces the computation of these values to the computation of the error function (essentially Γ(1/2, x)) and of the exponential integral (essentially Γ(0, x)). These computations are made using the same principles, but the formulas can be somewhat simplified and the continued fraction accelerated (see for example [Coh, Chapter 5]).

4. Computation of Hardy-Littlewood constants of quadratic polynomials

4.1. Computation of prime sums and Euler products with quadratic character

We now would like to compute prime sums and Euler products which, in addition

to involving regular functions over primes, also involve a quadratic character (\frac{D}{n} )

where D is a not necessarily fundamental discriminant.

A natural idea is to introduce, in addition to the Riemann ζ function, the Dirich-

let L-functions L(s, χ) for Dirichlet characters χ, in particular for χ(n) = (\frac{D}{n} ).

We can easily generalize the results of Section 2.1 to this case. If for any character χ we set

S_{m}(χ) = ∑_{p} \frac{χ(p)}{p^{m}}

then Möbius inversion gives us

S_{m}(χ) = ∑_{k≥1} \frac{μ(k)}{k} log(L(km, χ^{k})) ,

which as usual in practice is computed as

S_{m}(χ) = ∑_{p≤A} \frac{χ(p)}{p^{m}} + ∑_{k≥1} \frac{μ(k)}{k} log(L_{p>A}(km, χ^{k}))

with evident notations. Using the method of the preceding section to compute L(N, ψ), we may therefore easily compute the sums S_{m}(χ).

Let us specialize to the case where χ is a quadratic character (\frac{D}{n} ) for some (not

necessarily fundamental) discriminant D. Then χ^{k} = χ if k is odd and χ^{k} = χ_{0}, the trivial character modulo D (χ_{0}(n) = 1 if (n, D) = 1 and χ_{0}(n) = 0 otherwise) if k ≥ 2 is even.

Thus,

S_{m}(χ) = ∑_{χ(p)=1} \frac{1}{p^{m}} − ∑_{χ(p)=−1} \frac{1}{p^{m}} ,

hence

∑_{χ(p)=1} = \frac{1}{2}(S_{m}(χ_{0}) + S_{m}(χ))

and

∑_{χ(p)=−1} = \frac{1}{2}(S_{m}(χ_{0}) − S_{m}(χ)) ,

where

S_{m}(χ_{0}) = S_{m} − ∑_{p|D} \frac{1}{p^{m}} .

A little computation now shows that if

f_{j}(p) = ∑_{m≥1} \frac{u_{j}(m)}{p^{m}}

for j = 1, j = 0, j = −1 corresponding to the 3 possible values of χ(p), then

S_{f} = ∑_{χ(p)=1} f_{1}(p) + ∑_{χ(p)=0} f_{0}(p) + ∑_{χ(p)=−1} f_{−1}(p)

is given by the following formula.

For any arithmetic function u(m) and integer r ≥ 1, set

u^{∗,r}(N) = ∑_{k|N, (r,k)=1} μ(k)(N/k)u(N/k) ,

so that for example u^{∗,1}(N) = u^{∗}(N) with the notation of Section 2.2. Also, set

ζ_{p∤D}(s) = L(s, χ_{0}) = ζ(s) ∏_{p|D} (1 − \frac{1}{p^{s}} ) .

Then

S_{f} = \frac{1}{2} ∑_{N≥1} \frac{log(L(N, χ))}{N} (u_{1}^{∗,2}(N) − u_{−1}^{∗,2}(N))

+ \frac{1}{2} ∑_{N≥1} \frac{log(ζ_{p∤D}(N))}{N} (u_{1}^{∗,2}(N) + u_{−1}^{∗,2}(N) − 2u_{1}^{∗,2}(N/2))

+ ∑_{p|D} f_{0}(p) ,

where here and elsewhere, for any arithmetic function u, u(x) = 0 if x is not a postive integer.

In practice, as usual, we use L_{p>A}(N, χ) and ζ_{p>A, p∤D}(N) to obtain faster convergence.

It is clear that such formulas can be generalized to Dirichlet characters χ of higher order r, and the result involves in a simple manner the functions L(N, χ^{b}) for 0 ≤ b < r. The quadratic case is especially simple since it involves L(N, χ) and L(N, χ^{0}) = L(N, χ_{0}) which is essentially the Riemann ζ function.

4.2. Computation of Hardy-Littlewood constants of quadratic polynomials

The main application that we have in mind of the results of Section 4.1 is the computation of the Hardy-Littlewood constants of quadratic polynomials defined as follows.

Let A(X) be a polynomial with integer coefficients. Assume that A(X) is irreducible in ℚ[X] and has content 1. For any prime p, let ω(p) = ω(p, A) be the number of solutions in 𝔽_{p} of A(x) = 0. Then conjecturally, the number of integers n ≤ N such that A(n) is prime should be asymptotic to

\frac{H(A)}{deg(A)} \frac{N}{log N} ,

where

H(A) = ∏_{p} \frac{1 − \frac{ω(p)}{p}}{1 − \frac{1}{p}} .

The constant H(A) is the Hardy-Littlewood constant of the polynomial A. It is therefore interesting to compute this Euler product. Since the case of linear polynomials is trivial, in the present section we treat the next simplest case of quadratic polynomials. We will consider polynomials of larger degree later. Thus, let A(X) = aX^{2} + bX + c be an irreducible quadratic polynomial with (a, b, c) = 1, and let D = b^{2} − 4ac be its discriminant. It is easy to see that ω(p) is given by the following formulas.

• If p ∤ a, then ω(p) = 1 + (\frac{D}{p} ), including for p = 2.

• If p | a and p ∤ b, then ω(p) = 1.

• If p | a and p | b (hence p ∤ c), then ω(p) = 0.

A little computation shows that this implies that the Hardy-Littlewood constant H(A) is given by

H(A) = c_{2} ∏_{p>2} (1 − \frac{(\frac{D}{p} )}{p − 1}) ∏_{p|a,p∤2b} \frac{p − 1}{p − 2} ∏_{p|(a,b),p>2} \frac{p}{p − 1} ,

where c_{2} = 0 if c is even and a + b is even, and if c is odd, c_{2} = 1/2 if a + b is odd, c_{2} = 1 if a + b is even.

Write D = D_{0}f^{2} with D_{0} a fundamental discriminant. Since

∏_{p>2} (1 − \frac{(\frac{D}{p} )}{p − 1}) = ∏_{p>2} (1 − \frac{(\frac{D_{0}}{p} )}{p − 1}) ∏_{p|D,p>2} (1 − \frac{(\frac{D_{0}}{p} )}{p − 1})^{−1} .

Thus, we must compute the Euler product

C(D_{0}) = ∏_{p>2} (1 − \frac{(\frac{D_{0}}{p} )}{p − 1})

for a fundamental discriminant D_{0}.

As in Section 4.1, set χ(n) = (\frac{D_{0}}{n} ). Then

− log(C(D_{0})) = − ∑_{p,χ(p)=1} log \frac{p − 2}{p − 1} − ∑_{p,χ(p)=−1} log \frac{p}{p − 1} ,

hence with the notation of Section 4.1, C(D_{0}) = e^{−S_{f}} with f_{0}(p) = 0,

f_{1}(p) = − log (\frac{p − 2}{p − 1}) = ∑_{m≥1} \frac{2^{m} − 1}{mp^{m}} ,

f_{−1}(p) = − log (\frac{p}{p − 1}) = − ∑_{m≥1} \frac{1}{mp^{m}} ,

so that m u_{1}(m) = 2^{m} − 1 and m u_{−1}(m) = −1. Thus,

u_{1}^{∗,2}(N) = ∑_{k|N 2∤k} μ(k)(2^{N/k} − 1)

and

u_{−1}^{∗,2}(N) = − ∑_{k|N 2∤k} μ(k) .

Set

a(N) = \frac{u_{1}^{∗,2}(N) − u_{−1}^{∗,2}(N)}{2N}

and

b(N) = \frac{u_{1}^{∗,2}(N) + u_{−1}^{∗,2}(N) − 2u_{1}^{∗,2}(N/2)}{2N} .

so that

S_{f} = ∑_{N≥1} a(N) log(L(N, χ)) + ∑_{N≥1} b(N) log(ζ_{p∤D_{0}}(N)) .

Then by the above we have

a(N) = \frac{1}{2N} ∑_{k|N 2∤k} μ(k)2^{N/k}

and b(N) = a(N) − a(N/2) if N > 1, while b(1) = 0 (recall that we set a(x) = 0 if x is not integral). We can thus compute C(D_{0}), using in practice L_{p>A}(N, χ) and ζ_{p>A,p∤D_{0}}(N).

For example, in this way we compute

H(X^{2} + X + 41) = 3.319773177471421665323556857649887966468554585653 …

H(X^{2} + X + 75) = 0.310976679925987170004356287429628414529121902600 …

The main limitation of this method is our ability to compute L(N, χ) for reasonable values of N (typically N ≤ 30) to the required degree of accuracy. As mentioned in the preceding section, although the method explained there using the functional equation is efficient, it is a O(D_{0}^{1/2+ϵ}) method, and when D_{0} becomes large, it is quite a difficult computation. We have computed (to 20 decimal digits, say) the Hardy-Littlewood constants of quadratic polynomials having discriminants of the order of 10^{13} or 10^{14}, but this has required several hours if not days of CPU time. If D_{0} is much larger (say 10^{25}), I do not see any reasonable method to compute C(D_{0}) apart from computing the Euler product directly without using any tricks.

5. Computation of Hardy-Littlewood constants of general polynomials

In this section, we generalize the methods of Section 4 to the case of general polynomials in ℤ[X], of course assumed to be irreducible and content equal to 1.

For this, we need in addition to L-functions of Dirichlet characters, the Artin L-functions of the irreducible representations of the Galois group of the polynomial. Although the results will be stated in terms of these L-functions, I emphasize that we do not need the validity of the unproved Artin’s conjecture on the holomorphy of the L-functions for the results to be valid, since we may reduce to the computation of Dedekind ζ functions. We will however need this conjecture (which, in the cases we need, is proved) for the explicit computation of the necessary L-functions.

5.1. The case of cyclic cubic polynomials

We first look at the case of cubic polynomials, and will mention briefly in the next section how this (easily) generalizes to higher degrees. A cubic polynomial A(X) can have 2 possible Galois groups: the cyclic group C_{3} when disc(A) is a square, and the symmetric group S_{3} otherwise.

Let A(X) be an irreducible content-free polynomial of degree 3 in ℤ[X] with cyclic Galois group. Then, apart from a small finite number of primes (those dividing the leading term of A(X) and those dividing the discriminant D of the polynomial A), the quantity ω(p) can immediately be read from the splitting of p in the number field K defined by A over ℚ. For simplicity, we assume A monic (otherwise, we would also have to exclude primes dividing the leading term).

Consider first the cyclic case. A prime can then either be inert or totally split. Thus, excluding implicitly all primes dividing D in the sums and products, we have

ζ_{K,p∤D}(s) = ∏_{ω(p)=0} (1 − \frac{1}{p^{3s}} )^{−1} ∏_{ω(p)=3} (1 − \frac{1}{p^{s}} )^{−3} .

A slightly modified Möbius inversion gives

∑_{ω(p)=0} \frac{1}{p^{m}} = −\frac{1}{3} ∑_{3∤k} \frac{μ(k)}{k} log(ζ_{K,p∤D}(km)) + ∑_{3∤k} \frac{μ(k)}{k} log(ζ_{p∤D}(km))

and

∑_{ω(p)=3} \frac{1}{p^{m}} = \frac{1}{3} ∑_{3∤k} \frac{μ(k)}{k} log(ζ_{K,p∤D}(km)) + ∑_{3|k} \frac{μ(k)}{k} log(ζ_{p∤D}(km)) .

This allows us to compute sums involving ω(p), and in particular the Hardy-Littlewood constant. This is equivalent to computing sums involving the values of the characters of the group C_{3}, which are all of degree 1.

If we set

H_{D}(A) = ∏_{p∤D} \frac{1 − \frac{ω(p)}{p}}{1 − \frac{1}{p}}

(which differs only by an easily computed finite Euler product from the Hardy-Littlewood constant H(A)), then we obtain

− log(H_{D}(A)) = log lim_{s→1^{+}} (\frac{ζ_{K,p∤D}(s)}{ζ_{p∤D}(s)}) + ∑_{N≥2} a(N) log(ζ_{K,p∤D}(N))

− ∑_{N≥2} a(N/3) log(ζ_{p∤D}(N))

with

a(N) = \frac{1}{3N} ∑_{k|N,3∤k} μ(k)3^{N/k} .

If we let L_{K}(s) = ζ_{K}(s)/ζ(s), which is in fact a product of two L-series of Hecke characters, this can also be written

− log(H_{D}(A)) = ∑_{N≥1} a(N) log(L_{K,p∤D}(N)) + ∑_{N≥1} b(N) log(ζ_{p∤D}(N))

with a(N) as above and b(N) = a(N) − a(N/3) if N > 1, while b(1) = 0. Not surprisingly, this is completely analogous to the formula obtained in the quadratic case.

To compute the function ζ_{K}(N) we use Theorem 3.1 applied to the Hecke Lseries whose product is equal to L_{K}. More precisely, set χ(p) = 1 if ω(p) = 3, χ(p) = e^{2iπ/3} if ω(p) = 0 and χ(p) = 0 if p | D. Then clearly L_{K}(s) = L(s, χ)L(s, \overline{χ}) and on the other hand L(s, χ) can be considered as a Dirichlet L-function with (by abuse of notation) χ a character of order 3 modulo D. The advantage is that these Dirichlet L-functions have a single gamma factor Γ(s/2) in their functional equation (since K is totally real), hence the functions F_{i} which occur are again incomplete gamma functions. We leave the (easy) details to the reader.

5.2. The case of non-cyclic cubic polynomials

We now assume that the Galois group of the polynomial A is non-cyclic, i.e. is isomorphic to S_{3}. Once again, up to a finite number of primes, we need only to look at the splitting of primes in the number field K. In addition to the cases occuring in the cyclic case (ω(p) = 0 if p is inert, ω(p) = 3 if p is totally split), we also have the case where ω(p) = 1 if p = 𝔭_{1}𝔭_{2} where 𝔭_{i} is a prime ideal of degree i. Since we have 3 cases, we need 3 functions. As in the cyclic case, we can either use ζ functions or L-functions. If we use ζ functions, we take ζ_{K}, ζ and ζ_{k}, the Dedekind ζ function of the quadratic field ℚ(\sqrt{D}), or in other words, of the unique quadratic subfield of the Galois closure of K. If we use L-functions, we take the functions ζ, L_{D} = ζ_{k}/ζ and L_{K} = ζ_{K}/ζ which is the Artin L-function of the irreducible character of degree 2 of S_{3}.

Proceeding as before, and noting that (\frac{D}{p} ) = −1 if and only if ω(p) = 1, and

as usual omitting primes dividing D, we have

ζ_{K,p∤D}(s) = ∏_{ω(p)=0} (1 − \frac{1}{p^{3s}} )^{−1} ∏_{ω(p)=1} (1 − \frac{1}{p^{2s}} )^{−1} (1 − \frac{1}{p^{s}} )^{−1} ∏_{ω(p)=3} (1 − \frac{1}{p^{s}} )^{−3} ,

ζ_{k,p∤D}(s) = ∏_{ω(p)=0} (1 − \frac{1}{p^{s}} )^{−2} ∏_{ω(p)=1} (1 − \frac{1}{p^{2s}} )^{−1} ∏_{ω(p)=3} (1 − \frac{1}{p^{s}} )^{−2} .

A slightly modified Möbius inversion gives

∑_{ω(p)=0} \frac{1}{p^{m}} = \frac{1}{3} ∑_{3∤N} \frac{μ(N)}{N} log (\frac{ζ_{k}ζ}{ζ_{K}} _{p∤D}(Nm))

∑_{ω(p)=1} \frac{1}{p^{m}} = \frac{1}{2} ∑_{2∤N} \frac{μ(N)}{N} log (\frac{ζ^{2}}{ζ_{k}} _{p∤D}(Nm))

∑_{ω(p)=3} \frac{1}{p^{m}} = ∑ \frac{μ(N)}{N} log(ζ_{p∤D}(Nm) − \frac{1}{2} ∑_{2∤N} \frac{μ(N)}{N} log (\frac{ζ^{2}}{ζ_{k}} _{p∤D}(Nm))

− \frac{1}{3} ∑_{3∤N} \frac{μ(N)}{N} log (\frac{ζ_{k}ζ}{ζ_{K}} _{p∤D}(Nm)) .

As usual, this allows us to compute sums involving ω(p), and in particular the Hardy-Littlewood constants. With the same notations as above, we obtain after some calculation

− log(H_{D}(A)) = log lim_{s→1^{+}} (\frac{ζ_{K,p∤D}(s)}{ζ_{p∤D}(s)}) + ∑_{N≥2} a(N) log(ζ_{K,p∤D}(N))

+ ∑_{N≥2} b(N) log(ζ_{k,p∤D}(N)) − ∑_{N≥2}(a(N/3) + 2b(N)) log(ζ_{p∤D}(N))

with

a(N) = \frac{1}{3N} ∑_{k|N,3∤k} μ(k)3^{N/k}

and

b(N) = \frac{1}{2N} ∑_{k|N,2∤k} μ(k)(3^{N/k} − 1) − a(N) .

If we let L_{K}(s) = ζ_{K}(s)/ζ(s), this is now the Artin L-function of the unique irreducible representation of degree 2 of Gal(K/ℚ) ≃ S_{3}, and if we set L_{k}(s) = ζ_{k}(s)/ζ(s), which is equal, up to a finite number of Euler factors to L(s, (D/n)), the above formula can be rewritten

− log(H_{D}(A)) = ∑_{N≥1} a(N) log(L_{K,p∤D}(N)) + ∑_{N≥1} b(N) log(L_{k,p∤D}(N))

+ ∑_{N≥1} c(N) log(ζ_{p∤D}(N))

with a(N) and b(N) as above and c(N) = a(N) − b(N) − a(N/3) if N > 1, while c(1) = 0.

We have already seen how to compute L_{k}(N) and ζ(N). To compute L_{K}(s), we proceed as follows. Once again we want to apply Theorem 3.1. To obtain the necessary Dirichlet series coefficients is easy, since they follow directly from the splitting behavior of the primes in the extension K/ℚ. On the other hand, the functions F_{i}(s, x) (both equal to a single function F (s, x)) may be more difficult to compute, since the gamma factor has more than a single gamma function. In fact, this is the very reason why it is simpler to compute L_{K}(s) than ζ_{K}(s) = L_{K}(s)ζ(s).

We consider first the case where D < 0, i.e. the field K is a complex cubic field. In this case, the duplication formula of the gamma function simplifies things considerably since we still have a single gamma function in the gamma factor for L_{K} (this would not be the case for ζ_{K}). More precisely, Theorem 3.1 can be applied with L_{1} = L_{2} = L_{K},

γ_{1}(s) = γ_{2}(s) = γ(s) = (2π)^{−s}|D|^{s/2}Γ(s) ,

hence

F_{i}(s, x) = (2π)^{−s}|D|^{s/2}Γ (s, \frac{2πx}{\sqrt{|D|}}) .

Writing L_{K}(s) = ∑_{n≥1} u(n)n^{−s}, where the u(n) are easily computed from the Dirichlet series expansion of ζ_{K}(s), we obtain

Γ(s)L_{K}(s) = ∑_{n≥1} \frac{u(n)}{n^{s}}Γ (s, \frac{2πn}{\sqrt{|D|}})+(2π)^{2s−1}|D|^{1/2−s} ∑_{n≥1} \frac{u(n)}{n^{1−s}} Γ (1 − s, \frac{2πn}{\sqrt{|D|}}) .

Thus, computing L_{K}(s) in this case is almost as easy as computing L_{k}(s). Note that the convergence is slightly slower, since the series converges like exp(−An/\sqrt{|D|}) instead of exp(−An^{2}/|D|), but the limitation is still the same, i.e. we need to take at least O(|D|^{1/2}) terms before having any convergence at all.

Consider now the case D > 0. Set A = π/\sqrt{D}. In this case the function L_{K} has a gamma factor A^{−s}Γ(s/2)^{2}, and this does not lead to a function F (s, x) as simple as an incomplete gamma function. The function W (t) is essentially a K-bessel function, and we find

F (s, x) = A^{−s/2}Γ(s/2)^{2} − 4 ∑_{n≥0} \frac{A^{2n}x^{2n+s}}{(2n + s)n!^{2}} (H_{n} − γ − log(Ax) + \frac{1}{2n + s})

where γ = 0.57721 … is Euler’s constant and H_{n} = ∑_{1≤k≤n} 1/k. As usual, this series converges for all x but suffers from severe cancellation problems when x is large.

In that case, we have two solutions. The first one is to use an accuracy which is slighly more than twice the needed accuracy. This enables us to compute F (s, x) for all x. The other solution is to find a generalization of the method using continued fractions that we used for the incomplete gamma function. For this, we would need to give recursions of higher order satisfied by F (s, x), and this can be done but would carry us too far for the present paper.

5.3. The case of general polynomials

We mention here briefly the case of polynomials of arbitrary degree. Let A be an irreducible content-free polynomial of degree n and Galois group G defining a number field K. Since there are exactly as many irreducible representations of G as conjugacy clases, and these conjugacy classes exactly correspond to the splitting behavior of primes in K/ℚ, as in the preceding cases we will have a system of linear equations for the quantities ∑_{ω(p)=r} p^{−m} in terms of sums of the form ∑_{k≥1} μ(k)/k log(L(k, ρ)). It can be shown that this system always has a solution, from which we deduce explicit formulas for the Hardy-Littlewood constants. We must then compute the values of the Artin L-functions L(N, ρ), and this is done using Theorem 3.1 and generalizations of the formulas found above for F (s, x). There is no problem of principle in performing all these computations, but the details become tedious, especially for the primes dividing the discriminant D of the polynomial A, since these have to be considered separately.

Laboratoire A2X, UMR 9936 du C.N.R.S., Université Bordeaux I, 351 Cours de la Libération, 33405 TALENCE Cedex, FRANCE