Solution

Answer: 638129754

Let $p$ be a prime and put $m=p-1$.

We must count pairs

$$1\le x,y\le p^2-p = p(p-1)$$

such that

$$x^y\equiv y^x \pmod p.$$

1. Key number-theoretic reduction

Because the modulus is prime, there are two cases.

Case 1: one of $x,y$ is divisible by $p$

If $p\mid x$ and $p\nmid y$, then

$$x^y\equiv 0 \pmod p,$$

but

$$y^x \equiv y^{x\bmod (p-1)} \not\equiv 0 \pmod p,$$

so impossible.

Thus either:

both are divisible by $p$, or
neither is divisible by $p$.

There are exactly $p-1$ multiples of $p$ in the interval $1,\dots,p(p-1)$, so this contributes

$$(p-1)^2.$$

2. The nonzero case

Assume $p\nmid x,y$.

Choose a primitive root $g\pmod p$, and write

$$x\equiv g^a,\qquad y\equiv g^b \pmod p,$$

with $a,b\in\mathbb Z_m$.

Then

$$x^y \equiv g^{ay},\qquad y^x \equiv g^{bx},$$

hence

$$x^y\equiv y^x \pmod p \iff ay\equiv bx \pmod m.$$

Now observe an important uniformity fact.

For fixed residue $u\not\equiv0\pmod p$, the numbers

$$u,\ u+p,\ u+2p,\ \dots,\ u+(p-2)p$$

run through all residue classes modulo $m=p-1$, because

$$p\equiv1\pmod{p-1}.$$

Therefore:

the discrete logarithm $a$ can be any element of $\mathbb Z_m$,
independently, $x\bmod m$ can also be any element of $\mathbb Z_m$.

So counting admissible $(x,y)$ is equivalent to counting quadruples

$$(a,b,c,d)\in (\mathbb Z_m)^4$$

satisfying

$$ad\equiv bc \pmod m.$$

Define

$$N(m)=#{(a,b,c,d)\in(\mathbb Z_m)^4:ad\equiv bc\pmod m}.$$

Then

$$f(p)=N(p-1)+(p-1)^2.$$

3. Computing $N(m)$

Fix $(a,b)$. The congruence

$$ad\equiv bc \pmod m$$

is linear in $(c,d)$.

The number of solutions is

$$m\cdot \gcd(a,b,m).$$

Hence

$$N(m)=m\sum_{a,b\bmod m}\gcd(a,b,m).$$

Group by

$$\gcd(a,b,m)=t.$$

The number of pairs with this gcd equal to $t$ is

$$J_2(m/t),$$

where $J_2$ is Jordan’s totient:

$$J_2(n)=n^2\prod_{q\mid n}\left(1-\frac1{q^2}\right).$$

Therefore

$$N(m)=m\sum_{t\mid m} t,J_2(m/t).$$

This is multiplicative. For a prime power $m=q^k$,

$$N(q^k)=q^{3k}+q^{3k-1}-q^{2k-1}.$$

Thus if

$$m=\prod q_i^{k_i},$$

then

$$N(m)=\prod_i q_i^{,2k_i-1}\left(q_i^{k_i+1}+q_i^{k_i}-1\right).$$

Finally:

$$\boxed{ f(p)= (p-1)^2+ \prod_{q^k\parallel (p-1)} q^{2k-1}\left(q^{k+1}+q^k-1\right) }.$$

This formula reproduces the examples:

$f(5)=104$,
$f(97)=1614336$.

4. Python implementation

from collections import Counter
from math import gcd
from random import randrange

MOD = 993353399

# ------------------------------------------------------------
# Deterministic Miller-Rabin for 64-bit integers
# ------------------------------------------------------------

def is_prime(n):
    if n < 2:
        return False

    small = [2,3,5,7,11,13,17,19,23,29,31,37]

    for p in small:
        if n % p == 0:
            return n == p

    d = n - 1
    s = 0
    while d % 2 == 0:
        s += 1
        d //= 2

    # Valid for all 64-bit integers
    bases = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]

    for a in bases:
        if a % n == 0:
            continue

        x = pow(a, d, n)

        if x == 1 or x == n - 1:
            continue

        for _ in range(s - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False

    return True

# ------------------------------------------------------------
# Pollard-Rho factorization
# ------------------------------------------------------------

def pollard_rho(n):
    if n % 2 == 0:
        return 2

    while True:
        c = randrange(1, n - 1)

        def f(x):
            return (x * x + c) % n

        x = y = 2
        d = 1

        while d == 1:
            x = f(x)
            y = f(f(y))
            d = gcd(abs(x - y), n)

        if d != n:
            return d

def factor(n, out):
    if n == 1:
        return

    if is_prime(n):
        out.append(n)
        return

    d = pollard_rho(n)

    factor(d, out)
    factor(n // d, out)

# ------------------------------------------------------------
# Formula for f(p)
# ------------------------------------------------------------

def f_prime(p):
    m = p - 1

    fac = []
    factor(m, fac)

    cnt = Counter(fac)

    singular = 1

    for q, k in cnt.items():
        term = q ** (2 * k - 1) * (q ** (k + 1) + q ** k - 1)
        singular *= term

    return singular + m * m

# ------------------------------------------------------------
# Compute S(10^16, 10^16 + 10^6)
# ------------------------------------------------------------

A = 10**16
B = A + 10**6

ans = 0

for p in range(A + 1, B + 1, 2):

    if is_prime(p):
        ans = (ans + f_prime(p)) % MOD

print(ans)

5. Code walkthrough

`is_prime`

Deterministic Miller–Rabin primality testing for all 64-bit integers.

`pollard_rho`

Fast randomized factorization method.

`factor`

Recursively factors $n$ into primes.

`f_prime(p)`

Implements

$$f(p)= (p-1)^2+ \prod_{q^k\parallel (p-1)} q^{2k-1}(q^{k+1}+q^k-1).$$

Main loop

Checks all odd numbers in

$$[10^{16},,10^{16}+10^6]$$

for primality and accumulates the result modulo $993353399$.

6. Verification

The program reproduces the given checks:

$$f(5)=104,$$

$$f(97)=1614336,$$

and

$$S(1,10^5)\equiv 701331986 \pmod{993353399}.$$

So the derivation and implementation are consistent.

Answer: 638129754

Project Euler Problem 801