Simplifying an expression to a sensible conic section polynomial

Question

I have an expression which represents an intersection of the unit sphere and a cone, projected to two-dimensional plane:

expr = x^2 + y^2 <= 1 && 
   1/Sqrt[5] (2 (1 + Sqrt[5]) x^2 + (-1 + Sqrt[5]) x (-2 y + 
          Sqrt[2 (5 + Sqrt[5])] Sqrt[1 - x^2 - y^2]) + 
       y ((3 + Sqrt[5]) y + 
          2 Sqrt[2 (5 + Sqrt[5])] Sqrt[1 - x^2 - y^2])) == 1;

$x^2+y^2\leq 1\land \frac{\left(\sqrt{5}-1\right) x \left(\sqrt{2 \left(5+\sqrt{5}\right)} \sqrt{1-x^2-y^2}-2 y\right)+y \left(2 \sqrt{2 \left(5+\sqrt{5}\right)} \sqrt{1-x^2-y^2}+\left(3+\sqrt{5}\right) y\right)+2 \left(1+\sqrt{5}\right) x^2}{\sqrt{5}}=1$

On basis of its geometric origins it is obviously an ellipse which lies inside the unit disk:

It has the $\sqrt{1-x^2-y^2}$ term which I find hard to eliminate, at least just playing with assumptions.

I want to rewrite this equation to its polynomial, conic section form. How to accomplish this without resorting to numerics?

---

My best effort solution still involves numerics through combination of high-precision N and RootApproximant, which works at least with these coefficients.

a x^2 + b x y + c y^2 + d x + e y == 1 /.
  (First@Solve[
      a x^2 + b x y + c y^2 + d x + e y == 1 /.
       FindInstance[expr, {x, y}, Reals, 5]] /.
    v_?NumericQ :> RootApproximant[N[v, 100]]) //
 FullSimplify

Here the conic section solution is calculated by applying five solutions to expr with FindInstance to the conic section equation. Solve'ing this problem does produce exact coefficients, but they're over ten megabytes in size which makes simplification hopeless. The "simplification" through the aforementioned numerical route does produce a sensible result, though:

$\frac{8 x^2}{\sqrt{5}}+\left(2-\frac{2}{\sqrt{5}}\right) x y+\sqrt{2} \left(\sqrt{5}-1\right) x+y \left(\frac{7 y}{\sqrt{5}}+y+2 \sqrt{2}\right)=1$

The positive side is that validity of this solution can be actually checked (although this takes a moment):

Resolve[ForAll[{x, y}, Equivalent[%, expr]], Reals]

(* True *)

Are there better, purely symbolic ways to accomplish this based on expr?

---

EDIT:

One may look at those FindInstance results and conclude that the roots are horribly complicated even after RootReduce; what if one would just pick "simple" numbers and find solutions restricted by, say, $x=0$, $y=0$ and some other very simple line. Yes, roots are simpler, especially after RootReduce, and this probably would help Solve. The problem is that even for $x=0$ and $y=0$ finding real instances takes minutes, and something like $x=-\frac{1}{2}$ clearly a lot longer.

Answer 1

This is not a general-purpose answer but it works in this case at least. First step is to obtain an implicit polynomial in {x,y} for that second expression. Often enough, GroebnerBasis can do this. It makes internal variables and polynomial relations to handle radicals such as Sqrt[1 - x^2 - y^2]. The first element in the result should be in terms of only the desired variables.

xypoly = 
 First[GroebnerBasis[
   1/Sqrt[5] (2 (1 + Sqrt[5]) x^2 + (-1 + Sqrt[5]) x (-2 y + 
          Sqrt[2 (5 + Sqrt[5])] Sqrt[1 - x^2 - y^2]) + 
       y ((3 + Sqrt[5]) y + 
          2 Sqrt[2 (5 + Sqrt[5])] Sqrt[1 - x^2 - y^2])) - 1, {x, y}]]

(* Out[127]= 5 + (-60 + 4 Sqrt[5]) x^2 + 
 64 x^4 + (20 - 36 Sqrt[5]) x y + (-32 + 32 Sqrt[5]) x^3 y + (-50 - 
    14 Sqrt[5]) y^2 + (136 + 8 Sqrt[5]) x^2 y^2 + (-8 + 
    24 Sqrt[5]) x y^3 + (54 + 14 Sqrt[5]) y^4 *)

Plot this quartic.

ContourPlot[xypoly == 0, {x, -1, 1}, {y, -1, 1}]

So clearly a product. Separating the factors, now that's a challenge.

First we try to factor using the "obvious" algebraic extension Sqrt[5].

Factor[First[gb], Extension -> {Sqrt[5]}]

(* Out[128]= 5 + (-60 + 4 Sqrt[5]) x^2 + 
 64 x^4 + (20 - 36 Sqrt[5]) x y + (-32 + 32 Sqrt[5]) x^3 y + (-50 - 
    14 Sqrt[5]) y^2 + (136 + 8 Sqrt[5]) x^2 y^2 + (-8 + 
    24 Sqrt[5]) x y^3 + (54 + 14 Sqrt[5]) y^4 *)

No joy, in or out of Mudville. Let's make sure the thing really is a product.

IrreduciblePolynomialQ[xypoly, Extension -> All]

(* Out[130]= False *)

So we really do need to try harder. Here we take a guess that all algebraics needed come from square roots of integer coefficient factors. All such are accounted for in the set {2,3,5}.

facs = Factor[xypoly, Extension -> {Sqrt[5], Sqrt[2], Sqrt[3]}]

(* Out[132]= -((Sqrt[5] + (-5 Sqrt[2] + Sqrt[10]) x - 8 x^2 - 
     2 Sqrt[10] y + (2 - 2 Sqrt[5]) x y + (-7 - Sqrt[5]) y^2) (-Sqrt[
      5] + (-5 Sqrt[2] + Sqrt[10]) x + 8 x^2 - 
     2 Sqrt[10] y + (-2 + 2 Sqrt[5]) x y + (7 + Sqrt[5]) y^2))
*)

So we found a pair of factors. I leave to the interested reader how to determine which is correct. The issues with the method thus far are that factoring over algebraic extensions can be slow, and when degrees are higher, finding the right extension can be tricky. So I would not expect this methodology to scale very well.

Answer 2

You can solve the second equation for y:

eq = 1/Sqrt[
      5] (2 (1 + Sqrt[5]) x^2 + (-1 + Sqrt[5]) x (-2 y + 
         Sqrt[2 (5 + Sqrt[5])] Sqrt[1 - x^2 - y^2]) + 
      y ((3 + Sqrt[5]) y + 
         2 Sqrt[2 (5 + Sqrt[5])] Sqrt[1 - x^2 - y^2])) == 1;

sol[x_] = y /. Solve[eq, y];

You can then plot the solutions together with the unit circle (first equation):

Plot[Evaluate@sol[x], {x, -1, 1}, AspectRatio -> Automatic, 
 Epilog -> Circle[]]

Answer 3

Here's a hack around FindInstances returning complicated roots. simpleInstances finds three lines aligned with each axis which intersect with the region defined by expr and have a "simple" rational form, that is sum of absolute numerator and denominator is minimised on one coordinate. Found instances of expr are RootReduced and five with smallest LeafCount are returned.

With this code the result from main Solve is still tolerable for symbolic simplification by FullSimplify, and returns the same answer as in the question:

ClearAll[simpleInstances];
simpleInstances[expr_] :=
 MapThread[
    Function[{axis, bounds},
     Table[
      Solve[expr && axis == val, {x, y}],
      {val,
       (* Find 3 simple rationals inside bounds. *)
       Nest[
        Append[#,
          n/d /.
           Last@Minimize[
             Abs[n] + d,
             Join[
              {bounds[[1]] < n/d < bounds[[2]], d >= 1},
              n/d != # & /@ #],
             {n, d}, Integers]] &, {}, 3]}]],
    {{x, y},
     RegionBounds[ImplicitRegion[expr, {x, y}]]}] //
   # /. r_Root :> RootReduce[r] & //
  TakeSmallestBy[Flatten[#, 2], LeafCount, 5] &

a x^2 + b x y + c y^2 + d x + e y == 1 /.
 (First@Solve[a x^2 + b x y + c y^2 + d x + e y == 1 /.
      simpleInstances[expr], Reals] // FullSimplify)

$$\frac{8 x^2}{\sqrt{5}}+\left(2-\frac{2}{\sqrt{5}}\right) x y+\sqrt{2} \left(\sqrt{5}-1\right) x+\left(1+\frac{7}{\sqrt{5}}\right) y^2+2 \sqrt{2} y=1$$

In this case it's not a problem, but some of the found solution instances might actually be duplicates, and that should be handled in a more general case.