# Ideal.variety() when working with symbolic rings

I have a system of multivariate polynomials I would like to solve (the equations come from chemical reaction networks). I have many symbolic constants, but I've managed to generate a Groebner basis with the following code (using a toy example with solutions {{x: x0, y:0},{x: 0,y: x0}):

```
x0 = var('x0') # In practice, there are many constants
S.<x,y> = PolynomialRing(SR,order='lex'); # in practice, there are 3-10 variables
I = S.ideal(
x**2 + y**2 - x0**2,
x + y - x0
)
print(I.dimension()) # dimension should always be zero
# >> 0
print(I.groebner_basis())
# >> [x + y - x0,
# >> y^2 + (-x0)*y]
```

In practice, I get many resulting equations which are extremely long, and it is numerically challenging to sequentially find the roots of these polynomials (in downstream code, I will get values for the constant parameters and then need to return solutions to these equations on the fly). Using the `I.variety()`

method would be great to symbolically isolate all possible solutions, but I get confusing internal errors when I try to evaluate it.

```
print(I.variety()) # should be equivalent to: I.variety(ring=SR)
#/Users/ksb/SageMath/local/lib/python2.7/site-packages/sage/rings/polynomial/multi_polynomial_ideal.pyc in _variety(T, V, v)
# 2340
# 2341 variable = f.variable(0)
#-> 2342 roots = f.univariate_polynomial().roots(ring=ring, multiplicities=False)
# 2343
# 2344 for root in roots:
#
#/Users/ksb/SageMath/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_element.pyx in sage.rings.polynomial.polynomial_element.Polynomial.roots (build/cythonized/sage/rings/polynomial/polynomial_element.c:68811)()
# 7714 return l
# 7715 else:
#-> 7716 return [val for val,m in l]
# 7717 vname = 'do_not_use_this_name_in_a_polynomial_coefficient'
# 7718 var = SR(vname)
# TypeError: 'sage.symbolic.expression.Expression' object is not iterable
```

Is this expected behavior? When I try `f.univariate_polynomial().roots(ring=ring, multiplicities=False)`

(setting `multiplicities=True`

is what `I.variety()`

does and leads to the error) I get `[1/2*x0,2]`

for roots, which makes little sense (it seems like the roots should be `[0,-x0]`

).

I would be happy to convert the `MPolynomial_polydict`

(or the `sage.rings.polynomial.polynomial_ring.PolynomialRing_field_with_category.element_class`

generated by `univariate_polynomial()`

) into a symbolic expression that I could just set it equal to zero and use `solve`

on it, but I don't know how to do this type conversion.

copying my comments to tmoteil here:

A less trivial problem looks like this:

When I do this, I get an empty list as the result.

I suppose the problem may be that there are no solutions that are valid over all possible values of the constants? Here's an example of there existing a solution if we substitute 1.0 (or 0.5) for all the constants:

The solution is {1/3,1/3,1/3}

Secretly I know that the value for every constant will be a positive number. And, while there may be lots of possible (x9,x10,x11) triples that satisfy the equations, the only ones that are "physically meaningful" are those with all three values between 0 and 1 (inclusive), possibly irrational. The fact that I cannot find any physically meaningful solutions after taking the Groebner basis and doing the rootfinding myself is what alerted me to a problem.

I'm not sure if these physically-motivated assumptions of mine are at the root of my confusion.