\(x=0\) のとき、条件を満たすのは \(y=0\) のみです。 また、\(y=0\) となるのは \(x=0\) のときに限ります。
以下では \(1\leq x, y \leq P-1\) とします。
ここで、\(r\) を法 \(P\) に対する原始根とします。

\(1 \leq a,b \leq P-1\) に対して \(x \equiv r^a, y \equiv r^b \pmod{P}\) とすると、\(x\) と \(a\) 、\(y\) と \(b\) はそれぞれ一対一に対応します。
以下の式変形が成り立ちます。

\(x^n \equiv y \pmod{P}\)
\(\Leftrightarrow\) \(r^{an} \equiv r^{b} \pmod{P}\)
\(\Leftrightarrow\) \(an \equiv b \pmod {P-1}\)

よって、 \(an \equiv b \pmod {P-1}\) を満たす正整数 \(n\) が存在するような \((a, b)\) の組の数を数えることに帰着されました。
これは

\(\sum_{a=1}^{P-1} \frac{P-1}{gcd(P-1, a)}\)

を求めれば良いですが、この計算量は \(O(P)\) です。
高速化のためにシグマ計算を \(g=gcd(P-1, a)\) ごとにまとめて考えると、

\(\sum_{g} \frac{P-1}{g} \times f(g)\)

となります。ただし、\(f(g)\) は \(1 \leq a \leq P-1\) かつ \(g=gcd(P-1,a)\) を満たす \(a\) の数です。
\(g\) としてあり得るものの数は高々 \(P-1\) の約数の数です。\(P-1\) の約数の数を \(d\) とすると、\(f(g)\) は \(g\) の降順に計算することで \(O(d^2)\) で全て求めることができます。
よってこの問題は \(O(\sqrt{P} + d^2)\) で解くことができました。
\(10^{12}\) 未満の高度合成数で最大のものは \(963761198400 =2^6*3^4*5^2*7*11*13*17*19*23\) で、この約数は \(6720\) 個なので、十分高速です。