Detailed mathematical approach

Problem Summary

For a prime \(p\), let \(R(p)\) be the number of ordered pairs \((x,y)\) with \(1\le x,y\le p(p-1)\) such that

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

The full task is to evaluate

$S(M,N)\equiv \sum_{\substack{M \le p \le N \\ p\text{ prime}}} R(p)\pmod{993353399}.$

A direct search over all pairs is far too large, so the solution converts the congruence into a group-theoretic count that depends only on the factorization of \(p-1\).

Mathematical Approach

Fix one prime \(p\) and write \(n=p-1\). The multiplicative group \(\mathbb{F}_p^\times\) has order \(n\), and that is the structure the solution exploits.

Step 1: Separate the Multiples of \(p\)

The interval \(1,2,\dots,pn\) contains exactly \(n\) multiples of \(p\).

If exactly one of \(x\) or \(y\) is divisible by \(p\), then one side of \(x^y \equiv y^x \pmod p\) is \(0\) and the other is nonzero, so such pairs never work.

If both \(x\) and \(y\) are divisible by \(p\), then both sides are \(0\), so every such pair works. Therefore these pairs contribute

$n^2$

solutions immediately.

Everything else comes from pairs with \(p\nmid x\) and \(p\nmid y\).

Step 2: Encode Nonzero Residues by a Primitive Root

Choose a primitive root \(g\) modulo \(p\). Then every nonzero residue class can be written uniquely as \(g^a\) for some \(a\in \mathbb{Z}/n\mathbb{Z}\).

So for \(p\nmid x\) and \(p\nmid y\), write

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

with \(a,b\in \mathbb{Z}/n\mathbb{Z}\).

Now look at all integers in \(1\le x\le pn\) having the same residue modulo \(p\). They are

$r,\ r+p,\ r+2p,\ \dots,\ r+(n-1)p.$

Because \(p\equiv 1 \pmod n\), these numbers are congruent modulo \(n\) to

$r,\ r+1,\ r+2,\ \dots,\ r+(n-1),$

so each residue class modulo \(n\) occurs exactly once. Hence every pair

$\bigl(a,u\bigr)\in (\mathbb{Z}/n\mathbb{Z})^2$

corresponds to a unique integer \(x\) with

$x\equiv g^a \pmod p,\qquad x\equiv u \pmod n,$

and similarly every \((b,v)\) determines a unique \(y\).

Step 3: Count Exponent Pairs for Fixed Logarithms

Since \(g\) has order \(n\), we have

$x^y \equiv y^x \pmod p \iff g^{ay}\equiv g^{bx}\pmod p \iff ay\equiv bx\pmod n.$

Only the classes of \(x\) and \(y\) modulo \(n\) matter in the exponents, so with \(u\equiv x\pmod n\) and \(v\equiv y\pmod n\), the condition becomes

$av\equiv bu\pmod n.$

For fixed \(a\) and \(b\), define

$\Phi_{a,b}(u,v)=bu-av \pmod n.$

This is a homomorphism from \((\mathbb{Z}/n\mathbb{Z})^2\) to \(\mathbb{Z}/n\mathbb{Z}\). If

$d=\gcd(a,b,n),$

then the image of \(\Phi_{a,b}\) is exactly the subgroup of multiples of \(d\), which has size \(n/d\). Therefore the kernel has size

$\frac{n^2}{n/d}=nd.$

So for each fixed pair \((a,b)\), the number of admissible pairs \((u,v)\) is

$n\,\gcd(a,b,n).$

Summing over all \(a,b\in \mathbb{Z}/n\mathbb{Z}\) gives the nonzero contribution

$n\,H(n),\qquad H(n)=\sum_{a=0}^{n-1}\sum_{b=0}^{n-1}\gcd(a,b,n).$

Step 4: Turn the GCD Sum into a Divisor Sum

Use the standard identity

$\gcd(a,b,n)=\sum_{d\mid \gcd(a,b,n)} \varphi(d),$

where \(\varphi\) is Euler's totient function.

Substituting this into \(H(n)\) and swapping the order of summation yields

$H(n)=\sum_{d\mid n}\varphi(d)\left(\frac{n}{d}\right)^2.$

Indeed, for a fixed divisor \(d\mid n\), there are exactly \(n/d\) residue classes modulo \(n\) divisible by \(d\), both for \(a\) and for \(b\).

This formula already shows that \(H(n)\) is multiplicative in \(n\).

Step 5: Evaluate the Prime-Power Factor

Let

$n=\prod_{i=1}^t q_i^{e_i}.$

Because \(H\) is multiplicative, it is enough to evaluate one prime power:

$H(q^e)=\sum_{k=0}^{e}\varphi(q^k)\,q^{2e-2k}.$

Now \(\varphi(q^0)=1\) and \(\varphi(q^k)=q^k-q^{k-1}\) for \(k\ge 1\), so

$H(q^e)=q^{2e}+\sum_{k=1}^{e}(q^k-q^{k-1})q^{2e-2k}.$

This simplifies to

$H(q^e)=q^{2e}+q^{2e-1}-q^{e-1}=q^{e-1}(q^{e+1}+q^e-1).$

Therefore

$H(n)=\prod_{i=1}^{t} q_i^{e_i-1}(q_i^{e_i+1}+q_i^{e_i}-1),$

and the final prime-local count is

$\boxed{R(p)=n^2+n\prod_{q^e\parallel n} q^{e-1}(q^{e+1}+q^e-1),\qquad n=p-1.}$

Worked Example: \(p=5\)

Here \(n=p-1=4\), so the search range is \(1\le x,y\le 20\).

There are \(4\) multiples of \(5\) in that range, namely \(5,10,15,20\). Pairs with both entries divisible by \(5\) contribute

$4^2=16$

solutions.

Next compute

$H(4)=\sum_{d\mid 4}\varphi(d)\left(\frac{4}{d}\right)^2=\varphi(1)\cdot 4^2+\varphi(2)\cdot 2^2+\varphi(4)\cdot 1^2=16+4+2=22.$

So the nonzero part contributes

$4\cdot 22=88,$

and hence

$R(5)=16+88=104.$

A direct brute-force check over \(1\le x,y\le 20\) gives the same value.

How the Code Works

The C++, Python, and Java implementations scan the interval \([M,N]\), handle the even prime separately, and test odd candidates for primality. Only primes contribute to the final sum.

For each such prime \(p\), the implementation factors \(n=p-1\). Fast modular exponentiation is used throughout, and Pollard-Rho is used to split composite factors until the full prime factorization is known. Equal prime factors are then grouped to obtain the exponents \(e\).

Once the factorization is available, the implementation evaluates the local prime-power factors

$q^{e-1}(q^{e+1}+q^e-1)\pmod{993353399}$

and multiplies them to obtain \(H(n)\) modulo \(993353399\). Finally it adds

$n^2+nH(n)\pmod{993353399}$

to the running interval sum. The program never searches over \((x,y)\); all counting is done through the closed formula above.

Complexity Analysis

If \(W=N-M+1\), the interval scan examines \(O(W)\) candidates, skipping even numbers after \(2\). For each prime \(p\), the dominant cost is factoring \(p-1\); once the factorization is known, evaluating the product formula only needs a number of modular multiplications proportional to the number of distinct prime factors. In practice the method is fast because it replaces a quadratic search over pairs by arithmetic on the much smaller factorization of \(p-1\). Memory usage per prime stays small, essentially the temporary factor list and recursion stack.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=801
2. Primitive roots: Wikipedia - Primitive root modulo n
3. Cyclic groups: Wikipedia - Cyclic group
4. Chinese remainder theorem: Wikipedia - Chinese remainder theorem
5. Euler's totient function: Wikipedia - Euler's totient function
6. Miller-Rabin primality test: Wikipedia - Miller-Rabin primality test
7. Pollard's rho algorithm: Wikipedia - Pollard's rho algorithm