A more general proof:
Let Q and R be the points at which lines through $P=(x_1,y_1)$ touch a non degenerate conic $S(x,y) \equiv Ax^2+2Bxy+Cy^2+2Dx+2Ey+F=0$. In other words, lines PR and PQ are the tangents to this conic at points Q and R, and RQ is the chord of contact.
Let $PR(x,y)=0$, $PQ(x,y)=0$, $RQ(x,y)=0$ be the equation of these lines.
As RQ is polar of P in relation to this conic,
$$RQ(x,y)\equiv (Ax+By+D)x_1+(Bx+Cy+E)y_1+(Dx+Ey+F)=0$$
On the other hand, the equation $\lambda(PR(x,y).PQ(x,y))+\mu(RQ(x,y))^2=0$ represents all conics which are touched by lines PR and PQ at points R and Q. Therefore, for especific values of $\lambda$ and $\mu $ (none of which can be equal to zero, because otherwise S would be a degenerate conic):
$$S(x,y)\equiv \lambda(PR(x,y).PQ(x,y))+\mu(RQ(x,y))^2=0$$
Then,
$$S(x_1,y_1)=\lambda(PR(x_1,y_1).PQ(x_1,y_1))+\mu(RQ(x_1,y_1))^2,$$
$$S(x_1,y_1)=\mu(RQ(x_1,y_1))^2$$
Besides that,
$$RQ(x_1,y_1)=(Ax_1+By_1+D)x_1+(Bx_1+Cy_1+E)y_1+(Dx_1+Ey_1+F),$$
$$RQ(x_1,y_1)=S(x_1,y_1)$$
Thus
$$S(x_1,y_1)=\mu(S(x_1,y_1))^2,$$
$$\mu=\frac {1}{S(x_1,y_1)}$$
Therefore
$$S(x_1,y_1).S(x,y)\equiv S(x_1,y_1)\lambda(PR(x,y).PQ(x,y))+(RQ(x,y))^2,$$
$$S(x_1,y_1)\lambda(PR(x,y).PQ(x,y))\equiv S(x_1,y_1).S(x,y)-(RQ(x,y))^2$$
Finally, equating left and right members of this identity to zero, we get that the equation of tangents PR and PQ to conic S can be represented by equation
$$S(x_1,y_1).S(x,y)-(RQ(x,y))^2=0$$
The standard circle form is $(x-1)^2+(y-1)^2=r^2=(2\sqrt2)^2$. Let $O(1,1)$, $A(a,a)$ and the subtended angle $\theta$. Then, $|AO|=\sqrt{2(a-1)^2}=\pm\sqrt2(a-1)$ and
$$ \sin \frac{\theta}2 = \frac r{|AO|} = \frac{2\sqrt2}{\pm \sqrt2(a-1)} =\pm \frac2{a-1}$$
Given that
$$ \frac12 =\sin\frac\pi6< \sin \frac{\theta}2 < \sin \frac\pi2 =1$$
we have
$$ \frac12 < \pm \frac2{a-1}<1$$
which leads to the ranges $a\in (-3,-1)\cup (3,5)$.
Best Answer
Here are a few other ways to find the tangent lines via (more or less) direct computation. They might be overkill when working with circles, but these methods generalize to other types of conics.
Solve a simpler problem and transform:
It’s a fairly simple matter to find the tangents to the unit circle through the point $(d,0)$: You can use either the Pythagoran theorem and similar triangles, or use the fact that the polar of this point is the line $x=\frac1d$ to find that the points of tangency are $\left(\frac1d,\pm{\sqrt{d^2-1}\over d}\right)$, and the corresponding lines are $x\mp y\sqrt{d^2-1}=d$. The general problem can be transformed into this one via translation, scaling and rotation, so the solution to the general problem can be obtained by transforming these lines back into the original coordinate system.
If we represent points as homogeneous column vectors, we can also represent lines as homogeneous row vectors, so that the equation $ax+by+c=0$ of a line can be written as $\mathscr l\mathbf x = [a:b:c]\,[x:y:1]^T = 0$. Given the nonsingular point transformation $\mathbf x'=M\mathbf x$, we have $$\mathscr l\mathbf x=\mathscr l(M^{-1}M\mathbf x)=(\mathscr lM^{-1})(M\mathbf x)=0,$$ which shows that lines are covariant and transform as $\mathscr l'=\mathscr lM^{-1}$. So, to convert the unit circle tangent lines into the original coordinate system, apply the three transformations in reverse order, i.e., if $M=RST$ is the transformation that maps points into the unit circle’s coordinate system, then the tangents are $\mathscr lM = [1:\mp\sqrt{d^2-1}:d]RST$. For the circle $(x-h)^2+(y-k)^2=r^2$ and point $\mathbf p=[x_0:y_0:1]^T$, this expands into $$\begin{bmatrix}1&\mp\sqrt{d^2-1}&-d\end{bmatrix} \begin{bmatrix}{\xi\over\sqrt{\xi^2+\eta^2}} & {\eta\over\sqrt{\xi^2+\eta^2}}&0 \\ -{\eta\over\sqrt{\xi^2+\eta^2}} & {\xi\over\sqrt{\xi^2+\eta^2}} & 0 \\ 0&0&1\end{bmatrix} \begin{bmatrix}\frac1r&0&0\\0&\frac1r&0\\0&0&1\end{bmatrix} \begin{bmatrix}1&0&-h\\0&1&-k\\0&0&1\end{bmatrix},$$ where $(\xi,\eta)=(x-h,y-k)$ and $d=\frac1r\sqrt{\xi^2+\eta^2}$, i.e., the $x$-coordinate of $RST\mathbf p$. Since these are homogeneous matrices, you can simplify your computations a bit by multiplying $R$ through by $\sqrt{\xi^2+\eta^2}$ and using $S=\operatorname{diag}(1,1,r)$ instead. This gives us $$M = \begin{bmatrix}\xi&\eta&-h\xi-k\eta \\ -\eta&\xi&h\eta-k\xi \\ 0&0&r\sqrt{\xi^2+\eta^2}\end{bmatrix}.$$ Similarly, you can multiply the lines by $r$, producing the equivalent lines $[r:\mp\sqrt{\xi^2+\eta^2-r^2}:-\sqrt{\xi^2+\eta^2}]$.
For your example, we have $h=3$, $k=2$, $r=3$ and $\mathbf p=[4:-5:1]$, giving $[3:\pm\sqrt{41}:-5\sqrt2]$ for the unit circle tangents and finally $$\begin{bmatrix}3 & \pm\sqrt{41} & -5\sqrt2\end{bmatrix} \begin{bmatrix}1&-7&11\\7&1&-23\\0&0&15\sqrt2\end{bmatrix} = \begin{bmatrix}3\pm7\sqrt{41} & -21\pm\sqrt{41} & -117\mp23\sqrt{41}\end{bmatrix}$$ as the solution to the problem. These lines and the circle are plotted below:
This method works for ellipses as well if you use a non-uniform scaling and modify $\xi$ and $\eta$ accordingly. If the ellipse isn’t axis-aligned, then you might need to add another rotation before scaling, but the basic method is the same.
Degenerate conic:
If you don’t need the equations of the individual tangent lines, then you can produce a degenerate conic that consists of this pair of lines quite directly. This works for any conic. Let $C$ be its homogeneous matrix and $\mathbf p$ be the point through which the tangents are to be drawn. The dual conic $C^*$ consists of all of the lines tangent to $C$, and $[\mathbf p]_\times^T C^*[\mathbf p]_\times$ consists of only the tangents through $\mathbf p$. Here, $[\mathbf p]_\times$ is the “cross product matrix:” if $\mathbf p = [x:y:z]$, then $$[\mathbf p]_\times = \begin{bmatrix}0&-z&y\\w&0&-x\\-y&x&0\end{bmatrix}.$$ The dual conic matrix $C^*$ is the adjugate of $C$, or, if $C$ is nonsingular, $C^{-1}$. For your example, $$[\mathbf p]_\times = \begin{bmatrix}0&-1&-5\\1&0&-4\\5&4&0\end{bmatrix}$$ and $$C=\begin{bmatrix}1&0&-3\\0&1&-2\\-3&-2&4\end{bmatrix},$$ so, using the adjugate, $$[\mathbf p]_\times^T C^* [\mathbf p]_\times = \begin{bmatrix}40&7&-125\\7&-8&-68\\-125&-68&160\end{bmatrix}.$$ After multiplying it out and pulling out a factor of two, this gives the equation $20x^2+7xy-4y^2-125x-68y+80=0$. If you do need the individual lines, you can try to factor this equation or use the method described here to split the degenerate conic. The latter involves computing the adjugate of $[\mathbf p]_\times^T C^* [\mathbf p]_\times$, among other things, but this for this problem this is extra work because you can instead compute the...
Dual to the intersection of line and conic:
Raffaele’s answer includes a method that computes the intersection of the polar line of the point through which tangents are to be drawn with the circle to find the points of tangency, from which you can compute the equations of the tangent lines. If you don’t need those points, however, you can compute the lines directly by working with the dual conic: the intersections of this conic with the given point are the tangent lines through that point. For the circle $$C = \begin{bmatrix}1&0&-h\\0&1&-k\\-h&-k&h^2+k^2-r^2\end{bmatrix}$$ the dual conic is given by its adjugate $$C^* = \begin{bmatrix}h^2-r^2&hk&h\\hk&k^2-r^2&k\\h&k&1\end{bmatrix}.$$ Finding the tangent lines is then a matter of solving the system $$[\lambda:\mu:\tau]\mathbf p=0 \\ [\lambda:\mu:\tau]C^*[\lambda:\mu:\tau]^T=0.$$ You can also compute it directly using the method outlined here.
From above, we have $$[\mathbf p]_\times^TC^*[\mathbf p]_\times = \begin{bmatrix}40&7&-125\\7&-8&-68\\-125&-68&160\end{bmatrix}.$$ We need to find a coefficient $\alpha$ such that $[\mathbf p]_\times^TC^*[\mathbf p]_\times+\alpha [\mathbf p]_\times$ has rank $1$. There are a couple of ways to do this, and one of them produces $\alpha = \sqrt{-\tiny{\begin{vmatrix}40&7\\7&-8\end{vmatrix}}} = 3\sqrt{41}$ for $$[\mathbf p]_\times^TC^*[\mathbf p]_\times+\alpha[\mathbf p]_\times = \begin{bmatrix}40&7-3\sqrt{41}&-125-15\sqrt{41} \\ 7+3\sqrt{41}&-8&-68-12\sqrt{41} \\ -125+15\sqrt{41}&-68+12\sqrt{41}&160\end{bmatrix}.$$ We can take for our pair of tangent lines any row and column of this matrix.