Definitions
Dispersion Relation
In general, to find the phase velocity or other properties of a wave mode, one starts by deriving a dispersion relation, $D(\mathbf{k}, \omega)$. The dispersion relation is derived from the equation:
$$
\mathbf{n} \times \left( \mathbf{n} \times \mathbf{E} \right) + \overleftrightarrow{\mathbf{K}} \cdot \mathbf{E} = 0 \tag{1}
$$
where $\mathbf{n}$ is the index of refraction and $\overleftrightarrow{\mathbf{K}}$ is the dielectric tensor (e.g., see https://physics.stackexchange.com/a/138460/59023 for further discussion of derivation of Equation 1). We rewrite Equation 1 in the tensor form $\overleftrightarrow{\mathbf{D}} \cdot \mathbf{E} = 0$. If $\overleftrightarrow{\mathbf{D}}$ has a determinant that goes to zero, then there is a non-trivial solution for $\mathbf{E}$.
The solutions to the dispersion relation give the wave frequency, $\omega$, as a function of the wavevector, $\mathbf{k}$.
MHD Fast Mode
The simplest form of the fast mode dispersion relation for arbitrary propagation angle, $\theta$, with respect to $\mathbf{B}_{o}$ starts by linearizing Faraday's law, the mass continuity equation, the momentum equation, and an equation of state (e.g., our expression for the speed of sound in the Definitions section above). By linearize, I mean that we assume any relevant quantity can be rewritten as $Q = Q_{o} + \delta Q$ where the $Q_{o}$ is the quasi-static term (i.e., roughly constant or background term) and $\delta Q$ is the fluctuating term. We also assume that $\delta Q \propto e^{i \left( \mathbf{k} \cdot \mathbf{x} - \omega \ t \right)}$ (actually this comes from a Fourier transform but I it does not change much in way of the process).
The latter assumption leads to the following $\partial / \partial t \rightarrow -i \ \omega$ and $\nabla \rightarrow i \ \mathbf{k}$ and that either of these operators acting on any $Q_{o}$ term goes to zero.
To start, the four equations I mentioned are given, to first order, as:
$$
\begin{align}
\frac{ \partial \ \delta \rho }{ \partial t } + \nabla \cdot \left( \mathbf{U}_{o} \ \delta \rho + \delta \mathbf{U} \ \rho_{o} \right) & = 0 \tag{2a} \\
\rho_{o} \frac{ \partial \ \delta \mathbf{U} }{ \partial t } & = \delta \mathbf{j} \times \mathbf{B}_{o} + \mathbf{j}_{o} \times \delta \mathbf{B} - \nabla \left( \delta P \right) \tag{2b} \\
\frac{ \partial \ \delta \mathbf{B} }{ \partial t } & = \nabla \times \left( \delta \mathbf{U} \times \mathbf{B}_{o} + \mathbf{U}_{o} \times \delta \mathbf{B} \right) \tag{2c} \\
\delta P & = \gamma \left( \frac{ P_{o} }{ \rho_{o} } \right) \delta \rho \tag{2d}
\end{align}
$$
We can greatly simplify this in a few ways. The first is that we just move into the bulk flow rest frame, i.e., $\mathbf{U}_{o} \rightarrow 0$. The second is to recognize that there are no $\mathbf{j}_{o}$ terms, since those would imply that $\nabla \times \mathbf{B}_{o} \neq 0$ which would violate our assumptions. Then we have:
$$
\begin{align}
-i \ \omega \ \delta \rho + i \ \rho_{o} \ \mathbf{k} \cdot \delta \mathbf{U} & = 0 \tag{3a} \\
-i \ \omega \ \rho_{o} \ \delta \mathbf{U} & = \frac{ i }{ \mu_{o} } \left( \mathbf{k} \times \delta \mathbf{B} \right) \times \mathbf{B}_{o} - i \ \mathbf{k} \ \delta P \tag{3b} \\
-i \ \omega \ \delta \mathbf{B} & = i \ \mathbf{k} \times \left( \delta \mathbf{U} \times \mathbf{B}_{o} \right) \tag{3c} \\
\delta P & = C_{s}^{2} \ \delta \rho \tag{3d}
\end{align}
$$
We can solve for $\delta \rho$ in Equation 3a to find a new relationship for $\delta P$ from Equation 3d to then insert into Equation 3b. This new form of Equation 3b combined with 3c gives (after some symbol gymnastics):
$$
\omega^{2} \ \delta \mathbf{U} = \frac{ 1 }{ \mu_{o} \ \rho_{o} } \left\{ \mathbf{k} \times \left[ \mathbf{k} \times \left( \delta \mathbf{U} \times \mathbf{B}_{o} \right) \right] \right\} \times \mathbf{B}_{o} + C_{s}^{2} \ \mathbf{k} \left( \mathbf{k} \cdot \delta \mathbf{U} \right) \tag{4}
$$
To further simplify, we define the coordinate basis such that the following are satisfied:
$$
\begin{align}
\mathbf{B}_{o} & = B_{o} \ \hat{\mathbf{z}} \tag{5a} \\
\mathbf{k} & = \left( k \ \sin{\theta}, 0, k \ \cos{\theta} \right) \tag{5b}
\end{align}
$$
and then define $V_{ph} = \omega/k$ as the phase speed and rewrite Equation 4 in matrix form and after (significant) simplification we have:
$$
\begin{align}
0 & = \left[
\begin{array}{ c c c }
V_{ph}^{2} - C_{s}^{2} \ \sin^{2}{\theta} - V_{A}^{2} & 0 & - C_{s}^{2} \ \sin{\theta} \ \cos{\theta} \\
0 & V_{ph}^{2} - V_{A}^{2} \ \cos^{2}{\theta} & 0 \\
- C_{s}^{2} \ \sin{\theta} \ \cos{\theta} & 0 & V_{ph}^{2} - C_{s}^{2} \ \cos^{2}{\theta}
\end{array} \right] \cdot \left[
\begin{array}{ c }
\delta U_{x} \\
\delta U_{y} \\
\delta U_{z} \\
\end{array} \right]
\tag{6}
\end{align}
$$
The determinant of the first matrix gives the dispersion relation:
$$
D\left( \mathbf{k}, \omega \right) = \left( V_{ph}^{2} - V_{A}^{2} \ \cos^{2}{\theta} \right) \left[ V_{ph}^{4} - V_{ph}^{2} \left( V_{A}^{2} + C_{s}^{2} \right) + V_{A}^{2} \ C_{s}^{2} \ \cos^{2}{\theta} \right] = 0 \tag{7}
$$
It is then easy to show that there are three roots, one of which is the fast mode root given by:
$$
\begin{align}
2 \ V_{f}^{2} & = \left( V_{A}^{2} + C_{s}^{2} \right) + \sqrt{ \left( V_{A}^{2} + C_{s}^{2} \right)^{2} - 4 \ V_{A}^{2} \ C_{s}^{2} \ \cos^{2}{\theta} } \tag{8a} \\
& = \left( V_{A}^{2} + C_{s}^{2} \right) + \sqrt{ \left( V_{A}^{2} - C_{s}^{2} \right)^{2} + 4 \ V_{A}^{2} \ C_{s}^{2} \ \sin^{2}{\theta} } \tag{8b}
\end{align}
$$
where I have used $V_{f} = V_{ph}$ for the fast mode.
What am I doing wrong?
Nothing, just not catching that these are not truly equivalent expressions. For instance, the form shown after your link to the Wikipedia article is missing a square in the numerator (i.e., the parentheses should be squared), which just reduces it to the second expression you show.
Let's start with the form:
$$
V_{f} = c \ \sqrt{ \frac{ \sigma \left( 1 + \beta_{s}^{2} \right) + \beta_{s}^{2} }{ \left( 1 + \beta_{s}^{2} \right) \left( 1 + \sigma \right) } } \tag{9}
$$
Now let us take the limit where $\beta_{s}^{2} \ll 1$ and $\sigma \ll 1$, then Equation 9 goes to:
$$
\begin{align}
\lim_{\beta_{s}^{2} \ll 1 \\ \sigma \ll 1} \left[ \left( 1 + \beta_{s}^{2} \right) \left( 1 + \sigma \right) \right]^{-1/2} & \approx 1 - \frac{ \beta_{s}^{2} }{ 2 } - \frac{ \sigma }{ 2 } \tag{10a} \\
V_{f} & \simeq c \ \sqrt{ \left[ \sigma \left( 1 + \beta_{s}^{2} \right) + \beta_{s}^{2} \right] \left( 1 - \frac{ \beta_{s}^{2} }{ 2 } - \frac{ \sigma }{ 2 } \right) } \tag{10b} \\
& = c \ \sqrt{ \sigma + \beta_{s}^{2} - \frac{ 1 }{ 2 } \left( \sigma^{2} + \beta_{s}^{4} + \beta_{s}^{2} \ \sigma^{2} + \beta_{s}^{4} \ \sigma \right) } \tag{10c} \\
& \approx c \ \sqrt{ \sigma + \beta_{s}^{2} } \tag{10d} \\
V_{f} & \approx \sqrt{ V_{A}^{2} + C_{s}^{2} } \tag{10e}
\end{align}
$$
where Equation 10e is the form shown in your second and fourth expressions (and the result found in the limit as $\theta \rightarrow \pi/2$ for Equations 8a and 8b) and Equation 10d is your third expression.
Justify Assumptions
I wrote a detailed answer at https://physics.stackexchange.com/a/179057/59023 that provides several estimates for typical values of $V_{A}$ and $C_{s}$ in the interplanetary medium.
Under most conditions both $C_{s}$ and $V_{A}$ satisfy < 100 km/s. The typical values one often quotes for reference are in the ~20-50 km/s range. The only place where $V_{A}$ approaches the speed of light, of which I am aware, is in the auroral acceleration region near Earth. Similarly, $C_{s}$ is almost never near $c$, except perhaps within stars.
Side Note: As an aside, I have been working in space plasma physics for over a decade and have never seen the form for the fast mode speed given by Equation 9. This is not to say that it is not used or irrelevant, just interesting because I have read a substantial fraction of the literature on these modes.
