It seems you're homogenizing by introducing a first coordinate equal to one, instead of a last coordinate conventionally used. So your conic is essentially described as
$$
(1,x,y)\cdot\begin{pmatrix}
\tfrac12 & 1 & \tfrac12 \\
1 & 1 & 3 \\
\tfrac12 & 3 & 1
\end{pmatrix}\cdot\begin{pmatrix}1\\x\\y\end{pmatrix} = 0
$$
I'll stick to your notation for this post, but readers should be careful if translating between this and a different convention.
Tangents
The asymptotes are the tangents at the points at infinity. So first you intersect the conic with the line at infinity, e.g. by plugging a generic point at infinity into the equation:
$$
0 = (0,x,y)\cdot\begin{pmatrix}
\tfrac12 & 1 & \tfrac12 \\
1 & 1 & 3 \\
\tfrac12 & 3 & 1
\end{pmatrix}\cdot\begin{pmatrix}0\\x\\y\end{pmatrix} = x^2+6xy+y^2
$$
You can check that neither $x=0$ nor $y=0$ is a solution. As the result is only defined up to a scalar factor, you may assume $y=1$ and solve for $x$:
$$x_{1,2} = \pm\sqrt8-3$$
Now you compute the tangent, which is simply the matrix times the point at infinity.
$$g_{1,2} =
\begin{pmatrix}
\tfrac12 & 1 & \tfrac12 \\
1 & 1 & 3 \\
\tfrac12 & 3 & 1
\end{pmatrix}\cdot\begin{pmatrix}0\\\pm\sqrt8-3\\1\end{pmatrix}
=\begin{pmatrix}\pm\sqrt8-\tfrac52\\\pm\sqrt8\\\pm\sqrt{72}-8\end{pmatrix}
$$
Now you can read this as the equations of a pair of lines:
$$(\pm\sqrt8-\tfrac52) + (\pm\sqrt8)x + (\pm\sqrt{72}-8)y = 0$$
Note that the principal axes (the ones you computed using eigenvectors) are the perpendicular bisectors of these asymptotes. Computing the axes from the asymptotes is possible, but the converse is impossible as the axes alone don't contain sufficient information.
Foci
For foci you can do the following:
- Compute the polar lines of the ideal circle points which are $(0,1,\pm i)$ in your coordinates. That's the imaginary unit $i$ there, so you need to perform the computation using complex numbers. The polar line is matrix times the coordinates of a point, so in the special case of a point that lies on the conic the polar will be the tangent.
- Intersect these polar lines with the conic. That requires a conic-line intersection similar to the one above, except that you can't simply assume that one of the coordinates is zero. But it is still just a quadratic equation in a single parameter. Or you can follow this approach instead.
- Connect the points of intersection with the circle point from which you started, to get the tangents from this point. So here you have two ideal circle points, and from each of them you have two tangents to the conic.
- Intersect these lines, taking one line for each of the ideal circle points. There are four ways to intersect them, resulting in four foci. Two of them will be real, two complex. Finding the pairings is simple: if you have one real focus, then switching to the other tangent for each of the circle points will result in the other real focus.
This approach follows J. Richter-Gebert's Perspectives on Projective Geometry Section 19.4.
For your input you get
- For $I_1$ = $(0,1,i)$ you get $h_1:(1+\tfrac12i)+(1+3i)x+(3+1i)y=0$
and for $I_2$ = $(0,1,-i)$ you get $h_2:(1-\tfrac12i)+(1-3i)x+(3-1i)y=0$.
- For $h_1$ you get $P_{11}=(32,-2-(1-3i)\sqrt{6i},-10+(3-i)\sqrt{6i})$
and $P_{12}=(32,-2+(1-3i)\sqrt{6i},-10-(3-i)\sqrt{6i})$.
For $h_2$ you get $P_{21}=\overline{P_{11}}$ and $P_{22}=\overline{P_{12}}$ since $h_2=\overline{h_1}$.
- Connecting $I_a$ to $P_{ab}$ gives $t_{ab}$ like this:
$\begin{array}{rr}
t_{11}:& (10-2i-6\sqrt{6i}) - 32i\,x + 32\,y = 0 \\
t_{12}:& (10-2i+6\sqrt{6i}) - 32i\,x + 32\,y = 0 \\
t_{21}:& (10+2i-6\sqrt{-6i}) + 32i\,x + 32\,y = 0 \\
t_{22}:& (10+2i+6\sqrt{-6i}) + 32i\,x + 32\,y = 0
\end{array}$
- Intersecting $t_{11}$ with $t_{21}$ yields one real focus at
$F_1=(16,-1-3\sqrt3,-5+3\sqrt3)$,
so $t_{12}$ with $t_{22}$ yield the other at $F_2=(16,-1+3\sqrt3,-5-3\sqrt3)$.
The other foci form a pair of complex conjugate points at $F_{3,4}=(16,-1\pm3i\sqrt3,-5\pm3i\sqrt3)$.
Note that I just wrote $16$ in the first coordinate to avoid fractions. If you don't like working with homogeneous coordinates, simply divide by this number and you get plain coordinates $x$ and $y$ in the second and third coordinate. The results are $F_1\approx(-0.387,0.012)$ and $F_2\approx(0.262,-0.637)$.
