Computing solutions to the congruence $1^n+2^n+\cdots +n^n\equiv p(\mathrm{m}\mathrm{o}\mathrm{d}n)$

It is well-known that the congruence $\sum_{i=1}^{ n} i^{ n} \equiv 1 \pmod{n}$hasexactlyfivesolutions:${1,2,6,42,1806}$. In this work, we characterize the solutions to the congruence $1^n + 2^n + \dotsb + n^n\equiv p \pmod{n}$foreveryprime$p$. This characterization leads to an algorithm for computing all such solutions, when there is a finite number of them. More generally, our algorithm enables computing all the solutions below a much higher bound as compared to what can be achieved by a naive exhaustive search.