Not-too-slow computation of Euler products / singular series

Question

I'd like to compute, to at least a few digits of accuracy, the constants that arise in Hardy-Littlewood conjecture F / Bateman-Horn conjecture, in particular for just a single quadratic polynomial.

Specifically, consider (for instance) the following number, which arises when studying primes in the polynomial $n^2 + n + 1$:

$$ \begin{align} C_{1} &= \frac{2}{1} \frac{2}{2} \frac{5}{4} \frac{5}{6} \frac{11}{10} \frac{11}{12} \frac{17}{16} \frac{17}{18} \frac{23}{22} \frac{29}{28} \frac{29}{30} \frac{35}{36} \frac{41}{40} \cdots \\ &= \prod_{p \equiv 1 \pmod3}\left(\frac{p - 2}{p - 1}\right) \prod_{p \equiv 2 \pmod3}\left(\frac{p}{p - 1}\right) \\ &= \prod_{p \equiv 1 \pmod3}\left(\frac{1 - 2/p}{1 - 1/p}\right) \prod_{p \equiv 2 \pmod3}\left(\frac{1}{1 - 1/p}\right) \\ \end{align}$$ where the products run over all primes that leave remainder $1$ and $2$ when divided by $3$, respectively. (Note: we cannot really write it as two separate infinite products; the order matters, so excuse the notation.)

Now, if I try to compute this the naive way:

#include <cstdio>
#include "FJ64_16k.h" // For a fast `is_prime` function: http://ceur-ws.org/Vol-1326/020-Forisek.pdf

int main() {
  double ans = 0.5;
  for (uint64_t p = 0; ; ++p) {
    if (p % 1000000 == 0) {
      printf("%lld %.9f\n", p, ans);
    }
    if (!is_prime(p)) continue;
    if (p % 3 == 1) ans *= (p - 2.0) / (p - 1.0);
    if (p % 3 == 2) ans *= p / (p - 1.0);
  }
}

Then even after computing with primes up to 100 million, we seem to have only have about four decimal digits of accuracy:

99000000 1.120724721
100000000 1.120725012
101000000 1.120727310

For another example, if we start with the polynomial $n^2 + 5n + 1$ (discriminant $21$), then we need to compute (something like) the constant $$C_2 = \prod_{p \in P_1 \cap Q_1 \cup P_2 \cap Q_2} \frac{(p-2)}{(p-1)} \prod_{p \in P_1 \cap Q_2 \cup P_2 \cap Q_1} \frac{p}{(p-1)} $$ where $P_1 = \{p \equiv 1 \pmod 3\}$, $P_2 = \{p \equiv 2 \pmod 3\}$, $Q_1 = \{p \equiv \text{$1$, $2$, or $4$}\pmod 7\}$, $Q_2 = \{p \equiv \text{$3$, $5$ or $6$}\pmod 7\}$ which is messier: we do one thing for $p \equiv 1, 4, 5, 16, 17, 20\pmod 21$ and another for $p \equiv 2, 8, 10, 11, 13, 19$.

---

Question: How can we compute these constants reasonably quickly?

I have seen the related questions on MathOverflow, and tried to read their answers and some of the mentioned references:

• Calculating the infinite product from the Hardy-Littlewood Conjecture F (2010)

• Calculating the constant in the Bateman-Horn-Stemmler conjecture (2011)

• High precision computation of Hardy-Littlewood constants (DVI) by Henri Cohen

• Some Constants from Number Theory (2001?) (PS) by Pascal Sebah and Xavier Gourdon

• Approximation of singular series and automata (2001), by Pieter Moree

• New quadratic polynomials with high densities of prime values (2003), by Jacobson and Williams

• Quadratic Polynomials which have a High Density of Prime Values (1990) by Fung and Williams

• On the Conjecture of Hardy & Littlewood concerning the Number of Primes of the $n^2 + a$ (1960), by Shanks

• Supplementary data and remarks concerning a Hardy-Littlewood conjecture (1963), by Shanks

• The calculation of certain Dirichlet series (1963), by Shanks and Wrench

• Note on the number of primes of the form $n^2 + 1$ (1922) by A. E. Western

• Some problems of ‘Partitio numerorum’; III: On the expression of a number as a sum of primes (1923), Hardy and Littlewood (the source of Conjecture F)

But the papers take quite a lot of effort for me to follow (I've never encountered Dirichlet L-series before, for example), and I keep suspecting there should be something simpler, if we don't want so much. (The older papers are simpler, but Western's paper simply says “by the use of a transformation suggested by Mr. Littlewood” without a reference or explaining anything.) In particular, I don't care about getting thousands of digits, I'll be happy if I can get, say six (ten would be fantastic). (Or even 4 digits, if they can be computed a lot faster.) And I don't care about computing such series that arise out of really large numbers; like the $3$ and $7$ in the second example above I might have at most two-digit factors.

With these looser constraints, is there something simpler that suffices, to compute these constants? What I'm looking for is either:

• ideally, a method of computation with sufficient detail, and which can be translated into a few lines of code that compute it from scratch (as in the program above) (after some reasonable amount of work on paper if necessary),
• or, the same, but it's ok to rely on certain constants (like $\zeta(2)$ say) that I can just look up and hard-code into the program, or on some existing functions (such as the logarithmic integral say), if they are well-known functions commonly available (in software libraries that are easy to install or have an online interface).

I'm looking for an elementary exposition at undergraduate level, say. For instance, the answer by KConrad makes sense (though it's not clear how good the convergence is), but I don't know how to compute $L(1,\chi_D)$.

Answer 1

The intent here is to build off of Keith's answer, and thus discuss how to evaluate the $L$-function at $1$ and to discuss the rate of convergence.

Evaluating $L(1, \chi)$

We will use the following theorem.

Theorem

Let $\chi$ be a primitive Dirichlet character of modulus $q$, and let $j$ be any positive integer. Then $$ L(1-j,\chi) = - \frac{q^{j-1}}{j} \sum_{1 \leq a \leq q} \chi(a)B_{j}(\tfrac{a}{q}),$$ where $B_j(x)$ is the $j$th Bernoulli polynomial, defined by $$\frac{t e^{Xt}}{e^t-1} = \sum_{n=0}^\infty B_n(X) \frac{t^n}{n!}.$$

As a corollary, choosing $j = 1$ and evaluating that $B_1(x) = x - \tfrac{1}{2}$, we have that $$ L(0, \chi) = - \sum_{1 \leq a \leq q} \chi(a) \left(\frac{a}{q} - \frac{1}{2}\right).$$ If $\chi \neq 1$, then $\sum_a \chi(a) = 0$, so that this simplifies to $$ L(0, \chi) = - \sum_{1 \leq a \leq q} \chi(a) \frac{a}{q}.$$ This is tremendously straightforward to actually compute.

Recall the functional equation for a Dirichlet $L$-function, which states that $$ \Lambda(s, \chi) := \left( \frac{\pi}{q} \right)^{-(s + A)/2} \Gamma(\tfrac{s + A}{2})L(s, \chi),$$ where $$ A = \begin{cases} 0 & \text{if } \chi(-1) = 1 \\ 1 & \text{if } \chi(-1) = -1 \end{cases},$$ satisfies the functional equation $$ \Lambda(1-s, \overline{\chi}) = \frac{i^A \sqrt q}{\tau(\chi)} \Lambda(s, \overline{\chi}).$$

There are several pieces here to disentangle. Note that $\overline{\chi}$ is the character obtained by (complex) conjugating $\chi$. For the example in your post, $\chi = \overline{\chi}$ since $\chi$ is real-valued. The function $\tau(\chi)$ is a Gauss sum, which are a bit annoying in general but ultimately they are a short finite sum. Note that when $\chi$ is quadratic, as in your example problem in the post, the Gauss sum is already known and $\tau(\chi) = \sqrt q$ if $q \equiv 1 \bmod 4$ and $\tau(\chi) = i \sqrt q$ if $q \equiv 3 \bmod 4$.

An example application: your question

Let's specialize to your post now. The underlying character $\chi$ is the primitive character mod $3$, given by $$ \chi(a) = \begin{cases} 0 & \text{if } a \equiv 0 \bmod 3 \\ 1 & \text{if } a \equiv 1 \bmod 3 \\ -1 & \text{if } a \equiv 2 \bmod 3 \end{cases}.$$ For concreteness, from the Euler product representation, we see that $$ L(1, \chi) = \prod_p (1 - \chi(p)/p)^{-1} = \prod_{p \equiv 1 \bmod 3} \left(1 - \frac{1}{p}\right) \prod_{p \equiv 2 \bmod 3} \left(1 + \frac{1}{p}\right).$$

From the corollary to the theorem above, we see that $$L(0, \chi) = -\frac{1}{3} \left(1 - 2\right) = \frac{1}{3}.$$ Since $3 \equiv 3 \bmod 4$, we have that $\tau(\chi) = i \sqrt 3$, and $A = 1$ (as appears in the functional equation above). Note again that $\overline{\chi} = \chi$.

Thus $$\Lambda(1, \chi) = \frac{i \sqrt 3}{i \sqrt 3} \Lambda(0, \chi),$$ or rather $$ L(1, \chi) \left( \frac{\pi}{3} \right)^{-1} \Gamma(1) = L(0, \chi) \left( \frac{\pi}{3} \right)^{-1/2} \Gamma(1/2).$$ Rearranging, we find that $$ L(1, \chi) = L(0, \chi) \left( \frac{\pi}{3} \right)^{1/2} \Gamma(1/2) / \Gamma(1).$$ I note that $\Gamma(1/2) = \pi^{1/2}$, $\Gamma(1) = 0! = 1$, and (from above) $L(0, \chi) = 1/3$, so ultimately $$ L(1, \chi) = \frac{1}{3} \frac{\pi}{\sqrt 3} = \frac{\pi\sqrt 3}{9}.$$

(Aside: this can be quickly verified in sage, and this is the output of quadratic_L_function_exact(1, -3). The - in -3 comes from a choice of numerator and denominator in specifying Dirichlet characters).

Convergence of Keith's answer

Keith Conrad showed that one reduces this computation to computing $$ \frac{1}{L(1,\chi_D)}\prod_{p} \frac{1 - (D|p)/(p-1)}{1 - (D|p)/p} = \frac{1}{L(1,\chi_D)}\prod_{p} \left(1 - \left(\frac{\chi(p)}{(p-1)(p-\chi(p))} \right) \right)$$ (where I've kept the simplified RHS explicit, instead of replacing it with a big Oh term). Generically, the convergence of an infinite product, written in this form, is the same speed as the convergence of the corresponding infinite series $$ \sum_{p} \left( \frac{\chi(p)}{(p-1)(p-\chi(p))} \right),$$ which converges at a rate at worst the same as $$ \sum_n \frac{1}{n^2}.$$ But in fact the series is approximately alternating, and the rate of convergence will be much faster.

From heuristics coming from the alternating series test, you might expect that to compute the first $d$ digits, you should need to use the first $\sqrt d$ terms or so.

Proof of Theorem

I will later edit in some exposition of this theorem, but this follows from representing the Dirichlet $L$-function as a sum of Hurwitz zeta functions and then evaluating each of these Hurwitz zeta functions at negative integer arguments. This is roughly done in Apostol's Intro to Analytic Number Theory, chapter 12.

But the point is that almost each step is actually very elementary, with the possible exception of proving the functional equations. Apostol actually uses the functional equations of Hurwitz zeta functions to prove the functional equation for Dirichlet $L$-functions, but some complex analysis gets involved.