Problem 545: Faulhaber's Formulas
View on Project EulerProject 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