What is Brun’s Constant?
Brun’s constant $B$ is defined as the sum of the reciprocals of all twin primes:
$$B = \left(\frac{1}{3} + \frac{1}{5}\right) + \left(\frac{1}{5} + \frac{1}{7}\right) + \left(\frac{1}{11} + \frac{1}{13}\right) + \cdots \approx 1.9021605$$
Unlike the harmonic series over all primes (which diverges), this series converges — a deep and surprising result proved by Viggo Brun in 1919. But it converges extremely slowly, and optimizing how we estimate it numerically is a beautiful computational challenge.
The Problem: How Slowly Does It Converge?
The partial sum up to $N$ is:
$$B(N) = \sum_{\substack{p \leq N \ p,, p+2 \text{ both prime}}} \left(\frac{1}{p} + \frac{1}{p+2}\right)$$
The error is known to behave roughly as:
$$B - B(N) \sim \frac{C}{\log^2 N}$$
for some constant $C$. This means to gain one extra decimal digit of accuracy, you need to increase $N$ exponentially. Our goal is to compute $B(N)$ efficiently and study this convergence.
Concrete Example Problem
Problem: Compute $B(N)$ for $N = 10^3, 10^4, 10^5, 10^6, 10^7$. Estimate the convergence rate empirically, fit the error model $\varepsilon(N) \sim C / \log^\alpha N$, and visualize everything including a 3D surface of the partial sum landscape.
Source Code
1 | # ============================================================ |
Code Walkthrough
1. Segmented Sieve of Eratosthenes
1 | def segmented_sieve(limit): |
This is the computational backbone. Instead of calling isprime() per number (which is $O(N \sqrt{N})$ in the naive case), we generate all primes up to $N+2$ in one pass using numpy slice assignment: is_prime[i*i::i] = False. Total cost: $O(N \log \log N)$.
2. Vectorized Brun Sum
1 | twin_mask = is_prime[primes + 2] & (primes + 2 <= N + 2) |
Rather than looping, we mask the prime array to find $p$ such that $p+2$ is also prime, then compute the entire sum with np.sum — single vectorized operation, no Python-level loop.
3. Error Model Fitting
The theoretical error behaves as:
$$\varepsilon(N) = B - B(N) \sim \frac{C}{\log^\alpha N}$$
Taking logarithms:
$$\log \varepsilon = \log C - \alpha \cdot \log(\log N)$$
This is a linear regression in $\log$–$\log$ space, handled by np.polyfit. The fitted $\alpha$ tells us empirically how fast the tail decays.
4. Cumulative Tracking
We record $B$ not just at milestone $N$ values, but at every twin prime pair, giving a smooth convergence curve. This reveals the “staircase” structure — $B(N)$ is piecewise constant between twin prime pairs and jumps at each new pair.
Graph Explanations
Plot 1 — $B(N)$ vs $N$ (semi-log): The partial sum creeps toward the dashed yellow line ($B \approx 1.9022$) but never reaches it in finite computation. The approach is visibly slow even on a log scale.
Plot 2 — Error Decay (log–log): The magenta dots are observed errors; the orange dashed line is the fitted model $\varepsilon \sim C / \log(N)^\alpha$. A straight line in log–log confirms the power-law decay in $\log N$.
Plot 3 — Cumulative Build-up: Each jump corresponds to one twin prime pair contributing $1/p + 1/(p+2)$. Early jumps (small $p$) are large; later jumps become vanishingly small — a beautiful visualization of the series thinning out.
Plot 4 — Gaps Between Twin Prime Pairs: The gaps grow roughly linearly (in a statistical sense), consistent with the Hardy–Littlewood conjecture that twin primes thin out as $\sim 2C_2 N / \log^2 N$.
Plot 5 — Successive Error Ratios: If $\varepsilon(N) \sim C/\log^\alpha N$, then going from $N = 10^k$ to $N = 10^{k+1}$ should reduce error by a factor of $((k+1)/k)^\alpha$. The gold line shows these empirical ratios — modest improvement each decade, confirming the brutal slowness.
Plot 6 — 3D Error Surface: The surface $\varepsilon = C / \log^\alpha N$ is plotted over $(\log_{10} N,, \alpha)$ space. The cyan dots mark our actual observations projected onto the fitted $\alpha$ plane. The steep slope in the $\alpha$ direction shows how sensitively the convergence speed depends on the true exponent.
Key Takeaway
$$\boxed{B - B(N) \approx \frac{C}{\log^\alpha N}, \quad \alpha \approx 2}$$
To gain one more decimal digit of $B$, you need $N$ to grow by a factor of roughly $e^{\sqrt{10}} \approx 23,000$. This is the fundamental obstacle — Brun’s constant is numerically accessible only through careful asymptotic extrapolation, not brute-force summation.
N B(N) Error Twin pairs Time(s)
-----------------------------------------------------------------
1,000 1.5180324636 0.3841281195 35 0.001
10,000 1.6168935574 0.2852670257 205 0.000
100,000 1.6727995848 0.2293609983 1,224 0.001
1,000,000 1.7107769308 0.1913836523 8,169 0.005
10,000,000 1.7383570439 0.1638035392 58,980 0.065
Fitted error model: ε(N) ≈ 2.6572 / log(N)^1.0025
Theoretical prediction: α ≈ 2 (found α ≈ 1.0025)
[Figure: brun_constant_analysis.png]
