Problem 978: Bernoulli Number Modular

Mathematical Foundation

Definition 1 (Bernoulli Numbers). The Bernoulli numbers $\{B_n\}_{n \ge 0}$ are the rational numbers defined by $B_0 = 1$ and the recurrence:

Equivalently, they are defined by the exponential generating function $\frac{t}{e^t - 1} = \sum_{n=0}^{\infty} B_n \frac{t^n}{n!}$.

Theorem 1 (Vanishing of Odd Bernoulli Numbers). For all odd $n \ge 3$, $B_n = 0$.

Proof. Consider $f(t) = \frac{t}{e^t - 1} + \frac{t}{2} = \frac{t(e^t + 1)}{2(e^t - 1)} = \frac{t}{2}\coth\frac{t}{2}$. Since $\coth$ is an even function, $f(t)$ is an even function of $t$. Therefore $f(t) = \sum_{n=0}^{\infty} B_n \frac{t^n}{n!} + \frac{t}{2}$ is even, which means all odd-indexed coefficients of $\frac{t}{e^t-1}$ vanish except $B_1 = -1/2$ (the $t/2$ term). Hence $B_n = 0$ for odd $n \ge 3$. $\square$

Theorem 2 (Von Staudt—Clausen). For even $n \ge 2$:

where $A_n \in \mathbb{Z}$. Consequently, the denominator of $B_n$ (in lowest terms) is:

Proof. We use the fact that for a prime $p$ and even $n \ge 2$:

which follows from the Clausen—von Staudt congruence, provable via the power sum formula $S_n(p) = \sum_{j=0}^{p-1} j^n \equiv -1 \pmod{p}$ when $(p-1) \mid n$ (by Fermat’s little theorem) and $S_n(p) \equiv 0 \pmod{p}$ otherwise. Since $B_n = \frac{1}{n+1}\sum_{j=0}^{n} \binom{n}{j} S_j(m)$ for suitable $m$, one extracts the $p$-adic valuation of $B_n$ to obtain the stated formula. The full proof is in Ireland & Rosen, A Classical Introduction to Modern Number Theory, Chapter 15. $\square$

Lemma 1 (Sign Pattern). For even $n \ge 2$: $(-1)^{n/2+1} B_n > 0$. That is, $B_2 > 0$, $B_4 < 0$, $B_6 > 0$, $B_8 < 0$, $\ldots$

Proof. From the functional equation of $\zeta(s)$ and the identity $B_n = (-1)^{n/2+1} \frac{2 \cdot n!}{(2\pi)^n} \zeta(n)$ for even $n \ge 2$. Since $\zeta(n) > 0$ for $n \ge 2$, the sign is $(-1)^{n/2+1}$. $\square$

Theorem 3 (Growth of Bernoulli Numerators). For even $n$,

as $n \to \infty$. In particular, $|p_n|$ grows super-exponentially.

Proof. From $B_n = (-1)^{n/2+1} \frac{2 \cdot n!}{(2\pi)^n} \zeta(n)$ and $\zeta(n) \to 1$ as $n \to \infty$. $\square$

Editorial

Compute sum_{n=0}^{200} (|p_n| + q_n) mod (10^9 + 7), where B_n = p_n / q_n is the n-th Bernoulli number in lowest terms. Bernoulli numbers are defined by the exponential generating function: t / (e^t - 1) = sum B_n * t^n / n! Key properties: denom(B_{2n}) = product of primes p where (p-1) | 2n. We compute B_0, B_1, …, B_N using exact rational arithmetic. We then reduce B[n] to lowest terms. Finally, sum |p_n| + q_n modulo MOD.

Pseudocode

Compute B_0, B_1, ..., B_N using exact rational arithmetic
Reduce B[n] to lowest terms
Sum |p_n| + q_n modulo MOD

Complexity Analysis

• Time: $O(N^2 \cdot M(N))$ where $M(N)$ is the cost of arithmetic on $O(N \log N)$-bit integers. The recurrence has $O(N)$ nonzero terms (only even indices plus $B_1$), each requiring a sum of $O(N)$ terms. Each rational arithmetic operation involves numbers with $O(N \log N)$ bits (by Theorem 3), so $M(N) = O(N \log N \cdot \log(N \log N))$ using fast multiplication. Total: roughly $O(N^3 \log^2 N)$.
• Space: $O(N^2 \log N)$ bits to store all $N$ Bernoulli numbers.

For $N = 200$, this is entirely tractable with Python’s fractions.Fraction.

Answer

#include <iostream>

// Reference executable for problem_978.
// The mathematical derivation is documented in solution.md and solution.tex.
static const char* ANSWER = "952347340";

int main() {
    std::cout << ANSWER << '\n';
    return 0;
}

"""
Problem 978: Bernoulli Number Modular Sum

Compute sum_{n=0}^{200} (|p_n| + q_n) mod (10^9 + 7),
where B_n = p_n / q_n is the n-th Bernoulli number in lowest terms.

Bernoulli numbers are defined by the exponential generating function:
    t / (e^t - 1) = sum B_n * t^n / n!

Key properties:
    - B_0 = 1, B_1 = -1/2, B_n = 0 for odd n >= 3
    - |B_{2n}| grows roughly as 2 * (2n)! / (2*pi)^{2n}
    - Denominators follow the Von Staudt-Clausen theorem:
      denom(B_{2n}) = product of primes p where (p-1) | 2n

Answer: computed via exact Fraction arithmetic

Methods:
    - compute_bernoulli(n_max): Compute B_0..B_{n_max} via recurrence
    - von_staudt_clausen_denominator(n): Predict denominator via theorem
    - verify_small_values(): Check known Bernoulli numbers
"""

from fractions import Fraction
from math import comb, log10

MOD = 10**9 + 7


def compute_bernoulli(n_max):
    """Compute B_0 to B_{n_max} using the standard recurrence."""
    B = [Fraction(0)] * (n_max + 1)
    B[0] = Fraction(1)
    for n in range(1, n_max + 1):
        B[n] = -sum(Fraction(comb(n + 1, k)) * B[k] for k in range(n)) / (n + 1)
    return B


def von_staudt_clausen_denominator(n):
    """
    For even n >= 2, denom(B_n) = product of primes p where (p-1) | n.
    """
    if n == 0:
        return 1
    if n == 1:
        return 2
    if n % 2 == 1:
        return 1  # B_n = 0 for odd n >= 3
    # Find all primes p where (p-1) divides n
    result = 1
    # Check primes up to n+1
    for p in range(2, n + 2):
        # Simple primality check
        if p < 2:
            continue
        is_prime = True
        for d in range(2, int(p**0.5) + 1):
            if p % d == 0:
                is_prime = False
                break
        if is_prime and n % (p - 1) == 0:
            result *= p
    return result


def verify_small_values(B):
    """Check known Bernoulli numbers."""
    assert B[0] == Fraction(1)
    assert B[1] == Fraction(-1, 2)
    assert B[2] == Fraction(1, 6)
    assert B[4] == Fraction(-1, 30)
    assert B[6] == Fraction(1, 42)
    assert B[8] == Fraction(-1, 30)
    assert B[10] == Fraction(5, 66)
    # Odd Bernoulli numbers (>= 3) are zero
    for n in range(3, 20, 2):
        assert B[n] == 0
    return True

# Computation

B = compute_bernoulli(200)
verify_small_values(B)

# Verify Von Staudt-Clausen for even indices
for n in range(0, 50, 2):
    if n >= 2:
        assert B[n].denominator == von_staudt_clausen_denominator(n), \
            f"VSC failed at n={n}: {B[n].denominator} vs {von_staudt_clausen_denominator(n)}"

answer = 0
for n in range(201):
    p, q = B[n].numerator, B[n].denominator
    answer = (answer + (abs(p) + q) % MOD) % MOD

print(answer)