Best way to simulate star sizes to scale in celestial sphere
Your maths looks ok, bar the fact that $1/$parallax is a biased estimate of the distance (but that can be forgiven so long as you are using data where the parallax uncertainties are much smaller than the parallax).

Your main problem is that stars do indeed have a vast range of sizes. Thus if you really do want to show the relative sizes of the stars you have a dynamic range problem.

The conventional way to deal with this would be to use a logarithmic scale, such that each increment in size corresponds to a multiple of the one below. e.g. Make your plotted radii proportional to $\log_{10} R_{\rm rel}$.

EDIT: In response to altered emphasis of question.

It seems you no longer want to scale the stars according to their actual radius, but rather according to their brightness. The logarithmic scaling is therefore already handled by the apparent magnitude (which is on a logarithmic scale).

Therefore all you need to do is decide what your smallest and largst radius are and compare that with the brightest and faintest magnitudes you want to plot.

e.g. If your biggest star is 5 units and your smallest is 0 units, and your bright and faint limits are $m_{\rm bright}$ and $m_{\rm faint}$ respectively, then the size of an arbitrary star of magnitude $m$ is
$$ r =\frac{5(m -m_{\rm faint})}{m_{\rm bright} -  m_{\rm faint}} \ .$$

And if you want your faintest stars to be size 1 and brightest to be 5, this modifies to
$$ r =\frac{4(m -m_{\rm faint})}{m_{\rm bright} -  m_{\rm faint}} +1 \ .$$

Yet another alternative would be to make the area proportional to magnitude. Again, scaling between radii of 5 and 1:

$$ r^2 =\frac{24(m -m_{\rm faint})}{m_{\rm bright} -  m_{\rm faint}} +1 \ .$$