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

Project Euler Problem 545

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:
  1. all $dt+1$ for $d\mid 308$ are composite, and
  2. 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