Project Euler Problem 971

Solution

Answer: 33626723890930

Let

$$p=5k-4 \qquad\Longrightarrow\qquad k=\frac{p+4}{5},$$

$$k-1=\frac{p-1}{5}.$$

Define

$$d=\frac{p-1}{5}.$$

For nonzero $x\in\mathbb F_p$,

$$f_p(x)=x^k+x=x\left(x^{k-1}+1\right)=x\left(x^d+1\right).$$

Now $x^d$ is always a 5th root of unity in $\mathbb F_p$, since

$$(x^d)^5=x^{p-1}=1.$$

Thus the nonzero elements split into the five cosets determined by the value of $x^d$.

For a fixed 5th root $r$,

$$x^d=r$$

implies

$$f_p(x)=x(1+r).$$

Hence the induced action on the five cosets is completely determined by the map

$$r \mapsto r(1+r)^d.$$

Every element in a periodic coset is periodic, and every element in a non-periodic coset eventually leaves it forever. Therefore:

each periodic coset contributes exactly $d=(p-1)/5$ points,
$x=0$ is always fixed.

So if $c$ is the number of periodic cosets among the five nonzero cosets,

$$C(p)=1+c\cdot \frac{p-1}{5}.$$

The following Python program computes this efficiently.

import math

def S(N):
    # odd-only sieve
    sieve = bytearray(b'\x01') * (N // 2)

    limit = int(math.isqrt(N))
    for i in range(3, limit + 1, 2):
        if sieve[i // 2]:
            start = i * i // 2
            sieve[start::i] = b'\x00' * (((N // 2 - 1 - start) // i) + 1)

    total = 0

    # only primes p == 1 (mod 5) matter
    for p in range(11, N + 1, 10):
        if not sieve[p // 2]:
            continue

        d = (p - 1) // 5

        # find a primitive 5th root of unity
        a = 2
        while True:
            t = pow(a, d, p)
            if t != 1:
                break
            a += 1

        t2 = (t * t) % p
        t3 = (t2 * t) % p
        t4 = (t2 * t2) % p

        roots = (1, t, t2, t3, t4)
        idx = {roots[i]: i for i in range(5)}

        # induced map on the 5 cosets
        T = [0] * 5
        for i, r in enumerate(roots):
            nxt = (r * pow(1 + r, d, p)) % p
            T[i] = idx[nxt]

        # count periodic states in the 5-state graph
        periodic = 0
        for i in range(5):
            j = i
            seen = 0
            while not ((seen >> j) & 1):
                seen |= 1 << j
                j = T[j]
            if j == i:
                periodic += 1

        total += 1 + periodic * d

    return total

print(S(100))      # 127
print(S(10**8))

The program reproduces the example:

$$S(100)=127.$$

Running it for $10^8$ gives:

$$S(10^8)=33626723890930.$$

Answer: 33626723890930