OFFSET
1,1
LINKS
Chai Wah Wu, Table of n, a(n) for n = 1..558
Project Euler, Problem 801. x^y = y^x.
FORMULA
a(n) = Sum_{k=0..p-1} b(k)^2, where b(k) is the number of solutions to x^y == k (mod p) with p = prime(n), 0 < x <= p and 0 < y <= p-1. - Jinyuan Wang, Jun 19 2022
EXAMPLE
For p=2:
p x y x^y mod p y^x mod p
- - - --------- ---------
2 1 1 1 1
2 2 2 0 0
Solutions: 2
.
For p=3:
p x y x^y mod p y^x mod p
- - - --------- ---------
3 1 1 1 1
3 1 4 1 1
3 2 2 1 1
3 2 4 1 1
3 3 3 0 0
3 3 6 0 0
3 4 1 1 1
3 4 2 1 1
3 4 4 1 1
3 4 5 1 1
3 5 4 1 1
3 5 5 2 2
3 6 3 0 0
3 6 6 0 0
Solutions: 14
PROG
(Python)
from sympy import prime
def f(n):
S = 0
n2n=(n*n) - n
for x in range(1, n2n + 1):
for y in range(x + 1 , n2n + 1):
if ((pow(x, y, n) == pow(y, x, n))):
S += 2
return S + n2n
def a(n): return f(prime(n))
(Python)
from sympy import prime
from sympy.ntheory.residue_ntheory import nthroot_mod
def A355069(n):
p = prime(n)
return sum(sum(len(nthroot_mod(k, y, p, True)) for y in range(1, p))**2 for k in range(p)) # Chai Wah Wu, Aug 31 2022
(PARI) a(n) = my(p=prime(n), v=vector(p)); for(x=1, p, for(y=1, p-1, v[1+lift(Mod(x, p)^y)]++)); sum(i=1, p, v[i]^2); \\ Jinyuan Wang, Jun 19 2022
CROSSREFS
KEYWORD
nonn
AUTHOR
Darío Clavijo, Jun 17 2022
EXTENSIONS
a(10)-a(26) from Michel Marcus, Jun 18 2022
More terms from Jinyuan Wang, Jun 19 2022
STATUS
approved