Problem 801: $x^y \equiv y^x \pmod{p}$

Mathematical Foundation

Theorem 1 (Discrete Logarithm Reduction). Let $p$ be an odd prime and $g$ a primitive root modulo $p$. For $x, y \in \{1, \ldots, p-1\}$, write $x = g^a$ and $y = g^b$ where $a = \operatorname{ind}_g(x)$ and $b = \operatorname{ind}_g(y)$. Then

Proof. Since $g$ is a primitive root, the map $k \mapsto g^k$ is a bijection $\mathbb{Z}/(p-1)\mathbb{Z} \to (\mathbb{Z}/p\mathbb{Z})^*$. We have $x^y = g^{ay}$ and $y^x = g^{bx}$, so $x^y \equiv y^x \pmod{p}$ if and only if $g^{ay} \equiv g^{bx} \pmod{p}$, which holds if and only if $ay \equiv bx \pmod{p-1}$ by the order of $g$. $\square$

Lemma 1 (Collision via Hash Function). Define $h : \{x \in \{1, \ldots, p-1\} : \gcd(x, p-1) = 1\} \to \mathbb{Z}/(p-1)\mathbb{Z}$ by

Then for $x, y$ both in the domain of $h$, we have $x^y \equiv y^x \pmod{p}$ if and only if $h(x) = h(y)$.

Proof. The condition $ay \equiv bx \pmod{p-1}$ can be rewritten (when $\gcd(x, p-1) = 1$ and $\gcd(y, p-1) = 1$) as $a \cdot x^{-1} \equiv b \cdot y^{-1} \pmod{p-1}$, i.e., $h(x) = h(y)$. $\square$

Lemma 2 (Handling $\gcd(x, p-1) > 1$). When $\gcd(x, p-1) > 1$, the inverse $x^{-1} \pmod{p-1}$ does not exist. For such $x$, the congruence $ay \equiv bx \pmod{p-1}$ must be checked directly: it holds iff $\gcd(x, p-1) \mid (ay - bx)$, which requires case-by-case enumeration grouped by the residue class of $a/\gcd(a, p-1)$.

Proof. This follows from the solvability condition for linear congruences: $ay \equiv bx \pmod{p-1}$ has solutions iff $\gcd(y, p-1) \mid bx - ay$, but here $x, y$ are both given. The condition is simply $p - 1 \mid ay - bx$. $\square$

Editorial

For a prime $p$, let $f(p) = \sum (x+y)$ where the sum is over all pairs $(x,y)$ with $1 \le x < y < p$ such that $x^y \equiv y^x \pmod{p}$. Find $\sum_{p \le 10^6} f(p) \bmod 1\,000\,000\,007$. We iterate over each prime p in primes. We then compute discrete logarithm table via BSGS or direct power. Finally, build hash map: key = h(x), value = list of x.

Pseudocode

for each prime p in primes
Compute discrete logarithm table via BSGS or direct power
Build hash map: key = h(x), value = list of x
Sum over collision pairs
for each bucket B in buckets
sum of (x+y) over pairs = (|B|-1)*S
Handle gcd(x, p-1) > 1 cases separately (direct check)

Complexity Analysis

• Time: For each prime $p$, computing the discrete logarithm table takes $O(p)$. Building the hash map and counting collisions is $O(p)$. Summing over all primes $p \le N$: $O\!\left(\sum_{p \le N} p\right) = O(N^2 / \ln N)$. With BSGS for discrete logs: $O(\sqrt{p} \log p)$ per prime, giving $O(N^{3/2}/\ln N)$ total.
• Space: $O(N)$ for the sieve and the largest discrete log table.

Answer

#include <bits/stdc++.h>
using namespace std;
const int MOD = 1e9 + 7;
long long power(long long b, long long e, long long m) {
    long long r = 1; b %= m;
    while (e > 0) { if (e & 1) r = r * b % m; b = b * b % m; e >>= 1; }
    return r;
}
int main() {
    const int LIM = 200;
    vector<bool> sieve(LIM, true);
    sieve[0] = sieve[1] = false;
    for (int i = 2; i * i < LIM; i++)
        if (sieve[i]) for (int j = i*i; j < LIM; j += i) sieve[j] = false;
    long long ans = 0;
    for (int p = 2; p < LIM; p++) {
        if (!sieve[p]) continue;
        for (int x = 1; x < p; x++)
            for (int y = x+1; y < p; y++)
                if (power(x, y, p) == power(y, x, p))
                    ans = (ans + x + y) % MOD;
    }
    cout << ans << endl;
    return 0;
}

"""
Problem 801: x^y ≡ y^x (mod p)
For a prime $p$, let $f(p) = \sum (x+y)$ where the sum is over all pairs $(x,y)$ with $1 \le x < y < p$ such that $x^y \equiv y^x \pmod{p}$. Find $\sum_{p \le 10^6} f(p) \bmod 1\,000\,000\,007$.
"""

def solve():
    MOD = 10**9 + 7
    def f_small(p):
        """Brute force for small primes."""
        total = 0
        for x in range(1, p):
            for y in range(x+1, p):
                if pow(x, y, p) == pow(y, x, p):
                    total += x + y
        return total
    # Demo for small primes
    result = 0
    sieve = [True] * 200
    sieve[0] = sieve[1] = False
    for i in range(2, 15):
        if sieve[i]:
            for j in range(i*i, 200, i):
                sieve[j] = False
    for p in range(2, 200):
        if sieve[p]:
            result = (result + f_small(p)) % MOD
    return result

answer = solve()
print(answer)