I'm polishing this post on Stack Overflow with proper math formatting and a nicer example.
The following description follows the German Wikipedia article Hauptachsentransformation. Its English counterpart, according to inter-wiki links, is principal component analysis. I find the former article a lot more geometric than the latter. The latter has a strong focus on statistical data, which might be useful for some applications.
Rotation
Your ellipse is described as
$$(x, y)\cdot\begin{pmatrix}A&C/2\\C/2&B\end{pmatrix}\cdot\begin{pmatrix}x\\y\end{pmatrix}+(D,E)\cdot\begin{pmatrix}x\\y\end{pmatrix}+F=0\tag1\\
(x, y)\cdot\left(\begin{array}{rr}4.36&-2.52\\-2.52&2.89\end{array}\right)\cdot\begin{pmatrix}x\\y\end{pmatrix}+(30.8,-0.6)\cdot\begin{pmatrix}x\\y\end{pmatrix}+81=0$$
First you identify the rotation.
Using eigenvectors
You do this by identifying the eigenvalues and eigenvectors of this $2\times2$ matrix. These eigenvectors will, after normalization, form an orthogonal matrix describing the rotation. I'll not go into the details of how to find eigenvectors of a matrix.
Using the example, you find that the eigenvalues of the matrix are $\lambda_1=1$ and $\lambda_2=6.25$. (Taking the eigenvalue with smaller absolute value first will give you the semimajor axis first in the final result.) For the former, an eigenvector is $(3,4)^T$, for the latter $(-4,3)^T$ is one. Normalizing them to unit length yields $(0.6,0.8)^T$ resp. $(-0.8,0.6)^T$. So the diagonal decomposition is
$$\left(\begin{array}{rr}4.36&-2.52\\-2.52&2.89\end{array}\right)=
\left(\begin{array}{rr}0.6&-0.8\\0.8&0.6\end{array}\right)\cdot
\left(\begin{array}{ll}1&0\\0&6.25\end{array}\right)\cdot
\left(\begin{array}{rr}0.6&0.8\\-0.8&0.6\end{array}\right)\tag2$$
The first of the three factors has the eigenvectors as columns, the last one is the transpose of the first and therefore has them as rows. The diagonal matrix in the center has the eigenvalues, in matching order.
Using rotation matrices
If you want to avoid computing eigenvalues and eigenvectors, you can also use the fact that you know what a rotation matrix looks like, and start with the equation
$$\begin{pmatrix}\cos\theta&\sin\theta\\-\sin\theta&\cos\theta\end{pmatrix}\cdot\begin{pmatrix}A&C/2\\C/2&B\end{pmatrix}\cdot\begin{pmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{pmatrix}\overset!=\begin{pmatrix}\lambda_1&0\\0&\lambda_2\end{pmatrix}$$
Now take one of the entries which has to be zero, and perform the matrix multiplication to obtain the corresponding left hand side. You will get
\begin{align*}
(B-A)\sin\theta\cos\theta &+ \tfrac12C\bigl(\cos^2\theta-\sin^2\theta\bigr) = 0 \tag3
\\
2(B-A)\sin\theta\cos\theta &= C\bigl(\sin^2\theta-\cos^2\theta\bigr)
\\
4(B-A)^2\sin^2\theta\cos^2\theta &= C^2\Bigl(\bigl(\sin^2\theta-\cos^2\theta\bigr)^2 - \bigl(\sin^2\theta+\cos^2\theta\bigr)^2 + 1\Bigr)
\\
4(B-A)^2\sin^2\theta\cos^2\theta &= C^2\Bigl(1-4\sin^2\theta\cos^2\theta\Bigr)
\\
u := 4\sin^2\theta\cos^2\theta &= \frac{C^2}{(B-A)^2+C^2}
\\
4\sin^2\theta\bigl(1-\sin^2\theta\bigr) &= u
\\
\sin^4\theta - \sin^2\theta + \tfrac14u &= 0
\\
\sin^2\theta &= \frac{1\pm\sqrt{1-u}}{2}
\end{align*}
So for our example this gives us
$$ u = \frac{C^2}{(B-A)^2+C^2} = \frac{5.04^2}{(2.89 - 4.36)^2+5.04^2} = 0.9216 $$
\begin{align*}
\sin^2\theta_1 &= \frac{1-\sqrt{0.0784}}2 = 0.36 &
\lvert\sin\theta_1\rvert &= \sqrt{0.36} = 0.6 &
\lvert\cos\theta_1\rvert &= \sqrt{1-0.36} = 0.8 \\
\sin^2\theta_2 &= \frac{1+\sqrt{0.0784}}2 = 0.64 &
\lvert\sin\theta_2\rvert &= \sqrt{0.64} = 0.8 &
\lvert\cos\theta_2\rvert &= \sqrt{1-0.64} = 0.6
\end{align*}
At this moment there appear to be $8$ solutions: $\theta_1$ and $\theta_2$, each with $2\cdot2=4$ possible sign combinations. But half of these will not satisfy equation $(3)$. The other half corresponds to rotations which differ by multiples of $90°$. So choosing between $\theta_1$ and $\theta_2$ you essentially choose which axis of the ellipse gets aligned with which axis of the coordinate system, and by choosing among the remaining sign combinations you can rotate things by $180°$ which does not change the ellipse.
For the sake of consistency with what I wrote above, I'l choose $\sin\theta=0.8$ and $\cos\theta=0.6$. Then the diagonal matrix can be computed (equivalent to $(2)$ above) as
$$\left(\begin{array}{rr}0.6&0.8\\-0.8&0.6\end{array}\right)\cdot
\left(\begin{array}{rr}4.36&-2.52\\-2.52&2.89\end{array}\right)\cdot
\left(\begin{array}{rr}0.6&-0.8\\0.8&0.6\end{array}\right)=
\left(\begin{array}{ll}1&0\\0&6.25\end{array}\right)$$
Applying the rotation
If you multiply the vector $(x,y)^T$ with the last matrix from $(2)$, then you will change the coordinate system in such a way that the mixed term vanishes, i.e. the $x$ and $y$ axes are parallel to the main axes of your ellipse. I'll write coordinates in that system as $(x',y')^T$:
$$\begin{pmatrix}x'\\y'\end{pmatrix}=\left(\begin{array}{rr}0.6&0.8\\-0.8&0.6\end{array}\right)\cdot\begin{pmatrix}x\\y\end{pmatrix}\tag4$$
The rotation of the ellipse can be read from that rotation matrix. The matrix used in $(3)$ transforms a point on the rotated ellipse into a point on the axis-aligned ellipse. So the direction is opposite to what you'd use when describing the rotation of the ellipse, and you best compute the angle from the first row of that matrix:
$$\theta=\arctan\frac{\sin\theta}{\cos\theta}=\arctan\frac{0.8}{0.6}\approx53.13°$$
Translation
If you multiply the row vector $(D,E)$ in equation $(1)$ with the first of the three matrices from $(2)$, then this effect will exactly cancel the multiplication expressed in $(4)$. Therefore in that changed coordinate system, you use the central diagonal matrix for the quadratic term, and this product for the linear term.
$$(30.8,-0.6)\cdot\left(\begin{array}{rr}0.6&-0.8\\0.8&0.6\end{array}\right)=(18,-25)\\
(x', y')\cdot\left(\begin{array}{ll}1&0\\0&6.25\end{array}\right)\cdot\begin{pmatrix}x'\\y'\end{pmatrix}+(18,-25)\cdot\begin{pmatrix}x'\\y'\end{pmatrix}+81=0\\
(1\,x'^2 + 18\,x')+(6.25\,y'^2-25\,y')+81=0$$
Now you have to complete the square independently for $x'$ and $y'$, and you end up with a form from which you can read the center coordinates.
\begin{align*}
x'^2 + 18\,x' &= \phantom{1.00}(x' + 9)^2 - 81 & x_c' &= -9 \\
6.25\,y'^2 - 25\,y' &= 6.25(y' - 2)^2 - 25 & y_c' &= \phantom+2 \\
\end{align*}
Note that the center coordinates $x_c'$ and $y_c'$ describe the center in the already rotated coordinate system. To obtain the original center you'd multiply again with the first matrix from $(2)$:
$$\begin{pmatrix}x_c\\y_c\end{pmatrix}=
\left(\begin{array}{rr}0.6&-0.8\\0.8&0.6\end{array}\right)\cdot
\left(\begin{array}{r}-9\\2\end{array}\right)=
\begin{pmatrix}-7\\-6\end{pmatrix}$$
You can also do a second change of coodinate system, this time into one where the ellipse is not only axis-aligned but also centered at the origin.
Scaling
The completed squares above contributed some more terms to the constant factor $F$, since the equation by now has become
$$(x' + 9)^2 - 81 + 6.25(y' - 2)^2 - 25 + 81 = 0$$
The combined constant term here is $81 - 81 - 25 = -25$.
Now you move this to the right side of the equation, then divide the whole equation by this number so that you get the $=1$ from the desired form. Then you can deduce the radii $a$ and $b$.
\begin{align*}
(x' + 9)^2 + 6.25(y' - 2)^2 &= 25 \\
\frac{(x' + 9)^2}{25} + \frac{6.25(y' - 2)^2}{25} &= 1
\end{align*}
\begin{align*}
\frac1{a^2}&=\frac{1}{25} & a &= \sqrt{\frac{25}{1}}\phantom{.0} = 5 \\
\frac1{b^2}&=\frac{6.25}{25} & b &= \sqrt{\frac{25}{6.25}} = 2
\end{align*}
Verifying the result
Now let's check that we didn't make any mistakes. With the parameters we found, you can piece together the equation
\begin{align*}
\left(\frac{0.6(x-(-7))+0.8(y-(-6))}{5}\right)^2+
\left(\frac{0.6(y-(-6))-0.8(x-(-7))}{2}\right)^2 &= 1
\\
\frac{\bigl(0.6(x+7)+0.8(y+6)\bigr)^2}{25}+
\frac{\bigl(0.6(y+6)-0.8(x+7)\bigr)^2}{4} &= 1
\\
4\bigl(0.6(x+7)+0.8(y+6)\bigr)^2+
25\bigl(0.6(y+6)-0.8(x+7)\bigr)^2 - 100 &= 0
\\ &\;\vdots \\
17.44\,x^2 + 11.56\,y^2 - 20.16\,xy + 123.2\,x - 2.4\,y + 324 &= 0
\\
4.36\,x^2 + 2.89\,y^2 - 5.04\,xy + 30.8\,x - 0.6\,y + 81 &= 0
\end{align*}
So this is indeed the right equation, equivalent to the one given in the original example.
Generalizations
With a few sign changes, the above works for hyperbolas as well.
The approach using eigenvectors will generalize to ellipsoids in space. The approach using the form of the rotation matrix is harder in space.
Best Answer
Translate everything by $-c$ so that the circle center goes to the origin and $p$ goes to $p'$. Now take the cross product of the vectors $n$ and $p'$, let $q$. Normalize $p'$ and $q$ to be unit vectors. Now the matrix $R$ formed by $p',q,n$ is a rotation matrix that applies the plane $oxy$ to the plane $op'q$ and the axis $ox$ to $op'$.
Take the transpose of this orthogonal matrix, which is its inverse and perform the inverse rotation.
You now have a circle in the plane $oxy$, centered at the origin, and a point $p''$ on the $ox$ axis. By trigonometry, it is not a big deal to find the tangency points. (Check https://en.wikipedia.org/wiki/Tangent_lines_to_circles.)
Apply the rotation followed by the reverse translation to obtain the coordinates in the original system of axis.