Simplifying an expression to a sensible conic section polynomial

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.

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.

• One can generate five solutions to the original expression with FindInstances and assign them to factors above, and the factor which evaluates to zero for all of them (with help from RootReduce) is the correct one. But yes, choosing a suitable factoring extension probably makes this solution approach more or less manual... Feb 27 at 4:25
• ... or one can Resolve which factor is Equivalent with the original expression. That's likely to be slower, though. Feb 27 at 12:02
• Related: One can pick rational (say) values for x, solve for y, and interpolate the correct quadratic implicit polynomial. This might avoid the guesswork and potentially slow factorization. Feb 27 at 15:14
• …which I (belatedly) see is similar to what was already done in the response by @kirma a couple of days ago. Feb 28 at 13:54

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[]]


• I'm positively surprised that Solve finds implicit solutions. There are three issues, though: 1. these are implicit solutions (maybe I wasn't, ehm, explicit enough that I'm interested on the implicit polynomial solution), 2. they come in four parts, and 3. two of these four solutions are not actually real solutions (but at least Solve warns of this...). Feb 26 at 11:22
• Ehm, I meant that Solve doesn't find purely polynomial solution, like $A x^2+B x y+C y^2+D x+E y=1$. Feb 26 at 11:35
• There are 4 solutions because , as a function of x, there are 4 single valued functions. To get a polynomial solution, you write a solution: y= ... Sqrt[..], isolate the Sqrt on one side and square. Feb 26 at 11:35
• Assigning solution instances from green and red curves (apart from intersections with the orange curve) to the original equation actually evaluate False. Feb 26 at 11:44
• I am not sure what is going on here. It would be interesting what Wolfram has to say about this. Further, are 2 equal ellipses feasible for your setup of the cone and sphere? Feb 26 at 15:02

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_] :=
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.