For how many residues $x \pmod{p}$ is $x^{-1} \in [1, (p-1)/2]$?

Let me try to explain why this problem is hard (concurring with Alexey Ustinov here). Fourier analysis shows that $$ N_p = \frac{1}{p^2} \sum_{m=0}^{p-1} \sum_{n=0}^{p-1} \left( \sum_{x=1}^{p-1} e_p( nx + mx^{-1} )\right) \left(\sum_{x=1}^{(p-1)/2} e_p(-nx) \right) \left(\sum_{x=1}^{(p-1)/2} e_p(-mx) \right) $$

where $ \left( \sum_{x=1}^{p-1} e_p( nx + mx^{-1} \right)$ is a Kloosterman sum $K(mn;c)$ and $ \left(\sum_{x=1}^{(p-1)/2} e_p(-nx) \right) \left(\sum_{x=1}^{(p-1)/2} e_p(-mx) \right)$ is large for $n$, $m$ close to zero and small elsewhere.

The term $n=0, m=0$ gives the main term, where $n=0$ and $m \neq 0$ or vice versa gives a small and easily computable term, and the interesting error term comes from the range where $n$ and $m$ are nonzero. This can be represented as a linear combination of Kloosterman sums $K(mn;p)$, with the largest contributions coming from $mn = \pm 1$.

We expect the sums $K(a;p)/\sqrt{p} $ for any $a$ to be equidistributed as random variables with a Sato-Tate (semicircular) distribution. It's reasonable to also conjecture independence over different $a$, which would enable calculating the distribution of any finite sum of Kloosterman sums. (In fact this follows from a conjecture I made in my PhD thesis). This is an infinite sum of Kloosterman sums, but it is probably reasonable to make a similar independence statement, which should give a Gaussian distribution or one very similar to it.

However, almost nothing is known about the distribution of $K(a;p)$ for fixed $a$ and varying $p$. (The case of fixed $p$ and varying $a$ is well-known by work of Deligne and Katz). Even the first moment is not known, although we do have a bound for the average of $K(a;c)$ over all numbers $c$, proved using the Kusnetsov formula.

This sum of Kloosterman sums looks much harder to analyze than the distribution of a single Kloosterman sum, so it is probably also too hard for current techniques.

The (not very surprising) fact that the mean is $1/2$ is known, the derivation uses Kloosterman sums (again, not too surprising), and a good reference for this fact and its applications is Mike Rubinstein's paper in Integers:

Rubinstein, Michael O., The distribution of solutions to $XN=N \pmod a$ with an application to factoring integers, Integers 13, Paper A12, 20 p. (2013). ZBL1294.11140.

Probably asypmtotic normality of $S_n$ is an open problem. Numerical experiments indeed give nice bell curves not only for the square $[1, (p-1)/2]$ but for arbitrary rectangles $\subset [0,p)^2$ as well. It is natural to conjecture that for small squares $S_1$, $S_2$ (of area $\asymp p^2$) numbers of points $N(S_i)$ $(i=1,2)$ of modular hyperbola $xy\equiv 1 \pmod p$ in these squares are independent normal random variables with expectations $\mathrm{Area}(S_i)/p$ and variances $\sim c\,\mathrm{Area}(S_i)p\log^ap$ for some $a$. (Variances must be proportional to the areas because both of them are additive.) I don't know te precise value of $a$, but I think that some power of $\log p$ in the standard bound $$\left|N(S)-\frac{\mathrm{Area}(S)}{p}\right|\ll\sqrt p\log^2p$$ is unremovable because normal distribution has unbounded support.