Project Euler Problem 512

Project Euler Problem 512

Solution

Answer: 50660591862310323

Let

$$f(n)=\left(\sum_{i=1}^{n}\varphi(n^i)\right)\bmod(n+1)$$

and

$$g(n)=\sum_{i=1}^{n} f(i).$$

We need $g(5\times 10^8)$.

1. Mathematical analysis

Step 1: Simplify $f(n)$

Recall the totient identity:

$$\varphi(n^k)=n^{k-1}\varphi(n)$$

for any integer $k\ge1$.

Thus

$$\sum_{i=1}^{n}\varphi(n^i) = \varphi(n)\sum_{i=0}^{n-1} n^i.$$

We only care modulo $n+1$. Since

$$n\equiv -1 \pmod{n+1},$$

we get

$$n^i \equiv (-1)^i \pmod{n+1}.$$

Hence

$$\sum_{i=0}^{n-1}n^i \equiv \sum_{i=0}^{n-1}(-1)^i \pmod{n+1}.$$

Now distinguish parity:

• If $n$ is even, the alternating sum has equally many $+1$ and $-1$, so:

$$\sum_{i=0}^{n-1}(-1)^i=0.$$

Thus

$$f(n)=0.$$

• If $n$ is odd, the alternating sum equals $1$, so:

$$f(n)\equiv \varphi(n)\pmod{n+1}.$$

Since for $n>1$,

$$0<\varphi(n)<n+1,$$

the residue is simply

$$f(n)=\varphi(n).$$

Also $f(1)=1=\varphi(1)$.

Therefore:

$$f(n)= \begin{cases} \varphi(n), & n\text{ odd},\[4pt] 0, & n\text{ even}. \end{cases}$$

So

$$g(N)=\sum_{\substack{n\le N\ n\text{ odd}}}\varphi(n).$$

---

Step 2: Convert to a summatory formula

Use the Möbius inversion formula:

$$\varphi(n) = n\sum_{d\mid n}\frac{\mu(d)}{d}.$$

Since $n$ is odd, all divisors are odd:

$$g(N) = \sum_{\substack{n\le N\ n\text{ odd}}} n \sum_{d\mid n}\frac{\mu(d)}{d}.$$

Write $n=dm$:

$$g(N) = \sum_{\substack{d\le N\ d\text{ odd}}} \mu(d) \sum_{\substack{m\le N/d\ m\text{ odd}}} m.$$

The sum of odd integers up to $M$:

$$1+3+\cdots+(2k-1)=k^2,$$

where

$$k=\left\lfloor \frac{M+1}{2}\right\rfloor.$$

Hence

$$g(N) = \sum_{\substack{d\le N\ d\text{ odd}}} \mu(d) \left( \left\lfloor \frac{\lfloor N/d\rfloor+1}{2} \right\rfloor \right)^2.$$

This can be computed in $O(\sqrt N)$ using divisor grouping and a memoized summatory Möbius function.

---

2. Python implementation

from functools import lru_cache

def compute_g(N):
    # Sieve limit for fast Möbius prefix values
    LIMIT = int(N ** (2 / 3)) + 10

    # Linear sieve for Möbius function
    mu = [0] * (LIMIT + 1)
    mu[1] = 1
    primes = []
    is_composite = [False] * (LIMIT + 1)

    for i in range(2, LIMIT + 1):
        if not is_composite[i]:
            primes.append(i)
            mu[i] = -1

        for p in primes:
            v = i * p
            if v > LIMIT:
                break

            is_composite[v] = True

            if i % p == 0:
                mu[v] = 0
                break
            else:
                mu[v] = -mu[i]

    # Prefix sum of Möbius
    mobius_prefix = [0] * (LIMIT + 1)
    for i in range(1, LIMIT + 1):
        mobius_prefix[i] = mobius_prefix[i - 1] + mu[i]

    @lru_cache(None)
    def mertens(n):
        """M(n) = sum_{k<=n} mu(k)"""
        if n <= LIMIT:
            return mobius_prefix[n]

        result = 1
        left = 2

        while left <= n:
            q = n // left
            right = n // q

            result -= (right - left + 1) * mertens(q)
            left = right + 1

        return result

    @lru_cache(None)
    def odd_mertens(n):
        """Sum of mu(k) over odd k <= n"""
        total = 0

        while n > 0:
            total += mertens(n)
            n //= 2

        return total

    answer = 0
    left = 1

    while left <= N:
        q = N // left
        right = N // q

        # Number of odd integers <= q
        k = (q + 1) // 2

        mobius_sum = odd_mertens(right) - odd_mertens(left - 1)
        answer += k * k * mobius_sum

        left = right + 1

    return answer

# Verification
print(compute_g(100))  # 2007

# Final result
print(compute_g(5 * 10**8))

---

3. Code walkthrough

Möbius sieve

We compute $\mu(n)$ up to about $N^{2/3}$ using a linear sieve.

This gives fast access to:

$$M(n)=\sum_{k\le n}\mu(k)$$

for smaller $n$.

mertens(n)

Uses the classic divide-and-conquer recurrence:

$$M(n) = 1 - \sum_{l=2}^{n} (r-l+1)M!\left(\left\lfloor \frac{n}{l}\right\rfloor\right)$$

grouping equal quotient ranges.

odd_mertens(n)

We only need odd $d$. Observe:

$$M(n)=O(n)-O(n/2),$$

where $O(n)$ is the odd-only Möbius sum.

So:

$$O(n)=M(n)+O(n/2).$$

This becomes:

$$O(n)=M(n)+M(n/2)+M(n/4)+\cdots.$$

Final summation

We group intervals where

$$\left\lfloor \frac{N}{d}\right\rfloor$$

is constant, reducing complexity to roughly $O(\sqrt N)$.

Verification

The code reproduces the given check:

$$g(100)=2007.$$

So the derivation is consistent.

Final result

Running the computation for

$$N=5\times 10^8$$

gives:

Answer: 50660591862310323