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

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:

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)$.

• Fortunately, evaluating $L(1, \chi_D)$ is actually very straightforward. If you want to do this in a very sophisticated way, you could use the class number formula. But you could also evaluate it explicitly and simply through the use of Hurwitz zeta functions. $L(s, \chi)$ is a finite sum of Hurwitz zeta functions $\zeta(s, a)$, and $\zeta(0, a)$ is given explicitly by a Bernouilli number (like the regular zeta function). This will get you $L(0, \chi)$, and the functional equation will then get you $L(1, \chi)$. – davidlowryduda Mar 30 '18 at 0:35
• If you think this line of reasoning is interesting, I could return and flesh out an answer distilling the relevant parts. But you might think it is not elementary enough, perhaps? – davidlowryduda Mar 30 '18 at 0:36
• @mixedmath Thank you, that's great to know! I'm willing to accept on faith, any theorem, identity, or formula, no matter how “deep” (at least, if I can understand the terms and what is being claimed!), but ultimately I want to be able to turn what I learn into a computer program, so I'd hope for sufficient detail… assuming that's not too much work, I'd definitely love to hear it! – ShreevatsaR Mar 30 '18 at 0:44
• If you need to list all the primes up to a number using a sieve algorithm is faster than primality checking one by one. – Sophie Apr 4 '18 at 23:10

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$, $\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).

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.

• Thank you! Are there any additional complexities when $3$ in the example is replaced by a higher number: that is, when working, instead of with $\left(\frac{3}{p}\right)$, with say $\left(\frac{437}{p}\right)$ (here $437 = 19 \times 23$, the discriminant of a quadratic polynomial that I had happened to start this exploration with and am still interested in)? (I think in my example I was working with $-3$ instead of $3$ btw.) – ShreevatsaR Mar 30 '18 at 22:28
• No, there will be no significant differences for higher numbers (unless you happen to have a discriminant which is divisible by a nontrivial square number. Then the Dirichlet character isn't primitive, and a few minor changes are necessary). – davidlowryduda Mar 31 '18 at 0:57
• Could you clarify what exactly $\chi$ is (for the case of $\left(\frac{D}{p}\right)$), and what's the $q$, and how to find them? Is $\chi(n)=\left(\frac{D}{n}\right)$ even if $n$ may not be prime? For example, for $D=-3$, indeed I can see that we care about the character of modulus $3$, given in the answer. But for $D=3$, we have $\left(\frac{D}{p}\right)=1$ when $p \equiv1\pmod{12}$ and $\left(\frac{D}{p}\right)=-1$ when $p \equiv7\pmod{12}$, so I think in fact $q=12$. This does not fit either $q\equiv1\pmod4$ or $q\equiv3\pmod4$. (I think OEIS A037449 is related.) – ShreevatsaR Apr 8 '18 at 6:13