Problem 951: Euler Product Approximation

Mathematical Foundation

Theorem 1 (Euler Product Formula). For $\operatorname{Re}(s) > 1$,

Proof. Each factor expands as a geometric series: $(1 - p^{-s})^{-1} = \sum_{k=0}^{\infty} p^{-ks}$, convergent for $\operatorname{Re}(s) > 0$. For any finite set of primes $P$,

by the Fundamental Theorem of Arithmetic: every positive integer factors uniquely into primes, so expanding the product yields exactly one term $n^{-s}$ for each $n$ whose prime factors all lie in $P$. As $P$ exhausts all primes, the right-hand side converges to $\zeta(s)$ by absolute convergence of the Dirichlet series for $\operatorname{Re}(s) > 1$. $\square$

Lemma 1 (Monotonicity). $E(s, N) \le \zeta(s)$ for all $N$, with $E(s, N) \nearrow \zeta(s)$ as $N \to \infty$.

Proof. Each factor $(1 - p^{-s})^{-1} > 1$ for $s > 0$. The partial product $E(s, N)$ accounts for all positive integers whose prime factors are at most $N$, a subset of all positive integers. Including additional primes strictly increases the product. $\square$

Theorem 2 (Error Bound). For $s = 3$ and $N \ge 2$,

Proof. We have $\zeta(3)/E(3,N) = \prod_{p > N}(1-p^{-3})^{-1}$. Taking logarithms:

By partial summation with the Prime Number Theorem $\pi(x) \sim x/\ln x$:

Since $\ln(1 + \varepsilon) \sim \varepsilon$ for small $\varepsilon$, we obtain $\zeta(3)/E(3,N) - 1 = O(1/(N^2 \ln N))$, and therefore $\zeta(3) - E(3,N) = \zeta(3) \cdot O(1/(N^2 \ln N)) = O(1/(N^2 \ln N))$. $\square$

Lemma 2 (Numerical Precision). For $N = 10^6$, the relative error $(\zeta(3) - E(3,N))/\zeta(3)$ is approximately $3 \times 10^{-14}$. Standard IEEE 754 double-precision arithmetic (53-bit mantissa, $\sim$15.9 decimal digits) suffices for extracting $\lfloor E(3,10^6) \times 10^{15} \rfloor$.

Proof. By Theorem 2, the relative error is $\approx 3/(2 \cdot 10^{12} \cdot 6\ln 10) \approx 1.1 \times 10^{-14}$. Working in log-space (accumulating $L = \sum_{p \le N} -\ln(1-p^{-3})$), each term has relative error $\varepsilon_{\mathrm{mach}} \approx 2^{-53}$. Over $\pi(10^6) = 78498$ terms, the accumulated error in $L$ is at most $78498 \cdot 2^{-53} \approx 8.7 \times 10^{-12}$, which propagates to an absolute error of $\sim 10^{-11}$ in $E$. This is well within the tolerance for the 15th decimal digit. $\square$

Editorial

Compute floor(E(3, 10^6) * 10^15) where E(3, N) = prod_{p <= N} (1 - p^{-3})^{-1}. The Euler product for zeta(s) converges to Apery’s constant zeta(3) = 1.2020569… We accumulate the product in log-space for numerical stability: L = sum_{p <= N} -ln(1 - p^{-3}) E = exp(L) Key identity: zeta(s) = prod_p (1 - p^{-s})^{-1} (Re(s) > 1) Error bound: zeta(3) - E(3,N) ~ 4*zeta(3) / (N^2 * ln(N)). We sieve of Eratosthenes. We then log-space accumulation. Finally, recover product and extract answer.

Pseudocode

Sieve of Eratosthenes
Log-space accumulation
Recover product and extract answer

Complexity Analysis

• Time: $O(N \log \log N)$ for the Sieve of Eratosthenes, plus $O(\pi(N)) = O(N/\ln N)$ for the log-space accumulation. Total: $O(N \log \log N)$.
• Space: $O(N)$ for the sieve array.

For $N = 10^6$: the sieve dominates at $\sim 10^6$ operations; the product accumulation over $78498$ primes is negligible.

Answer

/*
 * Problem 951: Euler Product Approximation
 *
 * Compute floor(E(3, 10^6) * 10^15) where
 *   E(3, N) = prod_{p <= N, p prime} (1 - p^{-3})^{-1}
 *
 * Algorithm:
 *   1. Sieve primes up to N = 10^6
 *   2. Accumulate L = sum -ln(1 - p^{-3}) in log-space (long double)
 *   3. E = exp(L), answer = floor(E * 10^15)
 *
 * The partial Euler product converges to zeta(3) = 1.20205690315959...
 * (Apery's constant) with error O(1/(N^2 ln N)).
 *
 * Complexity: O(N log log N) for sieve, O(pi(N)) for product.
 */

#include <bits/stdc++.h>
using namespace std;

int main() {
    const int N = 1000000;

    // Sieve of Eratosthenes
    vector<bool> is_prime(N + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; (long long)i * i <= N; i++) {
        if (is_prime[i]) {
            for (int j = i * i; j <= N; j += i) {
                is_prime[j] = false;
            }
        }
    }

    // Accumulate Euler product in log-space using long double
    long double log_prod = 0.0L;
    int prime_count = 0;
    for (int p = 2; p <= N; p++) {
        if (is_prime[p]) {
            long double term = -logl(1.0L - powl((long double)p, -3.0L));
            log_prod += term;
            prime_count++;
        }
    }

    long double E = expl(log_prod);

    // Verification: zeta(3) = 1.2020569031595942854...
    long double zeta3 = 1.2020569031595942854L;
    long double error = zeta3 - E;

    cout << fixed << setprecision(18);
    cout << "E(3, " << N << ") = " << (double)E << endl;
    cout << "zeta(3)     = " << (double)zeta3 << endl;
    cout << "Error       = " << (double)error << endl;
    cout << "Primes used = " << prime_count << endl;
    cout << endl;

    // Final answer
    long long answer = (long long)(E * 1e15L);
    cout << answer << endl;

    return 0;
}

"""
Problem 951: Euler Product Approximation

Compute floor(E(3, 10^6) * 10^15) where E(3, N) = prod_{p <= N} (1 - p^{-3})^{-1}.

The Euler product for zeta(s) converges to Apery's constant zeta(3) = 1.2020569...
We accumulate the product in log-space for numerical stability:
    L = sum_{p <= N} -ln(1 - p^{-3})
    E = exp(L)

Key identity: zeta(s) = prod_p (1 - p^{-s})^{-1}  (Re(s) > 1)
Error bound: zeta(3) - E(3,N) ~ 4*zeta(3) / (N^2 * ln(N))
"""

from math import log, exp, isqrt

ZETA3 = 1.2020569031595942854169

def sieve_of_eratosthenes(n: int) -> list[int]:
    """Return all primes up to n using the Sieve of Eratosthenes."""
    is_prime = bytearray(b'\x01') * (n + 1)
    is_prime[0] = is_prime[1] = 0
    for i in range(2, isqrt(n) + 1):
        if is_prime[i]:
            is_prime[i * i::i] = bytearray(len(is_prime[i * i::i]))
    return [i for i in range(2, n + 1) if is_prime[i]]

def euler_product(s: float, N: int) -> tuple[float, list[tuple[int, float]]]:
    """Compute E(s, N) = prod_{p<=N} (1 - p^{-s})^{-1} in log-space.

    Returns:
        (E, checkpoints): the product value and a list of (prime, running_product)
    """
    primes = sieve_of_eratosthenes(N)
    log_prod = 0.0
    checkpoints = []
    for p in primes:
        log_prod += -log(1 - p ** (-s))
        if p <= 500 or p % 5000 < 10:
            checkpoints.append((p, exp(log_prod)))
    return exp(log_prod), checkpoints

def solve(N: int = 10**6) -> int:
    """Compute floor(E(3, 10^6) * 10^15)."""
    E, _ = euler_product(3.0, N)
    return int(E * 10**15)

# --- Compute answer ---
N = 10**6
E_full, checkpoints = euler_product(3.0, N)
answer = int(E_full * 10**15)

assert abs(E_full - ZETA3) < 1e-12, f"E(3,10^6) = {E_full}, expected ~{ZETA3}"
print(answer)