Project Euler

Dirichlet Series Evaluation

Define D(s) = sum_(n=1)^N (mu(n))/(n^s) where mu is the Mobius function. Find floor(D(2) * 10^12) for N = 10^7.

Source sync Apr 19, 2026

Problem #0943

Level Level 39

Solved By 102

Languages C++, Python

Answer 1038733707

Length 175 words

number_theorylinear_algebrageometry

Given two unequal positive integers $a$ and $b$, we define a self-describing sequence consisting of alternating runs of $a$s and $b$s. The first element is $a$ and the sequence of run lengths is the original sequence.

For $a=2, b=3$, the sequence is: \[2, 2, 3, 3, 2, 2, 2, 3, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 3, 2, 2, 2, 3, 3, 3,...\] The sequence begins with two $2$s and two $3$s, then three $2$s and three $3$s, so the run lengths $2, 2, 3, 3, ...$ are given by the original sequence.

Let $T(a, b, N)$ be the sum of the first $N$ elements of the sequence. You are given $T(2,3,10) = 25$, $T(4,2,10^4) = 30004$, $T(5,8,10^6) = 6499871$.

Find $\sum T(a, b, 22332223332233)$ for $2 \le a \le 223$, $2 \le b \le 223$ and $a \neq b$. Give your answer modulo $2233222333$.

Problem 943: Dirichlet Series Evaluation

Mathematical Analysis

The Mobius Function

Definition. $\mu(n) = (-1)^k$ if $n$ is a product of $k$ distinct primes, $\mu(n) = 0$ if $n$ has a squared prime factor, and $\mu(1) = 1$ .

Connection to the Riemann Zeta Function

Theorem. $\sum_{n=1}^{\infty} \frac{\mu(n)}{n^s} = \frac{1}{\zeta(s)}$ for $\text{Re}(s) > 1$ .

For $s = 2$ : $\frac{1}{\zeta(2)} = \frac{6}{\pi^2} \approx 0.607927...$

Truncated Sum

$D(2) = \sum_{n=1}^{N} \frac{\mu(n)}{n^2} \approx \frac{6}{\pi^2} - \sum_{n>N} \frac{\mu(n)}{n^2}$ .

The tail is $O(1/N)$ , so for $N = 10^7$ , $D(2) \approx 6/\pi^2$ to about 7 decimal places.

Editorial

Compute D(s) = sum_{n=1}^{N} mu(n) / n^s for s = 2, which converges to 6/pi^2 (the reciprocal of the Riemann zeta function at s=2). The Mobius function mu(n) is: mu(1) = 1 mu(n) = (-1)^k if n is a product of k distinct primes mu(n) = 0 if n has a squared prime factor The target is floor(D(2) * 10^12) for large N. Results:. We sieve $\mu(n)$ for $n \leq N$ using a linear sieve. Finally, compute $\sum \mu(n)/n^2$ using floating-point or exact rational arithmetic.

Pseudocode

Sieve $\mu(n)$ for $n \leq N$ using a linear sieve
Compute $\sum \mu(n)/n^2$ using floating-point or exact rational arithmetic
Floor the result times $10^{12}$

The Mobius Sieve

Initialize $\mu(1) = 1$ . For each prime $p$ : for multiples $m$ of $p$ , multiply $\mu(m)$ by $-1$ . For multiples of $p^2$ , set $\mu(m) = 0$ .

Proof of Correctness

Mobius sieve: Standard.
Sum: Floating-point with sufficient precision, or use integer arithmetic.

Complexity Analysis

Sieve: $O(N \log \log N)$ .
Sum: $O(N)$ .

Answer

\boxed{1038733707}

C++ project_euler/problem_943/solution.cpp

#include <bits/stdc++.h>
using namespace std;
int main(){
    const int N=1000000;
    vector<int> mu(N+1,0);
    mu[1]=1;
    vector<bool> is_p(N+1,true);
    vector<int> primes;
    for(int i=2;i<=N;i++){
        if(is_p[i]){primes.push_back(i);mu[i]=-1;}
        for(int p:primes){
            if((long long)i*p>N) break;
            is_p[i*p]=false;
            if(i%p==0){mu[i*p]=0;break;}
            else mu[i*p]=-mu[i];
        }
    }
    double D=0;
    for(int n=1;n<=N;n++) D+=mu[n]/(double)n/(double)n;
    cout<<fixed<<setprecision(12)<<D<<endl;
    return 0;
}

Python project_euler/problem_943/solution.py

"""
Problem 943: Dirichlet Series Evaluation

Compute D(s) = sum_{n=1}^{N} mu(n) / n^s  for s = 2, which converges to 6/pi^2
(the reciprocal of the Riemann zeta function at s=2).

The Mobius function mu(n) is:
    mu(1) = 1
    mu(n) = (-1)^k if n is a product of k distinct primes
    mu(n) = 0      if n has a squared prime factor

The target is floor(D(2) * 10^12) for large N.

Results:
    - D(2) -> 6/pi^2 = 0.607927101854...
    - Convergence is O(1/N) for partial sums

Methods:
    1. compute_mobius       -- linear sieve for mu(n)
    2. partial_dirichlet    -- partial sums of mu(n)/n^s
    3. mobius_distribution  -- counts of mu(n) = -1, 0, +1
    4. convergence_rate     -- |D_N(2) - 6/pi^2| analysis
"""
import numpy as np

def compute_mobius(N):
    """Compute mu(n) for n=0..N using a linear sieve."""
    mu = [0] * (N + 1)
    mu[1] = 1
    is_prime = [True] * (N + 1)
    primes = []
    for i in range(2, N + 1):
        if is_prime[i]:
            primes.append(i)
            mu[i] = -1
        for p in primes:
            if i * p > N:
                break
            is_prime[i * p] = False
            if i % p == 0:
                mu[i * p] = 0
                break
            else:
                mu[i * p] = -mu[i]
    return mu

def partial_dirichlet(mu, s, N):
    """Compute sum_{n=1}^{N} mu(n) / n^s."""
    return sum(mu[n] / n**s for n in range(1, N + 1))

def mobius_distribution(mu, N):
    """Count how many n in [1,N] have mu(n) = -1, 0, +1."""
    neg = sum(1 for n in range(1, N + 1) if mu[n] == -1)
    zero = sum(1 for n in range(1, N + 1) if mu[n] == 0)
    pos = sum(1 for n in range(1, N + 1) if mu[n] == 1)
    return neg, zero, pos

def convergence_errors(mu, checkpoints):
    """Compute |D_N(2) - 6/pi^2| at each checkpoint N."""
    target = 6.0 / np.pi**2
    errors = []
    s = 0.0
    idx = 0
    for n in range(1, checkpoints[-1] + 1):
        s += mu[n] / n**2
        if idx < len(checkpoints) and n == checkpoints[idx]:
            errors.append(abs(s - target))
            idx += 1
    return errors

# Verification
mu_small = compute_mobius(20)
assert mu_small[1] == 1
assert mu_small[2] == -1   # prime
assert mu_small[3] == -1   # prime
assert mu_small[4] == 0    # 4 = 2^2
assert mu_small[5] == -1   # prime
assert mu_small[6] == 1    # 2*3, two distinct primes
assert mu_small[8] == 0    # 8 = 2^3
assert mu_small[10] == 1   # 2*5

# Computation
N = 100000
mu = compute_mobius(N)
D2 = partial_dirichlet(mu, 2, N)
print(D2)