|
%I #45 Sep 01 2022 19:27:35
%S 2,14,104,366,1550,3048,6272,9774,14894,34664,48750,84456,108320,
%T 128814,128846,209768,255374,424680,479886,563150,700704,782574,
%U 712334,1068320,1614336,1649000,1721454,1527566,2299752,2328704,3654126,3428750,3834656,4201134,4596584
%N a(n) is the number of solutions to x^y == y^x (mod p) where 0 < x,y <= p^2 - p and p is the n-th prime.
%H Chai Wah Wu, <a href="/A355069/b355069.txt">Table of n, a(n) for n = 1..558</a>
%H Project Euler, <a href="https://projecteuler.net/problem=801">Problem 801. x^y = y^x</a>.
%F 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
%e For p=2:
%e p x y x^y mod p y^x mod p
%e - - - --------- ---------
%e 2 1 1 1 1
%e 2 2 2 0 0
%e Solutions: 2
%e .
%e For p=3:
%e p x y x^y mod p y^x mod p
%e - - - --------- ---------
%e 3 1 1 1 1
%e 3 1 4 1 1
%e 3 2 2 1 1
%e 3 2 4 1 1
%e 3 3 3 0 0
%e 3 3 6 0 0
%e 3 4 1 1 1
%e 3 4 2 1 1
%e 3 4 4 1 1
%e 3 4 5 1 1
%e 3 5 4 1 1
%e 3 5 5 2 2
%e 3 6 3 0 0
%e 3 6 6 0 0
%e Solutions: 14
%o (Python)
%o from sympy import prime
%o def f(n):
%o S = 0
%o n2n=(n*n) - n
%o for x in range(1, n2n + 1):
%o for y in range(x + 1 , n2n + 1):
%o if ((pow(x, y, n) == pow(y, x, n))):
%o S += 2
%o return S + n2n
%o def a(n): return f(prime(n))
%o (Python)
%o from sympy import prime
%o from sympy.ntheory.residue_ntheory import nthroot_mod
%o def A355069(n):
%o p = prime(n)
%o 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
%o (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
%Y Cf. A000040, A355419.
%K nonn
%O 1,1
%A _Darío Clavijo_, Jun 17 2022
%E a(10)-a(26) from _Michel Marcus_, Jun 18 2022
%E More terms from _Jinyuan Wang_, Jun 19 2022
|
