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.
Best Answer
You can switch over to the modified solution formula, according to binomial formulas, $$ x=\frac{-b+\sqrt{b^2-4ac}}{2a}=-\frac{2c}{b+\sqrt{b^2-4ac}}. $$ The cancellation case in the first formula will give a numerically stable computation in the second formula, giving a value close to $-\frac{c}{b}$, the solution of the linear approximation.
In general, most of the cases with numerical issues when computing the roots, in the case of real roots, are avoided by first determining the non-cancelling variant of the numerator $$ u=-b-{\rm sign}^+(b)\sqrt{b^2-4ac},~~~~{\rm sign}^+(0)=+1, $$ and then computing the two roots as $$ x_1=\frac{u}{2a},~~~ x_2=\frac{2c}{u}. $$ After that, the stabilization discussion turns to the square root, to be able to keep computing it in the case of floating-point over- or underflow.