Problem 487: Sums of Power Sums

Mathematical Foundation

Theorem 1 (Faulhaber’s formula). For any non-negative integer $k$ and positive integer $n$,

where $B_j$ are the Bernoulli numbers defined by $B_0 = 1$ and the recurrence

Proof. Consider the exponential generating function $\frac{t}{e^t - 1} = \sum_{n=0}^{\infty} B_n \frac{t^n}{n!}$. Multiplying both sides by $e^{nt} - 1)/(e^t - 1)$:

Expanding both sides in powers of $t$ and comparing the coefficient of $t^k/k!$ yields

Since $f_k(n) = \sum_{i=1}^n i^k = \sum_{i=0}^n i^k$, and the $i=0$ term vanishes for $k \ge 1$, the formula holds (with a minor adjustment for $k=0$). The Bernoulli recurrence follows from setting $n=1$ in the identity $\sum_{j=0}^k \binom{k+1}{j}B_j = (k+1) \cdot 0^k$ for $k \ge 1$ (since $f_k(0) = 0$), giving $\sum_{j=0}^k \binom{k+1}{j}B_j = 0$. $\square$

Theorem 2 (Power sums modulo a prime). For a prime $p$ and integer $m > 0$,

Proof. The multiplicative group $\mathbb{F}_p^*$ is cyclic of order $p - 1$. Let $g$ be a generator. If $(p-1) \mid m$, then $i^m = 1$ for all $i \in \mathbb{F}_p^*$, so the sum is $p - 1 \equiv -1 \pmod{p}$.

If $(p-1) \nmid m$, then $g^m \ne 1$. The sum $T = \sum_{i=1}^{p-1} i^m = \sum_{j=0}^{p-2} g^{jm}$ is a geometric series:

since $g^{(p-1)m} = (g^{p-1})^m = 1^m = 1$. $\square$

Theorem 3 (Evaluation of $S_k(p)$). Substituting Faulhaber’s formula and exchanging sums:

By Theorem 2, the inner sum $\sum_{i=1}^{p-1} i^{k+1-j}$ is $-1$ if $(p-1) \mid (k+1-j)$ and $k+1-j > 0$, and $0$ otherwise. (When $j = k+1$, the exponent is 0, and $\sum_{i=1}^{p-1} 1 = p - 1 \equiv -1$.) Therefore,

Proof. Direct substitution and application of Theorem 2. The factor $1/(k+1)$ is computed as the modular inverse $(k+1)^{-1} \bmod p$ (valid when $p \nmid k+1$). $\square$

Lemma 1 (Bernoulli number recurrence modulo $p$). The Bernoulli numbers $B_0, B_1, \ldots, B_k$ can be computed modulo a prime $p$ using the recurrence

for $m \ge 1$, starting from $B_0 = 1$. All arithmetic is in $\mathbb{F}_p$.

Proof. This is a direct rearrangement of the defining recurrence $\sum_{j=0}^m \binom{m+1}{j}B_j = 0$, solving for $B_m$. Division by $m+1$ is valid when $p \nmid m+1$; when $p \mid m+1$, use the Kummer congruences or handle as a special case. $\square$

Editorial

f_k(n) = sum_{i=1}^{n} i^k Compute sums of f_k evaluated at various points modulo primes, using Faulhaber’s formulas and Bernoulli numbers. We iterate over each prime p in primes. We then compute Bernoulli numbers B[0..k] mod p. Finally, evaluate S_k(p).

Pseudocode

for each prime p in primes
Compute Bernoulli numbers B[0..k] mod p
Evaluate S_k(p)

Complexity Analysis

• Time: $O(k^2)$ per prime for computing Bernoulli numbers, plus $O(k)$ per prime for evaluating $S_k(p)$. For $P$ primes total: $O(P \cdot k^2)$. If $k$ is fixed and Bernoulli numbers are precomputed once (when working over $\mathbb{Q}$ and then reducing mod $p$), the per-prime cost drops to $O(k)$, giving $O(k^2 + P \cdot k)$ total.
• Space: $O(k)$ for the Bernoulli number array.

Answer

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

typedef long long ll;

ll mod_pow(ll base, ll exp, ll mod) {
    ll result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

ll mod_inv(ll a, ll p) {
    return mod_pow(a, p - 2, p);
}

// Compute Bernoulli numbers B[0..k] mod p
vector<ll> bernoulli_mod(int k, ll p) {
    vector<ll> B(k + 1, 0);
    B[0] = 1;
    for (int m = 1; m <= k; m++) {
        ll s = 0;
        ll c = 1; // C(m+1, j)
        for (int j = 0; j < m; j++) {
            s = (s + c % p * B[j]) % p;
            c = c % p * ((m + 1 - j) % p) % p * mod_inv(j + 1, p) % p;
        }
        B[m] = (p - s % p) % p * mod_inv(m + 1, p) % p;
    }
    return B;
}

// Compute S_k(p) = sum_{i=1}^{p-1} f_k(i) mod p
// Using Faulhaber + Fermat's Little Theorem
ll sum_power_sums(int k, ll p) {
    auto B = bernoulli_mod(k, p);
    ll inv_k1 = mod_inv(k + 1, p);
    ll result = 0;
    ll c = 1; // C(k+1, j)

    for (int j = 0; j <= k; j++) {
        int m = k + 1 - j;
        ll inner;
        if (m == 0) {
            inner = (p - 1) % p;
        } else if (m % (p - 1) == 0) {
            inner = p - 1; // = -1 mod p
        } else {
            inner = 0;
        }
        result = (result + c % p * (B[j] % p) % p * inner) % p;
        c = c % p * ((k + 1 - j) % p) % p * mod_inv(j + 1, p) % p;
    }

    return result % p * inv_k1 % p;
}

int main() {
    // Demonstration: compute S_k(p) for small primes and k values
    vector<ll> primes = {5, 7, 11, 13, 17, 19, 23};
    for (ll p : primes) {
        for (int k = 1; k <= 4; k++) {
            ll val = sum_power_sums(k, p);
            cout << "S_" << k << "(" << p << ") = " << val << endl;
        }
    }
    return 0;
}

"""
Problem 487: Sums of Power Sums
f_k(n) = sum_{i=1}^{n} i^k
Compute sums of f_k evaluated at various points modulo primes,
using Faulhaber's formulas and Bernoulli numbers.
"""

def bernoulli_mod_p(k: int, p: int) -> list:
    """Compute Bernoulli numbers B_0, B_1, ..., B_k modulo prime p."""
    B = [0] * (k + 1)
    B[0] = 1
    for m in range(1, k + 1):
        # B_m = -1/(m+1) * sum_{j=0}^{m-1} C(m+1,j) * B_j
        s = 0
        c = 1  # binomial coefficient C(m+1, j)
        for j in range(m):
            s = (s + c * B[j]) % p
            c = c * (m + 1 - j) % p * pow(j + 1, p - 2, p) % p
        B[m] = (-(s) * pow(m + 1, p - 2, p)) % p
    return B

def power_sum_mod(k: int, n: int, p: int):
    """Compute f_k(n) = sum_{i=1}^{n} i^k mod p using Faulhaber's formula."""
    B = bernoulli_mod_p(k, p)
    inv_k1 = pow(k + 1, p - 2, p)
    result = 0
    c = 1  # C(k+1, j)
    n_pow = pow(n, k + 1, p)  # n^(k+1-j), starting with j=0
    for j in range(k + 1):
        result = (result + c * B[j] % p * n_pow) % p
        # Update for next j
        c = c * (k + 1 - j) % p * pow(j + 1, p - 2, p) % p
        n_pow = n_pow * pow(n, p - 2, p) % p  # divide by n
    return result * inv_k1 % p

def sum_of_power_sums_mod(k: int, p: int):
    """Compute S_k(p) = sum_{i=1}^{p-1} f_k(i) mod p."""
    total = 0
    for i in range(1, p):
        total = (total + power_sum_mod(k, i, p)) % p
    return total

def sum_of_power_sums_fast(k: int, p: int):
    """
    Compute S_k(p) = sum_{i=1}^{p-1} f_k(i) mod p using Faulhaber + Fermat.
    sum_{i=1}^{p-1} i^m = -1 mod p if (p-1)|m and m>0, else 0 mod p (for m>0).
    """
    B = bernoulli_mod_p(k, p)
    inv_k1 = pow(k + 1, p - 2, p)

    result = 0
    c = 1  # C(k+1, j)
    for j in range(k + 1):
        m = k + 1 - j
        # sum_{i=1}^{p-1} i^m mod p
        if m == 0:
            inner = (p - 1) % p
        elif m % (p - 1) == 0:
            inner = p - 1  # = -1 mod p
        else:
            inner = 0
        result = (result + c * B[j] % p * inner) % p
        c = c * (k + 1 - j) % p * pow(j + 1, p - 2, p) % p

    return result * inv_k1 % p

# Verify: compare brute-force vs fast for small cases
p = 17
for k in range(1, 10):
    bf = sum_of_power_sums_mod(k, p)
    fast = sum_of_power_sums_fast(k, p)
    assert bf == fast, f"Mismatch at k={k}, p={p}: {bf} vs {fast}"

print(f"Verification passed for p={p}, k=1..9")

# Show some values
for p in [5, 7, 11, 13]:
    for k in [1, 2, 3, 4]:
        val = sum_of_power_sums_fast(k, p)
        print(f"S_{k}({p}) = {val}")