This is what happens when misdirected comics discuss intersected conics.

That's useful to know. I guess it's 'cos FindInstance is not able to handle an overabundance of sin.

There are various ways to get information about real root constraints. I will show one that is more or less feasible, although the result might be too large to be of use and even so the result is incomplete.

The idea is to determine where the zeros hit multiple roots (this is a part of what is called,perhaps not surprisingly, the "discriminant variety"). It is at such values that pairs cross to/from two real values to complex conjugates. This happens precisely at points where the ellipses touch tangentially (found with a Lagrange multiplier). In order to make this viable I switched from your formulation as two explicit ellipses to a pair of "generic" guadratics.

quadratics = {a*x^2 + b*x*y + c*y^2 + d*x + e*y + f, 
   p*x^2 + q*x*y + r*y^2 + s*x + t*y + u};
grad1 = Grad[quadratics[[1]], {x, y}];
grad2 = Grad[quadratics[[2]], {x, y}];
polys = Flatten[{quadratics, grad1 - lambda*grad2}];
elim = {x, y, lambda};

Timing[gb = 
   GroebnerBasis[polys, Complement[Variables[polys], elim], elim, 
    MonomialOrder -> EliminationOrder];]

(* {2.377909, Null} *)

It contains one polynomial, and that polynomial (in the various parameters) vanishes precisely where we have the transition points. So the vanishing gives the condition for where these crossings happen. I will note that the polynomial has a leafcount in the tens of thousands.

Another way to approach this might be to mess with resultants and discriminants. A resultant of the two quadratics, say with respect to x, eliminates x and gives a (parametrized) univariate in y. The discriminant of this new polynomial, with respect to y of course, should then also provide the vanishing conditions. But this step is slow, most likely because it also spawns what are called spurious factors.

If one really wants inequality conditions in the parameter space that divide the regions of zero, two, or four real solutions (counting multiplicity), that requires CylindricalDecomposition. I doubt any formulation will run in reasonable time or produce a result of less than too-stupendous-to-use size.

Did I really chase you to ellipses, or was that just, err, hyperbole?

My example is (what a surprise) finding out if two ellipses overlap. As before, my equation for an ellipse with semi-major axes, a and b, centered at {xc ,yc} and rotated by an angle theta is eleq[{a,b,xc,yc,theta}][{x,y}] == 0 where

eleq[{a_, b_, xc_, yc_, \[Theta]_}][{x_, 
   y_}] = (((x/a)^2 + (y/b)^2 - 1 ) /. {x -> 
      Cos[\[Theta]]*x + Sin[\[Theta]] y, 
     y -> -Sin[\[Theta]]*x + Cos[\[Theta]]*y}) /. {x -> x - xc, 
   y -> y - yc}

-1 + ((y - yc) Cos[\[Theta]] - (x - xc) Sin[\[Theta]])^2/b^2 + ((x - 
      xc) Cos[\[Theta]] + (y - yc) Sin[\[Theta]])^2/a^2

If the two ellipses do not overlap, all 4 solutions to { eleq[{a1,b1,xc1, yc1, theta1}][{x,y}] == 0, eleq[{a2,b2,xc2,yc2,theta2}]{x,y}]==0} will be complex. What I've done so far is to ComplexExpand eleq to get two equations in the real and imaginary parts of x and y and test that they give the same solutions for pairs of ellipses as the original equations. Not sure how to proceed from there. Daniel, your concern about my obsession with packing circles has caused me to move on to ellipses, so it would be good if you can help out.

Attachments: