The trinomial
$$ax^2+bx+c=0$$
May be solved when $a,b,c$ are complex numbers, and it works fine with real $a,b,c$, since $\Bbb R\subset\Bbb C$.
However, we'll start with the real case.
Real resolution
We start with the equation
$$ax^2+bx+c=0$$
Where $a\neq0$, and $a,b,c\in\Bbb C$.
Then, rewrite the equation
$$\left(x+\frac{b}{2a}\right)^2-\frac{b^2-4ac}{4a^2}=0$$
Then, according to the sign of $\Delta=b^2-4ac$, there are three possibilities:
- If $\Delta=0$, the equation amounts to
$$\left(x+\frac{b}{2a}\right)^2=0$$
It has thus one double root, $x=-\frac{b}{2a}$.
But since $\Delta=0$, we can also write
$$x=\frac{-b\pm\sqrt{\Delta}}{2a}$$
- If $\Delta>0$, then $\frac{b^2-4ac}{4a^2}$ is positive, so it's the square of number $\frac{\sqrt{\Delta}}{2a}$, and we can write
$$\left(x+\frac{b}{2a}\right)^2-\left(\frac{\sqrt{\Delta}}{2a}\right)^2=0$$
$$\left(x+\frac{b}{2a}+\frac{\sqrt{\Delta}}{2a}\right)\left(x+\frac{b}{2a}-\frac{\sqrt{\Delta}}{2a}\right)=0$$
Hence
$$x=\frac{-b\pm\sqrt{\Delta}}{2a}$$
- If $\Delta<0$, then the equation has no real solution, since it would amount to
$$\left(x+\frac{b}{2a}\right)^2=\frac{\Delta}{4a^2}$$
Where the left-hand side is $\geq0$, whereas the right-hand side is $<0$.
However, we can still write that $\Delta$ is the square of complex number $i\sqrt{4ac-b^2}$, hence we can write
$$\left(x+\frac{b}{2a}\right)^2-\left(\frac{i\sqrt{-\Delta}}{2a}\right)^2=0$$
$$\left(x+\frac{b}{2a}+\frac{i\sqrt{-\Delta}}{2a}\right)\left(x+\frac{b}{2a}-\frac{i\sqrt{-\Delta}}{2a}\right)=0$$
$$x=\frac{-b\pm i\sqrt{-\Delta}}{2a}$$
But, with the convention that $\sqrt{t}=i\sqrt{-t}$ when $t<0$, we can again write
$$x=\frac{-b\pm\sqrt{\Delta}}{2a}$$
So, in all cases, the solutions are
$$x=\frac{-b\pm\sqrt{\Delta}}{2a}$$
Hence your equivalence holds.
Complex resolution
Again, we start with the equation
$$ax^2+bx+c=0$$
Where $a\neq0$, and $a,b,c\in\Bbb C$.
Then, rewrite the equation
$$\left(x+\frac{b}{2a}\right)^2-\frac{b^2-4ac}{4a^2}=0$$
Now, a complex number $z$ has always a square root (two distinct ones if $z\neq0$, see for example here). We go on with a square root $\alpha$ of $b^2-4ac$, that is, $\alpha^2=b^2-4ac$, so that the equation can be written
$$\left(x+\frac{b}{2a}\right)^2-\frac{\alpha^2}{4a^2}=0$$
It's a difference of squares, and we know that $A^2-B^2=(A+B)(A-B)$, so
$$\left(x+\frac{b}{2a}+\frac{\alpha}{2a}\right)\left(x+\frac{b}{2a}-\frac{\alpha}{2a}\right)=0$$
So the roots are
$$x=-\frac{b}{2a}\pm\frac{\alpha}{2a}$$
Or if you prefer
$$x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}$$
Important note: there is no canonical way to write that $\sqrt{t}$ is a specific square root of $t\in\Bbb C$, hence this symbol is meaningless. However, since it appears with a $\pm$, we will always work with both roots, so using the symbol here is harmless. Nevertheless, it should be avoided.
Best Answer
At heart, it appears that you’re reading a bit too much into the meaning of the sign of the discriminant of your equation $$(P_x^2+P_y^2)T_x^2-2R^2P_xT_x+R^2(R^2-P_y^2)=0. \tag{6b}$$ Let’s look at it a bit more closely. We have $$\Delta = 4 R^4 P_x^2-4R^2(P_x^2+P_y^2)(R^2-P_y^2) = 4 R^2 P_y^2 (P_x^2+P_y^2-R^2).$$ When $P_y\ne0$, the sign depends only on $P_x^2+P_y^2-R^2$, which in turn depends on whether $P$ is exterior to, interior to, or on the circle, as you’ve noted. However, when $P_y=0$ this term is masked: $\Delta = 0$ regardless of the relative position of $P$ vis-à-vis the circle, so in this case we can’t draw any conclusions about that from the discriminant. Backing up a bit, when $P_y=0$, your equation 4b is simply $P_xT_x=R^2$, which is the equation of a line, and squared is $(P_xT_x-R^2)^2=0$. The term $P_x^2+P_y^2-R^2$ doesn’t even appear in the discriminant of this equation, $4R^4P_x^2-4R^4P_x^2$.
The discriminant hasn’t failed. It tells you that the equation has only one solution, just as it should. The failure, if there is one, lies in the assumption that all of the solutions of equation 6b are valid $x$-coordinates of the points of tangency. This is clearly not the case when $P_y=0$. Equation 6b can overgenerate potential solutions, just as does plugging computed values of $T_x$ into equation 1a and solving for $T_y$. As zahbaz points out, when $P_y=0$ and $P_x \lt R$, we have $T_x\gt R$ and there are no real values of $T_y$ that satisfy equation 1a, so all is well.
This brings up an important caveat: the solutions to a system of equations also satisfy derived equations, but not, in general, vice-versa. This shouldn’t be much of a surprise. After all, not all of the solutions to the original individual equations satisfy the system, either.
That aside, I’ve got a few other notes on your derivation. I think you might be a bit hasty in solving for $T_y$ right away. If nothing else, you’ve potentially discarded a solution right off the bat. In the same vein, when you square both sides of the equation in step 5, you risk introducing extraneous solutions. In this case, I think these two operations more or less cancel each other out.
Equation 2c, which captures the condition $\overrightarrow{OT} \perp \overrightarrow{TP}$, describes the circle with diameter $OP$, and equation 1a is of course the equation of the original circle, so your method essentially computes the intersection of two circles. If you were to instead simply add these two equations you get $$P_xT_x+P_yT_y=R^2.\tag{*}$$ This is the equation of the radical line of the two circles, which passes through their intersection points, if any. (As an illustration of the earlier caveat, at most two solutions to this equation satisfy the original system; most don’t.) The problem is then reduced to finding the intersection of this line with either of the two circles. This can be done with a straightforward substitution for one of $T_x$ or $T_y$ into 1a, and then using the quadratic formula to solve for the other. Substitute back into (*) to find the corresponding values of the other variable. Because this equation is linear, doing so won’t overgenerate potential solutions the way back-substituting into either circle equation does. If $P_x=0$ or $P_y=0$, equation (*) gives you the value of one of $T_x$ or $T_y$ directly, which you can then substitute into 1a to get the corresponding values of the other coordinate.
Since it looks like you eventually want the tangent lines themselves instead of just the points of tangency, I’ll briefly discribe a way to compute them directly using homogeneous coordinates. The equation of any conic can be written as $\mathbf x^TC\mathbf x=0$, where $C$ is a symmetric $3\times3$ matrix. For any nondegenerate conic, every tangent line $\mathbf l$ to it satisfies the dual equation $\mathbf l^TC^{-1}\mathbf l=0$. Now, the line through a pair of points can be represented by their cross product, so for a fixed point $\mathbf p$, the equation $(\mathbf p\times\mathbf x)^TC^{-1}(\mathbf p\times\mathbf x)=0$ is satisfied by all points that lie on a tangent to $C$ through $\mathbf p$. If $M_{\mathbf p}$ is the “cross-product matrix” of $\mathbf p$, i.e., if $M_{\mathbf p}\mathbf x=\mathbf p\times\mathbf x$ for all $\mathbf x$, this means that the matrix $A = M_{\mathbf p}^TC^{-1}M_{\mathbf p}$ represents the degenerate conic that consists of the tangents to $C$ through $\mathbf p$. If these lines are $\mathbf l$ and $\mathbf m$, then this matrix is a multiple of $\mathbf l\mathbf m^T-\mathbf m\mathbf l^T$. We can “split” this conic by finding a value of $\alpha$ such that $B = A-\alpha M_{\mathbf p}$ has rank one. The matrix $B$ is a multiple of either $\mathbf l\mathbf m^T$ or $\mathbf m\mathbf l^T$, from which we can read the two lines. Basically, this is a way to factor the degenerate conic equation into a product of linear terms.
Omitting the details of the calculation, one of the matrices that you might end up with for this particular problem is $$\begin{bmatrix} R^2-P_y^2 & P_xP_y + R\sqrt{P_x^2+P_y^2-R^2} & -R\left(RP_x+P_y\sqrt{P_x^2+P_y^2-R^2}\right) \\ P_xP_y - R\sqrt{P_x^2+P_y^2-R^2} & R^2-P_x^2 & -R\left(RP_y-P_x\sqrt{P_x^2+P_y^2-R^2}\right) \\ -R\left(RP_x-P_y\sqrt{P_x^2+P_y^2-R^2}\right) & -R\left(RP_y+P_x\sqrt{P_x^2+P_y^2-R^2}\right) & R^2(P_x^2+P_y^2) \end{bmatrix}.$$ Any row/column pair of $B$ that has a common nonzero diagonal element gives us the two tangent lines. Assuming that $P\ne0$, we can take the last row and column to get, after a bit of tweaking, the equations $$\left(P_x \pm \frac{P_y}R\sqrt{P_x^2+P_y^2-R^2}\right)x + \left(P_y \mp \frac{P_x}R\sqrt{P_x^2+P_y^2-R^2}\right)y = P_x^2+P_y^2.$$ As is evident from the quantity under the radical, there are two, one or no tangent lines when $P$ is exterior to, on or interior to the circle, respectively.