Problem 545: Faulhaber's Formulas

Project Euler Problem 545 Solution

EulerSolve provides an optimized solution for Project Euler Problem 545, Faulhaber's Formulas, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For even $k\ge 2$, let $D(k)$ be the denominator of the Bernoulli number $B_k$ written in lowest terms. Define $F(n)$ as the $n$-th positive integer $k$ for which $$D(k)=20010.$$ The goal is to compute $F(10^5)$. The crucial observation is that the problem does not require explicit Bernoulli-number arithmetic: only the prime divisibility pattern of $k$ matters. Mathematical Approach The three implementations turn the denominator condition into a sieve over multipliers of a mandatory base value. Step 1: Apply the Von Staudt-Clausen theorem For every even $k$, the denominator of $B_k$ is $$D(k)=\prod_{p-1\mid k} p,$$ where the product runs over all primes $p$ such that $p-1$ divides $k$. So the denominator is determined entirely by which values $p-1$ occur as divisors of $k$. Once this formula is used, the Bernoulli numbers themselves disappear from the computation. Step 2: Force the primes that must appear The target denominator factors as $$20010=2\cdot 3\cdot 5\cdot 23\cdot 29.$$ Therefore $k$ must be divisible by the corresponding values $$1,\ 2,\ 4,\ 22,\ 28,$$ because those are the numbers $p-1$ attached to the required primes....

Detailed mathematical approach

Problem Summary

For even $k\ge 2$, let $D(k)$ be the denominator of the Bernoulli number $B_k$ written in lowest terms. Define $F(n)$ as the $n$-th positive integer $k$ for which

$$D(k)=20010.$$

The goal is to compute $F(10^5)$. The crucial observation is that the problem does not require explicit Bernoulli-number arithmetic: only the prime divisibility pattern of $k$ matters.

Mathematical Approach

The three implementations turn the denominator condition into a sieve over multipliers of a mandatory base value.

Step 1: Apply the Von Staudt-Clausen theorem

For every even $k$, the denominator of $B_k$ is

$$D(k)=\prod_{p-1\mid k} p,$$

where the product runs over all primes $p$ such that $p-1$ divides $k$. So the denominator is determined entirely by which values $p-1$ occur as divisors of $k$. Once this formula is used, the Bernoulli numbers themselves disappear from the computation.

Step 2: Force the primes that must appear

The target denominator factors as

$$20010=2\cdot 3\cdot 5\cdot 23\cdot 29.$$

Therefore $k$ must be divisible by the corresponding values

$$1,\ 2,\ 4,\ 22,\ 28,$$

because those are the numbers $p-1$ attached to the required primes. The nontrivial conditions combine into

$$\operatorname{lcm}(2,4,22,28)=308.$$

Hence every admissible term has the form

$$k=308m,\qquad m\ge 1.$$

This guarantees that the primes $2,3,5,23,29$ all appear in the denominator. The remaining work is to make sure that no other prime appears.

Step 3: Turn every extra prime into a forbidden divisor of the multiplier

Take a prime $p\notin\{2,3,5,23,29\}$. If $p-1\mid k$, then $p$ would also divide $D(k)$, which is forbidden. Write

$$d=p-1,\qquad g=\gcd(d,308),\qquad d=g\,u.$$

Since $308=g\cdot \frac{308}{g}$, we have

$$d\mid 308m \iff g\,u\mid g\cdot \frac{308}{g}\,m \iff u\mid \frac{308}{g}\,m.$$

Now remove the common factor between $u$ and $308/g$. Define

$$r=\frac{u}{\gcd\left(u,\frac{308}{g}\right)}.$$

Then

$$u\mid \frac{308}{g}\,m \iff r\mid m.$$

So each unwanted prime produces a forbidden divisor $r$: every multiple of $r$ gives an invalid multiplier $m$.

Step 4: Generate forbidden divisors with an arithmetic-form sieve

Every prime $p>2$ is odd, so every value $p-1$ is even. Therefore $g=\gcd(p-1,308)$ can only be one of the even divisors of $308$:

$$g\in\{2,4,14,22,28,44,154,308\}.$$

For a fixed $g$, the candidate primes are exactly the numbers of the form

$$p=g\,u+1.$$

Instead of primality-testing each candidate separately, the implementations sieve directly in the $u$-space. For every small prime $q\nmid g$, the congruence

$$g\,u+1\equiv 0 \pmod q$$

has the unique solution

$$u\equiv -g^{-1}\pmod q,$$

because $g$ has a multiplicative inverse modulo $q$. Every value in that arithmetic progression makes $g\,u+1$ composite, so it can be crossed out. The sieve starts far enough along the progression that the prime $q$ itself is not accidentally removed when $g\,u+1=q$.

Step 5: Count admissible multipliers

After the sieve for a fixed $g$, every surviving $u$ gives a prime $p=g\,u+1$. Each such prime is translated into its forbidden divisor $r$, and all multiples of $r$ are marked in a second sieve over the multiplier axis. If the unmarked multipliers are

$$m_1<m_2<m_3<\dots,$$

then the required sequence is simply

$$F(n)=308m_n.$$

Worked Example

The smallest multiplier is $m=1$, so $k=308$. The divisors of $308$ are

$$1,2,4,7,11,14,22,28,44,77,154,308.$$

Among these, the values $d$ for which $d+1$ is prime are $1,2,4,22,28$, yielding

$$2,3,5,23,29.$$

Therefore

$$D(308)=2\cdot 3\cdot 5\cdot 23\cdot 29=20010,$$

so $308$ is valid. In contrast, if $m=2$, then $k=616$, and $617$ is prime with $617-1=616\mid k$. That extra prime would appear in the denominator, so $616$ is invalid. A second useful example is $p=67$: here $67-1=66$, $\gcd(66,308)=22$, so $u=3$ and $308/22=14$. Thus

$$r=\frac{3}{\gcd(3,14)}=3,$$

which means every multiple of $3$ is invalid. The implementations also verify the checkpoint

$$F(10)=96404.$$

How the Code Works

The C++, Python, and Java implementations share the same structure. They first build a small-prime table up to $\sqrt{308X+1}$, where $X$ is the current search bound on the multipliers. That table supports both the tiny validation checks and the arithmetic-form sieves that generate primes of the form $g\,u+1$.

For each even divisor $g$ of $308$, the implementation allocates a boolean array indexed by $u=1,2,\dots,X$. Using each small prime $q$ that does not divide $g$, it computes the unique residue class modulo $q$ for which $g\,u+1$ is divisible by $q$, then clears that whole progression. Any index that survives corresponds to a prime candidate $g\,u+1$. The required primes $2,3,5,23,29$ are skipped, while every other surviving prime is converted into a forbidden divisor of the multiplier.

A second sieve marks every multiple of every forbidden divisor. Scanning the multiplier array from left to right then produces the admissible values in increasing order. If the requested index has not yet been reached, the search bound is doubled and the process is repeated. The implementations also include small consistency checks such as $D(4)=30$, $D(308)=20010$, and $F(10)=96404$.

Complexity Analysis

Let $X$ be the search bound on the multipliers $m$. Building the prime table up to $\sqrt{308X+1}$ costs $O(\sqrt{X}\log\log X)$ time and $O(\sqrt{X})$ memory. The eight arithmetic-form sieves over $u$ have total harmonic cost $O(X\log\log X)$ up to a constant factor. The final multiple-marking phase costs

$$O\left(X\sum_{r\in R}\frac{1}{r}\right),$$

where $R$ is the set of forbidden divisors discovered up to $X$; in the worst case this is $O(X\log X)$. Therefore the implemented method uses $O(X)$ memory and worst-case time $O(X\log X)$, while still being vastly faster in practice than recomputing Bernoulli denominators for every even $k$.

Footnotes and References

Problem page: https://projecteuler.net/problem=545
Von Staudt-Clausen theorem: Wikipedia - Von Staudt-Clausen theorem
Bernoulli numbers: Wikipedia - Bernoulli number
Faulhaber's formula: Wikipedia - Faulhaber's formula
Modular multiplicative inverse: Wikipedia - Modular multiplicative inverse

Problem 545 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <vector>
#include <cmath>
#include <functional>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

static constexpr u64 kDenTarget = 20010ULL;
static constexpr u64 kL = 308ULL;  // lcm(4,22,28)

static std::vector<int> sieve_primes(int n) {
    std::vector<bool> is_comp(static_cast<std::size_t>(n + 1), false);
    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (!is_comp[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
            if (i <= n / i) {
                for (int j = i * i; j <= n; j += i) is_comp[static_cast<std::size_t>(j)] = true;
            }
        }
    }
    return primes;
}

static bool is_prime_u64(u64 x, const std::vector<int>& primes) {
    if (x < 2) return false;
    for (int p : primes) {
        const u64 pp = static_cast<u64>(p);
        if (pp * pp > x) break;
        if (x % pp == 0) return x == pp;
    }
    return true;
}

static int mod_inv(int a, int mod) {
    // Extended Euclid, mod is prime in our use.
    int t = 0, newt = 1;
    int r = mod, newr = a;
    while (newr != 0) {
        const int q = r / newr;
        const int tmp_t = t - q * newt;
        t = newt;
        newt = tmp_t;
        const int tmp_r = r - q * newr;
        r = newr;
        newr = tmp_r;
    }
    if (r != 1) return 0;
    if (t < 0) t += mod;
    return t;
}

static u64 bernoulli_denominator_even_k(u64 k, const std::vector<int>& primes_small) {
    // Von Staudt–Clausen: denom(B_k) = Π_{p-1 | k} p for even k>=2.
    // Compute by enumerating divisors d|k and checking if d+1 is prime.
    std::vector<std::pair<u64, int>> fac;
    u64 x = k;
    for (int p : primes_small) {
        const u64 pp = static_cast<u64>(p);
        if (pp * pp > x) break;
        if (x % pp == 0) {
            int e = 0;
            while (x % pp == 0) {
                x /= pp;
                ++e;
            }
            fac.emplace_back(pp, e);
        }
    }
    if (x > 1) fac.emplace_back(x, 1);

    u64 den = 1;
    std::vector<u64> divs{1};
    for (const auto& [p, e] : fac) {
        const std::size_t cur = divs.size();
        u64 pe = 1;
        for (int i = 1; i <= e; ++i) {
            pe *= p;
            for (std::size_t j = 0; j < cur; ++j) divs.push_back(divs[j] * pe);
        }
    }
    for (u64 d : divs) {
        const u64 cand = d + 1;
        if (is_prime_u64(cand, primes_small)) den *= cand;
    }
    return den;
}

static u64 brute_F_10() {
    // Brute the first 10 k with D(k)=20010 by scanning multiples of 308.
    const int lim = 200000;  // enough to reach F(10)=96404.
    const auto primes = sieve_primes(static_cast<int>(std::sqrt(lim + 1)) + 5);

    int found = 0;
    u64 tenth = 0;
    for (u64 m = 1; kL * m <= static_cast<u64>(lim); ++m) {
        const u64 k = kL * m;
        const u64 den = bernoulli_denominator_even_k(k, primes);
        if (den == kDenTarget) {
            ++found;
            if (found == 10) {
                tenth = k;
                break;
            }
        }
    }
    return tenth;
}

static u64 kth_multiplier_with_den(u64 M, u64 kth) {
    // m is valid iff no prime p outside {2,3,5,23,29} has (p-1) | (308*m).
    // For each such p, let r=(p-1)/gcd(p-1,308). Then r | m implies invalid.
    //
    // We compute all forbidden r<=M by sieving primes of the form p=g*u+1 with g|308.

    const u64 pmax = kL * M + 1ULL;
    const int qmax = static_cast<int>(std::sqrt(static_cast<long double>(pmax))) + 2;
    const std::vector<int> primes = sieve_primes(qmax);

    std::vector<char> forbidden(static_cast<std::size_t>(M + 1), 0);

    auto is_target_prime = [](u64 p) -> bool {
        return p == 2ULL || p == 3ULL || p == 5ULL || p == 23ULL || p == 29ULL;
    };

    const int gs[] = {2, 4, 14, 22, 28, 44, 154, 308};
    for (int g : gs) {
        const u64 maxp_g = static_cast<u64>(g) * M + 1ULL;
        const u64 qlim = static_cast<u64>(std::sqrt(static_cast<long double>(maxp_g))) + 1ULL;

        std::vector<char> is_prime_u(static_cast<std::size_t>(M + 1), 1);
        is_prime_u[0] = 0;

        for (int q : primes) {
            if (static_cast<u64>(q) > qlim) break;
            if (g % q == 0) continue;

            const int a = g % q;
            const int inv = mod_inv(a, q);
            // u ≡ (-1) * g^{-1} (mod q)
            const int u0 = static_cast<int>((static_cast<long long>(q - 1) * inv) % q);

            const u64 qq = static_cast<u64>(q);
            const u64 u_start = (qq * qq - 1ULL + static_cast<u64>(g) - 1ULL) / static_cast<u64>(g);
            u64 u = static_cast<u64>(u0);
            if (u < u_start) {
                u += ((u_start - u + qq - 1ULL) / qq) * qq;
            }
            for (; u <= M; u += qq) is_prime_u[static_cast<std::size_t>(u)] = 0;
        }

        const u64 div = kL / static_cast<u64>(g);
        for (u64 u = 1; u <= M; ++u) {
            if (!is_prime_u[static_cast<std::size_t>(u)]) continue;
            const u64 p = static_cast<u64>(g) * u + 1ULL;
            if (is_target_prime(p)) continue;
            const u64 r = u / std::gcd(u, div);
            forbidden[static_cast<std::size_t>(r)] = 1;
        }
    }

    std::vector<char> invalid(static_cast<std::size_t>(M + 1), 0);
    for (u64 r = 2; r <= M; ++r) {
        if (!forbidden[static_cast<std::size_t>(r)]) continue;
        for (u64 x = r; x <= M; x += r) invalid[static_cast<std::size_t>(x)] = 1;
    }

    u64 cnt = 0;
    for (u64 m = 1; m <= M; ++m) {
        if (invalid[static_cast<std::size_t>(m)]) continue;
        ++cnt;
        if (cnt == kth) return m;
    }
    return 0ULL;
}

static u64 solve_F_1e5() {
    const u64 kth = 100'000ULL;
    u64 M = 4'000'000ULL;
    while (true) {
        const u64 m = kth_multiplier_with_den(M, kth);
        if (m != 0) return kL * m;
        M *= 2;
    }
}

}  // namespace

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    // Statement validations.
    {
        const auto primes = sieve_primes(1000);
        assert(bernoulli_denominator_even_k(4, primes) == 30ULL);
        assert(bernoulli_denominator_even_k(308, primes) == kDenTarget);
    }
    assert(brute_F_10() == 96404ULL);

    std::cout << solve_F_1e5() << '\n';
    return 0;
}

Python

import math

kDenTarget = 20010
kL = 308

def sieve_primes(n):
    is_comp = [False] * (n + 1)
    primes = []
    for i in range(2, n + 1):
        if not is_comp[i]:
            primes.append(i)
            if i <= n // i:
                for j in range(i * i, n + 1, i):
                    is_comp[j] = True
    return primes

def is_prime_u64(x, primes):
    if x < 2: return False
    for p in primes:
        if p * p > x: break
        if x % p == 0: return x == p
    return True

def mod_inv(a, mod):
    return pow(a, -1, mod)

def bernoulli_denominator_even_k(k, primes_small):
    fac = []
    x = k
    for p in primes_small:
        if p * p > x: break
        if x % p == 0:
            e = 0
            while x % p == 0:
                x //= p
                e += 1
            fac.append((p, e))
    if x > 1: fac.append((x, 1))
    
    divs = [1]
    for p, e in fac:
        cur_len = len(divs)
        pe = 1
        for i in range(1, e + 1):
            pe *= p
            for j in range(cur_len):
                divs.append(divs[j] * pe)
                
    den = 1
    for d in divs:
        cand = d + 1
        if is_prime_u64(cand, primes_small):
            den *= cand
    return den

def brute_F_10():
    lim = 200000
    primes = sieve_primes(int(math.sqrt(lim + 1)) + 5)
    found = 0
    m = 1
    while kL * m <= lim:
        k = kL * m
        den = bernoulli_denominator_even_k(k, primes)
        if den == kDenTarget:
            found += 1
            if found == 10:
                return k
        m += 1
    return 0

def kth_multiplier_with_den(M, kth):
    pmax = kL * M + 1
    qmax = int(math.sqrt(pmax)) + 2
    primes = sieve_primes(qmax)
    
    forbidden = bytearray(M + 1)
    
    gs = [2, 4, 14, 22, 28, 44, 154, 308]
    for g in gs:
        maxp_g = g * M + 1
        qlim = int(math.sqrt(maxp_g)) + 1
        
        is_prime_u = bytearray(b'\x01' * (M + 1))
        is_prime_u[0] = 0
        
        for q in primes:
            if q > qlim: break
            if g % q == 0: continue
            
            a = g % q
            inv = mod_inv(a, q)
            u0 = ((q - 1) * inv) % q
            
            u_start = (q * q - 1 + g - 1) // g
            u = u0
            if u < u_start:
                u += ((u_start - u + q - 1) // q) * q
            
            if u <= M:
                end_range = M + 1
                count = (end_range - u + q - 1) // q
                is_prime_u[u:end_range:q] = b'\x00' * count
                
        div = kL // g
        
        for u in range(1, M + 1):
            if not is_prime_u[u]: continue
            p = g * u + 1
            if p in (2, 3, 5, 23, 29): continue
            r = u // math.gcd(u, div)
            if r <= M:
                forbidden[r] = 1
                
    invalid = bytearray(M + 1)
    for r in range(2, M + 1):
        if not forbidden[r]: continue
        end_range = M + 1
        count = (end_range - r + r - 1) // r
        invalid[r:end_range:r] = b'\x01' * count
            
    cnt = 0
    for m in range(1, M + 1):
        if not invalid[m]:
            cnt += 1
            if cnt == kth:
                return m
    return 0
    
def solve_F_1e5():
    kth = 100000
    M = 4000000
    while True:
        m = kth_multiplier_with_den(M, kth)
        if m != 0: return kL * m
        M *= 2

def solve():
    return str(solve_F_1e5())

if __name__ == '__main__':
    print(solve())

Java

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

public class Euler545 {
    static final long kDenTarget = 20010L;
    static final long kL = 308L;

    static List<Integer> sievePrimes(int n) {
        boolean[] isComp = new boolean[n + 1];
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; i++) {
            if (!isComp[i]) {
                primes.add(i);
                if (i <= n / i) {
                    for (int j = i * i; j <= n; j += i) {
                        isComp[j] = true;
                    }
                }
            }
        }
        return primes;
    }

    static boolean isPrimeU64(long x, List<Integer> primes) {
        if (x < 2)
            return false;
        for (int p : primes) {
            long pp = p;
            if (pp * pp > x)
                break;
            if (x % pp == 0)
                return x == pp;
        }
        return true;
    }

    static int modInv(int a, int mod) {
        int t = 0, newt = 1;
        int r = mod, newr = a;
        while (newr != 0) {
            int q = r / newr;
            int tmpT = t - q * newt;
            t = newt;
            newt = tmpT;
            int tmpR = r - q * newr;
            r = newr;
            newr = tmpR;
        }
        if (r != 1)
            return 0;
        if (t < 0)
            t += mod;
        return t;
    }

    static long bernoulliDenominatorEvenK(long k, List<Integer> primesSmall) {
        List<long[]> fac = new ArrayList<>();
        long x = k;
        for (int p : primesSmall) {
            long pp = p;
            if (pp * pp > x)
                break;
            if (x % pp == 0) {
                long e = 0;
                while (x % pp == 0) {
                    x /= pp;
                    e++;
                }
                fac.add(new long[] { pp, e });
            }
        }
        if (x > 1)
            fac.add(new long[] { x, 1 });

        List<Long> divs = new ArrayList<>();
        divs.add(1L);
        for (long[] pair : fac) {
            long p = pair[0];
            long e = pair[1];
            int cur = divs.size();
            long pe = 1;
            for (int i = 1; i <= e; i++) {
                pe *= p;
                for (int j = 0; j < cur; j++) {
                    divs.add(divs.get(j) * pe);
                }
            }
        }

        long den = 1;
        for (long d : divs) {
            long cand = d + 1;
            if (isPrimeU64(cand, primesSmall)) {
                den *= cand;
            }
        }
        return den;
    }

    static long bruteF10() {
        int lim = 200000;
        List<Integer> primes = sievePrimes((int) Math.sqrt(lim + 1) + 5);
        int found = 0;
        for (long m = 1; kL * m <= lim; m++) {
            long k = kL * m;
            long den = bernoulliDenominatorEvenK(k, primes);
            if (den == kDenTarget) {
                found++;
                if (found == 10)
                    return k;
            }
        }
        return 0;
    }

    static long gcd(long a, long b) {
        while (b != 0) {
            long t = a % b;
            a = b;
            b = t;
        }
        return a;
    }

    static long kthMultiplierWithDen(int M, int kth) {
        long pmax = kL * M + 1L;
        int qmax = (int) Math.sqrt(pmax) + 2;
        List<Integer> primes = sievePrimes(qmax);

        boolean[] forbidden = new boolean[M + 1];

        int[] gs = { 2, 4, 14, 22, 28, 44, 154, 308 };
        for (int g : gs) {
            long maxpG = (long) g * M + 1L;
            long qlim = (long) Math.sqrt(maxpG) + 1L;

            boolean[] isPrimeU = new boolean[M + 1];
            Arrays.fill(isPrimeU, true);
            isPrimeU[0] = false;

            for (int q : primes) {
                if (q > qlim)
                    break;
                if (g % q == 0)
                    continue;

                int a = g % q;
                int inv = modInv(a, q);
                int u0 = (int) (((long) (q - 1) * inv) % q);

                long uStart = ((long) q * q - 1L + g - 1L) / g;
                long u = u0;
                if (u < uStart) {
                    u += ((uStart - u + q - 1L) / q) * q;
                }

                if (u <= M) {
                    int step = q;
                    for (int ui = (int) u; ui <= M; ui += step) {
                        isPrimeU[ui] = false;
                    }
                }
            }

            long div = kL / g;
            for (int u = 1; u <= M; u++) {
                if (!isPrimeU[u])
                    continue;
                long p = (long) g * u + 1L;
                if (p == 2 || p == 3 || p == 5 || p == 23 || p == 29)
                    continue;
                long r = u / gcd(u, div);
                if (r <= M) {
                    forbidden[(int) r] = true;
                }
            }
        }

        boolean[] invalid = new boolean[M + 1];
        for (int r = 2; r <= M; r++) {
            if (!forbidden[r])
                continue;
            for (int x = r; x <= M; x += r) {
                invalid[x] = true;
            }
        }

        int cnt = 0;
        for (int m = 1; m <= M; m++) {
            if (!invalid[m]) {
                cnt++;
                if (cnt == kth)
                    return m;
            }
        }
        return 0;
    }

    static long solveF1e5() {
        int kth = 100000;
        int M = 4000000;
        while (true) {
            long m = kthMultiplierWithDen(M, kth);
            if (m != 0)
                return kL * m;
            M *= 2;
        }
    }

    public static String solve() {
        return Long.toString(solveF1e5());
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}