Another relatively straightforward problem. This was because much simpler when doing the problems in logarithms rather than powers because computers are much faster at calculating large logarithms than they are in calculating large powers.
We are looking for all values where
\[800800^{800800}>p^qq^p\]which is equivalent to
\[800800\log {800800}>q\log p+p\log q\]I estimated the value for the largest prime required $p$ assuming that $q=2$.
\[\begin{alignat}{2} 800800^{800800}&=2^qq^2 \\ \log({800800^{800800}})&=q\log2+2\log q \\ \end{alignat}\]Since the value $q\log2$ will generally be larger than $2\log q$, I reduced the problem to the estimate
\[\begin{alignat}{2} \log({800800^{800800}})&=q\log2 \\ \frac{\log({800800^{800800}})}{\log2}&=q \\ \end{alignat}\]I used that value, which was roughly $15$ million, and used a sieve to get all required primes.
From here it was just a matter of doing the searching. I could have made it more efficient as this code searches for all approximately $1$ billion solutions. Code runs in 200 seconds using Python 3.10.8.
import math
def sieve(top):
sieve = [False, False] + [True] * (top - 1)
primes = []
for p, isprime in enumerate(sieve):
if not isprime:
continue
primes.append(p)
num = p * p
while num <= top:
sieve[num] = False
num += p
return primes
def solve(n):
logn = math.log(n)
top = int(logn / math.log(2)) + 1
primes = sieve(top)
lenprimes = len(primes)
count = 0
for i in range(lenprimes - 1):
p = primes[i]
logp = math.log(p)
for j in range(i + 1, lenprimes):
q = primes[j]
if q*logp + p*math.log(q) <= logn:
count += 1
elif j == i + 1:
return count
else:
break
if __name__ == '__main__':
import sys
n = eval(sys.argv[1])
print(solve(n))
import pytest
from problem import solve
_test_solve = (
pytest.param(800**1, 2, id='800**1'),
pytest.param(800**2, 5, id='800**2'),
pytest.param(800**3, 10, id='800**3'),
pytest.param(800**4, 13, id='800**3'),
pytest.param(800**800, 10790, id='800**800'),
)
@pytest.mark.parametrize('n,expect', _test_solve)
def test_solve(n, expect):
assert expect == solve(n)