We study the distribution of the position of the rightmost particle $x_{\max}$ in a $N$-particle Riesz gas in one dimension confined in a harmonic trap. The particles interact via long-range repulsive potential, of the form $r^{-k}$ with $-2<k<\infty$ where $r$ is the inter-particle distance. In equilibrium at temperature $O(1)$, the gas settles on a finite length scale $L_N$ that depends on $N$ and $k$. We numerically observe that the typical fluctuation of $y_{\max} = x_{\max}/L_N$ around its mean is of $O(N^{-η_k})$. Over this length scale, the distribution of the typical fluctuations has a $N$ independent scaling form. We show that the exponent $η_k$ obtained from the Hessian theory predicts the scale of typical fluctuations remarkably well. The distribution of atypical fluctuations to the left and right of the mean $\langle y_{\max} \rangle$ are governed by the left and right large deviation functions, respectively. We compute these large deviation functions explicitly $\forall k>-2$. We also find that these large deviation functions describe a pulled to pushed type phase transition as observed in Dyson's log-gas ($k\to 0$) and $1d$ one component plasma ($k=-1$). Remarkably, we find that the phase transition remains $3^{\rm rd}$ order for the entire regime. Our results demonstrate the striking universality of the $3^{\rm rd}$ order transition even in models that fall outside the paradigm of Coulomb systems and the random matrix theory. We numerically verify our analytical expressions of the large deviation functions via Monte Carlo simulation using an importance sampling algorithm.