Too easy for a P6 (in fact, I solved this instead of P4... I suck at geo) I will come to uploading my solution when I've got time.

Weird sol, but seems to work anyway: Fix any $\epsilon > 0$, and we will show that the number of solutions is $\le \epsilon p$. Let $S$ be the set of solutions to $x^{p-1} \equiv 1 \pmod {p^2}$ with $1 \le x \le p$. Let $T$ denote the set of solutions $x$ with $1 \le x \le p^2$. Assume $|S| > \epsilon p$ happens infinitely often. It is well known that (say, from the existence of a primitive root modulo $p^2$) $T$ contains exactly $p-1$ elements. Now note that if $x, y \in S$, then $xy \in T$, as $1 \le xy \le p^2$ and $(xy)^{p-1} \equiv 1 \pmod {p^2}$. For each $u \in T$, let $c_u$ be the number of solutions to $u = xy$ for $x, y \in S$. Then we have \[ \left(\sum_{s \in S} s \right)^2 = \sum_{t \in T} tc_t \le \sum_{t \in T} p^2c_t. \]Now note that the left hand side is $\ge \frac{|S|^2(|S|-1)^2}{4}$, so there exists a constant $C (= 10^{-10})$ and $t \in T$ such that $c_t \ge C \epsilon^4 p$. But this means that $\tau(t) \ge C\epsilon^4 \sqrt{t}$ holds infinitely often, which contradicts this. We’re done.

The above proof gives $\mathcal O(p^{3/4}),$ what is the best bound?

$a$ must be in the form of $g^{xp}$ mod $p^2$ while not exceeding $n$, and let such $a$s be $a_1,a_2,\ldots,a_k$. Clearly all $a_ia_j$ is also in the form of $g^{xp}$ mod $p^2$. Now write all $k^2$ values of $a_ia_j$. Then each number (call it $t$) can be written out at most $\sigma_0(t)$ times. It isn't hard to show that these are less than $100t^{0.49}$. (Better bounds exist but this is enough to solve the problem.) Thus same numbers can be written at max $100p^{0.98}$ times, where there are only $p-1$ numbers of the form $g^{xp}$ mod $p^2$. Thus $k^2<100p^{1.98}$ and $k<10p^{0.99}$. Set $K$ as $(10\times2^{2024})^{100}$ and we're done.

EthanWYX2009 wrote: The above proof gives $\mathcal O(p^{3/4}),$ what is the best bound? It might not be true, but I heard from my friend that we can prove the $o (p^{\frac{1}{2} + \varepsilon})$ bound

starchan wrote:

What is Erdos Szemedei theory？？？

lminsl wrote: Weird sol, but seems to work anyway: Fix any $\epsilon > 0$, and we will show that the number of solutions is $\le \epsilon p$. Let $S$ be the set of solutions to $x^{p-1} \equiv 1 \pmod {p^2}$ with $1 \le x \le p$. Let $T$ denote the set of solutions $x$ with $1 \le x \le p^2$. Assume $|S| > \epsilon p$ happens infinitely often. It is well known that (say, from the existence of a primitive root modulo $p^2$) $T$ contains exactly $p-1$ elements. Now note that if $x, y \in S$, then $xy \in T$, as $1 \le xy \le p^2$ and $(xy)^{p-1} \equiv 1 \pmod {p^2}$. For each $u \in T$, let $c_u$ be the number of solutions to $u = xy$ for $x, y \in S$. Then we have \[ \left(\sum_{s \in S} s \right)^2 = \sum_{t \in T} tc_t \le \sum_{t \in T} p^2c_t. \]Now note that the left hand side is $\ge \frac{|S|^2(|S|-1)^2}{4}$, so there exists a constant $C (= 10^{-10})$ and $t \in T$ such that $c_t \ge C \epsilon^4 p$. But this means that $\tau(t) \ge C\epsilon^4 \sqrt{t}$ holds infinitely often, which contradicts this. We’re done. A slightly different finish which gives a bound of $o(n^{\frac{1}{2}+\epsilon})$, for any $\epsilon > 0$: By counting the number of triples $(x,y,u)$ where $x,y \in S, u\in T$ and $xy = u$, we get that $$|S|^2 \leq \sum_{t\in T} d(s) \leq (p-1) o(p^{2\epsilon}) \implies |S| = o(p^{\frac{1}{2}+\epsilon})$$

@starchan, Archit isn't it more reasonable to use density to improve bounding conditions? I mean the bound O(p^3/4) can be improved i think! Woll, can i improve by using probability ?

@EthanWYX2009 Can't the bound be improved using probability?

So, let me try to clear things up a little bit: 1) There are exactly $p-1$ solutions to $a^{p-1} \equiv 1 \pmod{p^2}$ when we consider $a$ modulo $p^2$. These are exactly the numbers of the form $a=g^{kp}$ where $g$ is a primitive root modulo $p^2$ and they form a subgroup of $(\mathbb{Z}/p^2\mathbb{Z})^{\times}$ of order $p-1$. Moreover, there is a unique such $a$ in each invertible residue class modulo $p$. 2) The problem really is about the following: Could it happen that these $p-1$ values $a$ cluster abnormally among the "small" $a \le p$? 3) As above, the reason why this is impossible is as follows: Suppose that $k$ of the values $a$ are "small" i.e. $\le p$. Then all their $k^2$ products $aa'$ are also solutions (by the subgroup property) and since $a,a' \le p$, their residues modulo $p^2$ are equal to their product. But using the bound $\tau(n) \ll n^{\varepsilon}$ we see that each such number $\le p^2$ can occur at most $p^{2\varepsilon}$ times as such a product. Hence we get at least $\gg \frac{k^2}{p^{2\varepsilon}}$ distinct solutions $a$ in total. But we have only $p-1$, so we get $k \ll p^{1/2+\varepsilon}$. 4) If we use this with the classical $\varepsilon=\frac{1}{2}$ we get $k \ll p^{3/4}$ as noted in #5. But as noted in #10, the bound for $\tau(n)$ and hence our total bound is really true for any $\varepsilon>0$. 5) What can we expect? By basic probabilistic heuristics, it is clear that we should epect $k \ll p^{\varepsilon}$ to be true for any $\varepsilon>0$. But it is also quite clear that this will be impossibly hard to prove. Even getting some exponent less than $\frac{1}{2}$ could be quite difficult. Note that even an optimal version of the Erdös-Szemeredi Sum-Product Theorem would not help us, since we have essentially already used that the product set is as large as possible in our set-up (for this reason, using the sum-product theorem for this problem is really a complete overkill...).

I tried this problem again today with the goal $\mathcal O(p^{1/2+\varepsilon})$ (I have bad memory, usually can't remember solutions after 3 months), not so hard as I thought. The key is seperating $A$(set of solution) into $B,C$ where $B$ contains all number less than $\sqrt p$ and $C$ contains all number bigger than $\sqrt p.$ We only need to estimate $|C|.$ For $c_1,c_2\in C$ $(c_1c_2)^{p-1} \equiv 1 \pmod {p^2}.$ By $p<c_1c_2<p^2,$ using Hensel Lemma $c_1c_2\mod p\notin A.$ Note that $m=c_1c_2$ is calculated at most $\tau (m)<Cm^{\delta}$ where $\delta>0$ is real. Therefore $$\binom{|C|}2\le C(p^2)^{\delta}\cdot p$$Taking $\delta =\varepsilon/2$ we get $|C|\in \mathcal O(p^{1/2+\varepsilon}),$ done$.\Box$