I would strongly advice against using closed form solutions since they tend to be numerically very unstable. You need to take extreme care in the way and order of your evaluations of the discriminant and other parameters.
The classical example is the one for the quadratic equation $ax^2+bx+c=0$. Calculating the roots as $$x_{1,2} = \frac{-b \pm \sqrt{b^2-4ac}}{2a}$$ will get you into trouble for polynomials where $b\gg 4ac$ since then you get cancellation in the numerator. You need to calculate $$x_1 = \frac{-(b+sign(b)\sqrt{b^2-4ac})}{2a}; x_2 = \frac{c}{a}\frac{1}{x_1}$$.
Higham in his masterpiece "Accuracy and Stability of Numerical Algorithms" (2nd ed, SIAM) uses a direct search method to find coefficients of a cubic polynomial for which the classical analytical cubic solution gives very inaccurate results. The example he gives is $[a,b,c] = [1.732, 1, 1.2704]$. For this polynomial the roots are well-separated and hence the problem is not ill-conditioned. However, if he calculates the roots using the analytical approach, and evaluates the polynomial in these roots, he obtains a residue of $\mathcal{O}(10^{-2})$ while using a stable standard method (the companion matrix method), the residue is of order $\mathcal{O}(10^{-15})$. He proposes a slight modification to the algorithm, but even then, he can find a set of coefficients leading to residues of $\mathcal{O}(10^{-11})$ which is definitely not good. See p480-481 of the above mentioned book.
In your case, I would apply Bairstow's method. It uses an iterative combination of Newton iteration on quadratic forms (and then the roots of the quadratic are solved) and deflation. It is easily implemented and there are even some implementations available on the web.
