Find an orthonormal basis and the signature of the quadratic form 
Consider the quadratic form given by the matrix below (in the canonical basis)
\begin{pmatrix}
1 & 1 & -1\\
1 & 1 & 3\\
-1 & 3 & 1
\end{pmatrix}
Find an orthonormal basis of it and find its signature.


First I calculated the eigenvalues, which are $4, \frac{-1+ \sqrt{17}}{2}, \frac{-1-\sqrt{17}}{2}$. Then I calculated the eigenvectors associated to $4$ and $\frac{-1+ \sqrt{17}}{2}$ and normalized them, which gave me
\begin{align}
    e_1 &= \frac{1}{\sqrt{2}}\begin{bmatrix}
           0 \\
           1 \\
            1
         \end{bmatrix}, \quad
e_2 = \frac{\sqrt{2}}{\sqrt{17+3\sqrt{17}}}\begin{bmatrix}
           -\frac{3+\sqrt{17}}{2} \\
           -1 \\
            1
         \end{bmatrix}
  \end{align}
And the third vector of the basis I want to be orthogonal to $e_1$ and $e_2$, so
$$e_3 = \frac{1}{\sqrt{17+3\sqrt{17}}} e_1 \wedge e_2 =\frac{1}{\sqrt{17+3\sqrt{17}}} \begin{bmatrix}
           2 \\
           -\frac{3+\sqrt{17}}{2} \\
            \frac{3+\sqrt{17}}{2}
         \end{bmatrix} $$
I can't detail the calculations because they are very big. Perhaps someone can confirm the results. For the signature I know that the two possibilities are $(0,3)$ and $(2,1)$ but I don't know how to find the right one.

This business with eigenvalues and eigenvectors is not how you diagonalize a quadratic form. It will give a correct result if done correctly, but is way too long and computationally painful. Simply use a "complete the square" method. In the canonical coordinates, your quadratic form is
$$x^2 + y^2+z^2+2xy-2xz+6yz.$$
Now let us complete the squares:
$$\begin{align}
  & x^2 + y^2+z^2+2xy-2xz+6yz \\
&= (x^2+2xy +y^2) + z^2-2xz+6yz \\
 &=  (x+y)^2 + z^2 -2xz+6yz \\
&= (x+y)^2 + z^2 -2z(x-3y) \\
&= (x+y)^2 + (z-(x-3y))^2 - (x-3y)^2.
\end{align}$$
This means that the quadratic form is simply $(x')^2+(y')^2-(z')^2$ with the change of coordinates $$x'=x+y,\quad y'=-x+3y+z,\quad z'=x-3y.$$
This already clearly shows that the signature is $(2,1)$.
So to find your orthogonal basis, you just have to invert the matrix
$$\begin{pmatrix}
1 & 1 & 0 \\
-1 & 3 & 1 \\
1 & -3 & 0
\end{pmatrix}$$
and the basis will be given by the columns of the inverse. No need to compute eigenvalues or eigenvectors with complicated expressions. This gives
$$\frac{1}{4}\begin{pmatrix}
3 & 0 & 1 \\
1 & 0 & -1 \\
0 & 4 & 4
\end{pmatrix},$$
so an orthonormal basis is $e_1 = \begin{bmatrix} 3/4 \\ 1/4 \\ 0\end{bmatrix}$, $e_2 = \begin{bmatrix} 0 \\ 0 \\ 1\end{bmatrix}$ and $e_3 = \begin{bmatrix} 1/4 \\ -1/4 \\ 1\end{bmatrix}$.
