Well-distribution of square of an interval $[1,p^{1-\varepsilon}]$ modulo a prime $p$ Let $p$ be a large prime. I would like to say that the multi-set $[1,p^{1-\varepsilon}]^2 = \{ab \mod p: a, b \in [1,p^{1-\varepsilon}]\}$ is close to uniformly distributed, i.e. that every nonzero residue class mod $p$ occurs with almost equal frequency. 
Quantitatively, is it true that (almost?) every nonzero residue class appears at most $O(p^{1-2\varepsilon})$ times?
 A: Here is a more explicit version of Noam Elkies's comment.
Theorem. Let $p>2$ be a prime number. Let $\mathcal{U},\mathcal{V}\subseteq\{1,2,\dots,p-1\}$ be two intervals, and let $r\in\{1,2,\dots,p-1\}$
be a nonzero residue modulo $p$. Then
$$\Biggl|\sum_{\substack{u\in\mathcal{U},\ v\in\mathcal{V} \\ uv\equiv r \pmod{p}}} 1-\frac{|\mathcal{U}||\mathcal{V}|}{p-1}\Biggr|<2p^{1/2}(\log p)^2.$$
You can find a proof from scratch here (see Theorem 9).
A: If you're OK with "almost every", we can use a $L^2$ method and multiplicative character sum bounds to do better.
Let's estimate $$ \sum_{c \in \mathbb F_p^\times} \left( \left[\sum_{a,b \in I | ab=c } 1\right]  - p^{1-2\epsilon}\right)^2  $$
by Plancherel $$= \frac{1}{p-1} \sum_{\chi \in \left( \mathbb F_p^{\times} \right)^{\vee}, \chi \neq 1} \left | \sum_{a,b \in I}  \chi(ab) \right|^2 $$
and multiplicativity$$= \frac{1}{p-1} \sum_{\chi \in \left( \mathbb F_p^{\times} \right)^{\vee}, \chi \neq 1} \left | \sum_{a\in I}  \chi(a) \right|^4 $$
Now assume we have a nontrivial character sum bound $\left| \sum_{a \in I} \chi(a) \right| \leq N$ for all nontrivial characters $\chi$. Then we continue:
$$ \leq \frac{1}{p-1} \sum_{\chi \in \left( \mathbb F_p^{\times} \right)^{\vee}, \chi \neq 1} N^2 \left | \sum_{a\in I}  \chi(a) \right|^2  \leq \frac{1}{p-1} \sum_{\chi \in \left( \mathbb F_p^{\times} \right)^{\vee} } N^2 \left | \sum_{a\in I}  \chi(a) \right|^2  = N^2 p^{1-\epsilon} $$
Then we see that
$$\left[\sum_{a,b \in I | ab=c } 1\right] \leq 2 p^{1-2\epsilon} $$
for all but
$$\frac{ N^2 p^{1-\epsilon}}{ p^{2-4 \epsilon} } = N^2 p^{3\epsilon-1}$$
This is $\ll p$ as long as $N \ll p^{1 - 3\epsilon/2}$.  Using the Pólya-Vinagradov bound where $N  \approx p^{1/2}$, this works for $\epsilon <1/3$. Using the Burgess bound we can do slightly better, but not get all the way to $\epsilon = 1/2$, which would require $N= p^{1/4}$ - full square root cancellation, outside of reach.
