Project Euler Problem 545
The sum of the kth powers of the first n positive integers can be expressed as a polynomial of degree k+1 with rational
Solution
Answer: 921107572
Let
$$S_k(n)=\sum_{i=1}^n i^k$$
be the Faulhaber polynomial. We are interested in the coefficient of $n$, whose denominator is $D(k)$.
1. Bernoulli-number connection
Faulhaber’s formula is
$$\sum_{i=1}^n i^k = \frac{1}{k+1} \sum_{j=0}^{k} (-1)^j \binom{k+1}{j} B_j n^{k+1-j},$$
where $B_j$ are the Bernoulli numbers.
The coefficient of $n$ comes from the term $j=k$, so:
$$a_1 = B_k.$$
Hence:
$$D(k)=\text{denominator}(B_k).$$
For example,
$$B_4=-\frac{1}{30},$$
so $D(4)=30$, matching the problem statement.
2. Apply the von Staudt–Clausen theorem
The theorem says:
$$\operatorname{den}(B_k) = \prod_{p-1\mid k} p \qquad (k \text{ even}),$$
where the product runs over primes $p$ such that $p-1\mid k$.
We want:
$$D(k)=20010.$$
Factor:
$$20010 = 2\cdot 3\cdot 5\cdot 23\cdot 29.$$
Thus the only primes $p$ with $p-1\mid k$ must be
$$p\in{2,3,5,23,29}.$$
Equivalently,
$$1,2,4,22,28 \mid k,$$
but no other $p-1$ may divide $k$.
The least common multiple is
$$\mathrm{lcm}(1,2,4,22,28)=308.$$
So every valid $k$ has the form
$$k=308t.$$
Now the condition becomes:
For every divisor $d\mid 308$,
$$dt+1$$
must be composite (otherwise a new prime $q=dt+1$ would satisfy $q-1\mid k$, adding an unwanted factor to the denominator).
This gives a recursive multiplicative test:
- $t=1$ is valid ($k=308$, giving $F(1)=308$).
- A number $t$ is valid iff:
- all $dt+1$ for $d\mid 308$ are composite, and
- every $t/p$ (for prime divisors $p\mid t$) is valid.
This reproduces:
$$F(10)=96404,$$
matching the statement.
3. Python implementation
import sympy as sp
DIVS = sp.divisors(308)
def find_F(target):
N = 200000
while True:
# Smallest prime factor sieve
spf = list(range(N + 1))
for i in range(2, int(N**0.5) + 1):
if spf[i] == i:
for j in range(i * i, N + 1, i):
if spf[j] == j:
spf[j] = i
good = [False] * (N + 1)
good[1] = True
count = 1 # t = 1 corresponds to F(1)=308
for t in range(2, N + 1):
# Local primality condition:
# d*t + 1 must be composite for all d | 308
local_ok = True
for d in DIVS:
if sp.isprime(d * t + 1):
local_ok = False
break
if not local_ok:
continue
# Recursive divisor condition
n = t
ok = True
used = set()
while n > 1:
p = spf[n]
if p not in used:
used.add(p)
if not good[t // p]:
ok = False
break
n //= p
if ok:
good[t] = True
count += 1
if count == target:
return 308 * t
N *= 2
print(find_F(10)) # 96404
print(find_F(100000)) # final answer
4. Code walkthrough
DIVS = divisors(308)computes all divisors of $308$.- We search for numbers $t$ such that $k=308t$.
- For each $t$, we test whether every $dt+1$ is composite.
- We also enforce the recursive divisor property by checking all $t/p$, where $p$ runs over prime divisors of $t$.
- We count valid $t$'s in increasing order until the $100000$-th one.
Verification:
- $F(1)=308$ ✔
- $F(10)=96404$ ✔
The computation gives:
$$t_{100000}=2,!990,!609$$
and therefore
$$F(10^5)=308\times 2,!990,!609 =921,!107,!572.$$
Answer: 921107572