diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4c80a54..d3b0098 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,6 +11,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Custom `Monomial` and `Polynomial` classes, to better handle polynomials with coefficients in any of $\mathbb{C}, \mathbb{F}_p, \mathbb{Q}_p$.
 - DOI.
+- `Ideal.point_on_variety` now checks whether the provided `directions` are consistent with the given variety.
+- Added support for rational field, `Field('rational', 0, 0)`.
+- Tentatively introducing shorter aliases for field names, such as `C`, `Qp` and `Fp`.
 
 ### Changed
 
@@ -21,6 +24,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Improved logic for call to `Singular` with long commands requiring the writing of a file. As a consequence, parallelised calls should now have much better stability.
 
+### Deprecated
+
+- `Field.random_element` is now deprecated. Use `Field.random` instead.
+
 ## [0.2.1] - 2024-05-04
 
 ### Added
diff --git a/syngular/__init__.py b/syngular/__init__.py
index 9806f9a..6a2ab59 100644
--- a/syngular/__init__.py
+++ b/syngular/__init__.py
@@ -8,4 +8,4 @@
 from .qring import QuotientRing, QRing  # noqa
 from .tools import SingularException  # noqa
 from .field import Field  # noqa
-from .polynomial import Polynomial  # noqa
+from .polynomial import Polynomial, Monomial  # noqa
diff --git a/syngular/field.py b/syngular/field.py
index a261bc1..76a1b47 100644
--- a/syngular/field.py
+++ b/syngular/field.py
@@ -36,12 +36,14 @@ def __repr__(self):
 
     def __call__(self, other):
         """Cast to field."""
-        if self.name == "mpc":
+        if self.name in ["mpc", "C"]:
             return mpmath.mpc(mpmath.mpmathify(other))
-        elif self.name == "padic":
+        elif self.name in ["padic", "Qp"]:
             return PAdic(other, self.characteristic, self.digits)
-        elif self.name == "finite field":
+        elif self.name in ["finite field", "Fp"]:
             return ModP(other, self.characteristic)
+        elif self.name in ["rational", "Q"]:
+            return Q(other)
         else:
             raise NotImplementedError
 
@@ -51,7 +53,7 @@ def name(self):
 
     @name.setter
     def name(self, value):
-        if value not in ['mpc', 'gaussian rational', 'finite field', 'padic']:
+        if value not in ['rational', 'Q', 'mpc', 'C', 'gaussian rational', 'finite field', 'Fp', 'padic', 'Qp']:
             raise Exception("Field must be one of 'mpc', 'gaussian rational', 'finite field', 'padic'.")
         else:
             self._name = value
@@ -69,7 +71,7 @@ def characteristic(self, value):
 
     @property
     def digits(self):
-        if self.name == 'mpc':
+        if self.name in ['mpc', 'C']:
             return mpmath.mp.dps
         else:
             return self._digits
@@ -78,7 +80,7 @@ def digits(self):
     def digits(self, value):
         if value < 0:
             raise Exception("Digits must be positive.")
-        elif self.name == 'mpc':
+        elif self.name in ['mpc', 'C']:
             mpmath.mp.dps = value
         else:
             self._digits = value
@@ -88,7 +90,7 @@ def __contains__(self, other):
 
     @property
     def is_algebraically_closed(self):
-        if self.name == 'mpc':
+        if self.name in ['mpc', 'C']:
             return True
         else:
             return False
@@ -97,30 +99,30 @@ def is_algebraically_closed(self):
     def tollerance(self):
         if self.name in ['gaussian rational', 'finite field']:
             return 0
-        elif self.name == 'mpc':
+        elif self.name in ['mpc', 'C']:
             return mpmath.mpf('10e-{}'.format(int(min([0.95 * mpmath.mp.dps, mpmath.mp.dps - 4]))))
-        elif self.name == 'padic':
+        elif self.name in ['padic', 'Qp']:
             return PAdic(0, self.characteristic, 0, self.digits)
 
     @property
     def singular_notation(self):
-        if self.name == 'mpc':
+        if self.name in ['mpc', 'C']:
             return '(complex,{},I)'.format(self.digits - 5)
-        elif self.name in ['finite field', 'padic']:
+        elif self.name in ['finite field', 'Fp', 'padic', 'Qp']:
             return str(self.characteristic)
         else:
             return None
 
     def sqrt(self, val):
-        if self.name == "finite field":
+        if self.name in ["finite field", 'Fp']:
             if not isinstance(val, ModP):
                 val = self(val)
             return finite_field_sqrt(val)
-        elif self.name == "padic":
+        elif self.name in ["padic", 'Qp']:
             if not isinstance(val, PAdic):
                 val = self(val)
             return padic_sqrt(val)
-        elif self.name == "mpc":
+        elif self.name in ["mpc", 'C']:
             return mpmath.sqrt(val)
         else:
             raise Exception(f"Field not understood: {self.field.name}")
@@ -131,15 +133,15 @@ def j(self):
 
     i = I = j
 
-    def random_element(self, shape=(1, )):
+    def random(self, shape=(1, )):
         if shape == (1, ):
-            if self.name == "padic":
+            if self.name in ["padic", 'Qp']:
                 p, k = self.characteristic, self.digits
                 return PAdic(random.randrange(0, p ** k - 1), p, k)
-            elif self.name == "finite field":
+            elif self.name in ["finite field", 'Fp']:
                 p = self.characteristic
                 return ModP(random.randrange(0, p), p)
-            elif self.name == "mpc":
+            elif self.name in ["mpc", 'C']:
                 return mpmath.mpc(str(Q(random.randrange(-100, 101), random.randrange(1, 201))),
                                   str(Q(random.randrange(-100, 101), random.randrange(1, 201))))
             else:
@@ -147,6 +149,16 @@ def random_element(self, shape=(1, )):
         else:
             raise NotImplementedError
 
+    def random_element(self, *args, **kwargs):
+        import warnings
+        warnings.warn(
+            "random_element is deprecated and will be removed in a future version. "
+            "Use random instead.",
+            DeprecationWarning,
+            stacklevel=2
+        )
+        return self.random(*args, **kwargs)
+
     def epsilon(self, shape=(1, ), ):
         if shape == (1, ):
             if not hasattr(self, "_ε"):
@@ -165,7 +177,7 @@ def epsilon(self, shape=(1, ), ):
     @property
     def ε(self):
         if not hasattr(self, "_ε"):
-            self._ε = self.epsilon
+            self._ε = self.epsilon()
         return self._ε
 
     @ε.setter
diff --git a/syngular/ideal.py b/syngular/ideal.py
index 2704c94..74cf966 100644
--- a/syngular/ideal.py
+++ b/syngular/ideal.py
@@ -9,6 +9,7 @@
 from .qring import QuotientRing
 from .ideal_algorithms import Ideal_Algorithms
 from .variety import Variety_of_Ideal
+from .polynomial import Polynomial
 
 
 class Ideal(Ideal_Algorithms, Variety_of_Ideal, object):
@@ -192,8 +193,10 @@ def __contains__(self, other):
         """Implements ideal membership."""
         if isinstance(other, Ideal):
             return self.__ideal_contains__(other)
-        elif isinstance(other, sympy.Expr) or isinstance(other, str):
+        elif isinstance(other, (str, sympy.Expr, Polynomial)):
             return self.__poly_contains__(other)
+        else:
+            raise NotImplementedError
 
     def __ideal_contains__(self, other):
         qIdeal = self / other
@@ -203,7 +206,7 @@ def __ideal_contains__(self, other):
             return False
 
     def __poly_contains__(self, other):
-        if isinstance(other, sympy.Expr):
+        if not isinstance(other, str):
             other = str(other)
         assert isinstance(other, str)
         singular_commands = [f"ring r = {self.ring};",
diff --git a/syngular/polynomial.py b/syngular/polynomial.py
index a6608ea..711f518 100644
--- a/syngular/polynomial.py
+++ b/syngular/polynomial.py
@@ -2,6 +2,7 @@
 import numpy
 import operator
 import re
+import sympy
 
 from collections import defaultdict
 from multiset import FrozenMultiset
@@ -48,6 +49,9 @@ def __repr__(self):
     def __str__(self):
         return "*".join([f"{key}^{val}" for key, val in self.items()])
 
+    def tolist(self):
+        return [entry for key, val in self.items() for entry in [key, ] * val]
+
     def subs(self, values_dict):
         return functools.reduce(operator.mul, [values_dict[key] ** val for key, val in self.items()], 1)
 
@@ -77,12 +81,19 @@ class Polynomial(object):
     #     return self
 
     def __init__(self, coeffs_and_monomials, field):
-        if isinstance(coeffs_and_monomials, str):
-            coeffs_and_monomials = self.__rstr__(coeffs_and_monomials, field)
-        elif coeffs_and_monomials in field:
-            coeffs_and_monomials = [(coeffs_and_monomials, Monomial(''))]
+        if isinstance(coeffs_and_monomials, str) or isinstance(coeffs_and_monomials, sympy.Basic):
+            coeffs_and_monomials = self.__rstr__(str(coeffs_and_monomials), field)
         elif isinstance(coeffs_and_monomials, Polynomial):
             coeffs_and_monomials = coeffs_and_monomials.coeffs_and_monomials
+        elif isinstance(coeffs_and_monomials, list):
+            assert all([isinstance(entry, tuple) for entry in coeffs_and_monomials])
+            assert all([len(entry) == 2 for entry in coeffs_and_monomials])
+            assert all([entry[0] in field for entry in coeffs_and_monomials if entry[0] != 0])
+            assert all([isinstance(entry[1], Monomial) for entry in coeffs_and_monomials])
+        elif coeffs_and_monomials in field:
+            coeffs_and_monomials = [(coeffs_and_monomials, Monomial(''))]
+        else:
+            raise NotImplementedError(f"Received {coeffs_and_monomials} \n of type {type(coeffs_and_monomials)}")
         self.coeffs_and_monomials = coeffs_and_monomials
         self._field = field
 
@@ -91,7 +102,7 @@ def __str__(self, for_repr=False):
                            if hasattr(self.field, "name") and self.field.name in ["padic", "finite field", ] else
                            f"{str(coeff) if not for_repr else repr(coeff)}") +
                           (f"*{monomial}" if str(monomial) != "" else "")
-                          for coeff, monomial in self.coeffs_and_monomials).replace("+-", "-")
+                          for coeff, monomial in self.coeffs_and_monomials).replace("+-", "-").replace("+ -", "- ")
 
     def __repr__(self):
         return f"Polynomial(\"{self.__str__(for_repr=True)}\", {self.field})"
@@ -107,6 +118,8 @@ def coeffs_and_monomials(self):
     def coeffs_and_monomials(self, temp_coeffs_and_monomials):
         self._coeffs_and_monomials = [(coeff, monomial) for coeff, monomial in temp_coeffs_and_monomials if coeff != 0]
         self.reduce()
+        if self._coeffs_and_monomials == []:  # ensure there is at least an entry
+            self._coeffs_and_monomials = [(0, Monomial(''))]
 
     @property
     def field(self):
@@ -128,6 +141,10 @@ def coeffs(self):
     def coeffs(self, temp_coeffs):
         self.coeffs_and_monomials = [(temp_coeff, monomial) for temp_coeff, (coeff, monomial) in zip(temp_coeffs, self.coeffs_and_monomials)]
 
+    @property
+    def monomials(self):
+        return [monomial for _, monomial in self.coeffs_and_monomials]
+
     def rationalise(self):
         from pyadic.finite_field import vec_chained_FF_rationalize
         rat_coeffs = vec_chained_FF_rationalize([numpy.array(self.coeffs).astype(int), ], [self.field.characteristic, ]).tolist()
@@ -138,7 +155,7 @@ def rationalise(self):
 
     def reduce(self):
         """Merges equal monomials"""
-        unequal_monoms = set([monom for coeff, monom in self.coeffs_and_monomials])
+        unequal_monoms = set([monom for _, monom in self.coeffs_and_monomials])
         if len(unequal_monoms) == len(self.coeffs_and_monomials):
             return  # already reduced
         new_coeffs_and_monomials = []
@@ -149,10 +166,12 @@ def reduce(self):
 
     @staticmethod
     def __rstr__(polynomial, field):
+        polynomial = polynomial.replace(" ", "").replace("+-", "-")
+        polynomial = re.sub(r"(\+|\-)I\*{0,1}([\d\.]+)", r"\1\2j", polynomial)  # format complex nbrs
         parentheses = [(("(", ), (")", )), (("{", ), ("}", )), (("[", "⟨", "<", ), ("]", "⟩", ">"))]
         lopen_parentheses = [parenthesis[0] for parenthesis in parentheses]
         lclos_parentheses = [parenthesis[1] for parenthesis in parentheses]
-        parentheses_balance = [0 for parenthesis in parentheses]
+        parentheses_balance = [0 for _ in parentheses]
         next_match = ""
         coeffs_and_monomials_strings = []
         for char in polynomial:
@@ -201,12 +220,15 @@ def subs(self, base_point, field=None):
         new_poly.reduce()
         return new_poly
 
+    def __call__(self, *args, **kwargs):
+        return self.subs(*args, **kwargs)
+
     def __len__(self):
         return len(self.coeffs_and_monomials)
 
     def __eq__(self, other):
         if other == 0:
-            if self.coeffs_and_monomials == []:
+            if len(self) == 1 and self.coeffs[0] == 0 and self.monomials[0] == Monomial(''):
                 return True
             else:
                 return False
@@ -229,6 +251,9 @@ def __add__(self, other):
         else:
             raise NotImplementedError(f"Operation: __add__; self: {self}; self class {self.__class__}; other: {other}; other class {other.__class__}.")
 
+    def __sub__(self, other):
+        return self + (-1 * other)
+
     def __mul__(self, other):
         if isinstance(other, Polynomial):
             new_coeffs_and_monomials = []
@@ -245,6 +270,9 @@ def __mul__(self, other):
     def __rmul__(self, other):
         return self * other
 
+    def __neg__(self):
+        return -1 * self
+
     def __pow__(self, n):
         assert (isinstance(n, int) or n.is_integer())
         if n < 0:
diff --git a/syngular/variety.py b/syngular/variety.py
index df742ce..c55cf63 100644
--- a/syngular/variety.py
+++ b/syngular/variety.py
@@ -50,7 +50,8 @@ def wrapper(self, field, base_point={}, directions=None, valuations=tuple(), ind
 class Variety_of_Ideal:
 
     @retry_to_find_root(max_tries=100)
-    def point_on_variety(self, field, base_point={}, directions=None, valuations=tuple(), indepSetNbr=None, seed=None, verbose=False):
+    def point_on_variety(self, field, base_point={}, directions=None, valuations=tuple(), indepSetNbr=None, seed=None,
+                         verbose=False, directions_analytic_check=False):
         """Generate a representative point on or close to the variety associated to this ideal.
         The point is 'valuations' away from the exact variety, in the directions specified by 'directions'.
         If 'directions' are not provided, pick the first n=codim simplest generators from 'self'.
@@ -64,17 +65,26 @@ def point_on_variety(self, field, base_point={}, directions=None, valuations=tup
 
         random.seed(seed)
 
+        # do not modify directions, in case re-try is triggered, better the input is identical
+        directions = deepcopy(directions)
         # handle directions, i.e. the generators of the sub-ideal of maximal codimension
         if directions is None or directions == []:
             directions = []
-            if verbose:
-                print("Directions not provided, obtaining them from ideal generators.")
-            for poly in sorted(self.generators, key=lambda x: len(x)):
-                if Ideal(self.ring, directions + [poly, ]).codim > Ideal(self.ring, directions).codim:
-                    directions += [poly, ]
-                if len(directions) == self.codim:
-                    break
-            assert Ideal(self.ring, directions).codim == self.codim
+            if field.name != "finite field":
+                if verbose:
+                    print("Directions not provided, obtaining them from ideal generators.")
+                for poly in sorted(self.generators, key=lambda x: len(x)):
+                    if Ideal(self.ring, directions + [poly, ]).codim > Ideal(self.ring, directions).codim:
+                        directions += [poly, ]
+                    if len(directions) == self.codim:
+                        break
+                assert Ideal(self.ring, directions).codim == self.codim
+        elif directions_analytic_check:
+            # make sure the provided directions make sense, given the ideal (they need to belong to it)
+            # if this is not triggered, then the check is perfomed numerically later. The analytic check can be expensive.
+            for direction in directions:
+                if Polynomial(direction, Field("rational", 0, 0)) not in self:
+                    raise Exception(f"Invalid direction, {direction} was not in {self}.")
         # make sure provided directions are expanded strings
         for i, direction in enumerate(directions):
             if isinstance(direction, sympy.core.Basic):
@@ -87,7 +97,7 @@ def point_on_variety(self, field, base_point={}, directions=None, valuations=tup
         if field.name == "padic":
             prime, iterations = field.characteristic, field.digits
             if valuations == tuple():
-                valuations = tuple(field.digits for poly in directions)
+                valuations = tuple(field.digits for _ in directions)
         elif field.name == "finite field":
             prime, iterations = field.characteristic, 1
         elif field.name == "mpc":
@@ -120,7 +130,7 @@ def point_on_variety(self, field, base_point={}, directions=None, valuations=tup
 
         # handle the base point, i.e. the values of the independent variables
         if base_point == {}:
-            base_point = {indepSymbol: field.random_element() for indepSymbol in indepSymbols}
+            base_point = {indepSymbol: field.random() for indepSymbol in indepSymbols}
         else:
             base_point = {str(key): field(val) for key, val in base_point.items()}
         base_point |= {depSymbol: Polynomial(depSymbol, field) for depSymbol in depSymbols}
@@ -144,7 +154,7 @@ def point_on_variety(self, field, base_point={}, directions=None, valuations=tup
                 raise Exception("The dimension of the semi-numerical ideal was not zero.")
 
             root_dicts = lex_groebner_solve(oSemiNumericalIdeal.groebner_basis, prime=prime)
-            check_solutions(oSemiNumericalIdeal.groebner_basis, root_dicts, prime=prime)
+            check_solutions(oSemiNumericalIdeal.groebner_basis, root_dicts, field)  # they may be stricter then wanted for mpc.
 
             try:
                 root_dict = root_dicts[0]
@@ -167,6 +177,15 @@ def point_on_variety(self, field, base_point={}, directions=None, valuations=tup
             update_point_dict(base_point, root_dict, field)
             # print("updated point:", base_point)
 
+            if iteration == 0 and not directions_analytic_check:
+                # instead of analytically checking the directions are consistent with the ideal
+                # we check they vanish numerically at the constructed point
+                for direction in directions:
+                    num_poly = Polynomial(direction, field).subs(base_point).subs({key: 0 for key in depSymbols}).coeffs[0]
+                    # print("val:", float(abs(num_poly)))
+                    if abs(num_poly) > abs(field.ε):
+                        raise Exception(f"Invalid direction, {direction} was not in {self}. Numerical membership check failed.")
+
             if iteration < iterations - 1:
                 if prime is not None:
                     valuations = [valuation - 1 for valuation in valuations]
@@ -336,18 +355,14 @@ def lex_groebner_solve(equations, prime=None):
     return root_dicts
 
 
-def check_solutions(equations, root_dicts, prime=None):
+def check_solutions(equations, root_dicts, field):
     """Checks that all solutions in root_dicts solve the equations."""
+    field = field if field.name not in ["padic", "Qp"] else Field("finite field", field.characteristic, 1)
     for root_dict in root_dicts:
-        check = []
-        for equation in equations:
-            check += [sympy.sympify(equation)]
-            for key, value in root_dict.items():
-                if prime is None:
-                    check[-1] = sympy.simplify(check[-1].subs(key, value))
-                else:
-                    check[-1] = check[-1].subs(key, value) % prime
-        if prime is None:
-            assert all([numpy.isclose(complex(entry), 0) for entry in check])
-        else:
-            assert all([entry == 0 for entry in check])
+        check = [Polynomial(equation, field).subs(root_dict) for equation in equations]
+        assert all([len(entry) == 1 for entry in check])  # check these are elements of the field
+        if not all([abs(entry.coeffs[0]) <= field.tollerance for entry in check]):
+            if field.characteristic == 0:
+                raise RootPrecisionError
+            else:
+                raise AssertionError
diff --git a/tests/test_poly.py b/tests/test_poly.py
new file mode 100644
index 0000000..def70cf
--- /dev/null
+++ b/tests/test_poly.py
@@ -0,0 +1,44 @@
+from syngular import Monomial, Polynomial, Field
+
+from fractions import Fraction
+
+from pyadic import ModP
+
+
+coeffs_and_monomials = [
+    (Fraction(-215, 1), Monomial("wb^1*z^9*zb^7*w^5*X^10")),
+    (Fraction(1156, 1), Monomial("wb^3*z^3*zb^2*w^11*X^4")),
+    (Fraction(-125, 2), Monomial("z^4*zb^6*w^8*X^6")),
+    (Fraction(-9857, 2), Monomial("wb^1*z^6*zb^6*w^5*X^8")),
+    (Fraction(-954, 1), Monomial("wb^3*z^9*zb^6*w^4*X^9")),
+    (Fraction(2, 1), Monomial("z^9*zb^8*X^9")),
+    (Fraction(-4485, 2), Monomial("wb^4*z^2*zb^5*w^9*X^5")),
+    (Fraction(216, 1), Monomial("wb^4*z^8*zb^6*w^2*X^11")),
+    (Fraction(-7, 2), Monomial("wb^5*z^9*zb^1*w^9*X^6")),
+    (Fraction(625, 2), Monomial("wb^4*z^2*zb^7*w^8*X^6"))
+]
+
+point = {
+    'z': ModP('1978390662 % 2147481317'),
+    'zb': ModP('117802826 % 2147481317'),
+    'w': ModP('962111889 % 2147481317'),
+    'wb': ModP('1912527506 % 2147481317'),
+    'X': ModP('747011570 % 2147481317')
+}
+
+QQ = Field("rational", 0, 0)
+Fp = Field('finite field', 2147481317, 1)
+
+
+def test_eq_zero():
+    poly = Polynomial(coeffs_and_monomials, QQ)
+    assert poly - poly == 0
+
+
+def test_rstr():
+    poly = Polynomial(coeffs_and_monomials, QQ)
+    poly_reloaded = Polynomial(str(poly), QQ)
+    assert len(poly_reloaded) == len(poly)
+    assert Polynomial(poly_reloaded, QQ) == poly
+    assert Polynomial(poly_reloaded, QQ) - poly == 0
+    assert poly.subs(point, Fp) == poly_reloaded.subs(point, Fp)
diff --git a/tests/test_variety.py b/tests/test_variety.py
index fb4d3c7..cafafb2 100644
--- a/tests/test_variety.py
+++ b/tests/test_variety.py
@@ -1,8 +1,9 @@
 import numpy
+import sympy
 
 from copy import copy
 
-from syngular import Field, Ring, Ideal
+from syngular import Field, Ring, QRing, Ideal, Polynomial
 
 
 def test_finite_field_variety_point():
@@ -51,3 +52,24 @@ def test_complex_variety_point_with_base_point_and_unequal_valuations():
     assert numpy.isclose(abs(complex(eval('zb^2*wb*X^2-zb*w*wb*X-zb^2*X^2-zb*wb*X^2+zb*wb*X+zb*X^2-w*wb+wb'.replace("^", "**"), copy(point_dict)))), 10 ** -30)
     assert numpy.isclose(abs(complex(eval('z*zb*X-z*w+z+w'.replace("^", "**"), copy(point_dict)))), 10 ** -60)
     assert [point_dict[key] == base_point[key] for key in base_point.keys()].count(True) == 3
+
+
+def test_padic_variety_point_in_qring():
+    """This is ("Δ_14|23|56", "⟨1|2+3|4]" ) from lips at six point."""
+    ring = Ring(0, ('a1', 'b1', 'c1', 'd1', 'a2', 'b2', 'c2', 'd2', 'a3', 'b3', 'c3', 'd3', 'a4', 'b4', 'c4', 'd4', 'a5', 'b5', 'c5', 'd5', 'a6', 'b6', 'c6', 'd6'), 'dp')
+    ideal = Ideal(ring, ['b1*d1 + b2*d2 + b3*d3 + b4*d4 + b5*d5 + b6*d6', '-a1*d1 - a2*d2 - a3*d3 - a4*d4 - a5*d5 - a6*d6',
+                         '-b1*c1 - b2*c2 - b3*c3 - b4*c4 - b5*c5 - b6*c6', 'a1*c1 + a2*c2 + a3*c3 + a4*c4 + a5*c5 + a6*c6'])
+    qring = QRing(ring, ideal)
+    directions = [
+        '-((a1*c1 + a4*c4)*(b1*d1 + b4*d4) - (-a1*d1 - a4*d4)*(-b1*c1 - b4*c4))*((a2*c2 + a3*c3)*(b2*d2 + b3*d3) - (a2*d2 + a3*d3)*(b2*c2 + b3*c3)) + ((a1*c1 + a4*c4)*(b2*d2 + b3*d3)/2 + (-a1*d1 - a4*d4)*(b2*c2 + b3*c3)/2 + (a2*c2 + a3*c3)*(b1*d1 + b4*d4)/2 + (a2*d2 + a3*d3)*(-b1*c1 - b4*c4)/2)**2',  # noqa
+        '-c4*(-a1*(b2*d2 + b3*d3) + b1*(a2*d2 + a3*d3)) + d4*(-a1*(b2*c2 + b3*c3) + b1*(a2*c2 + a3*c3))'
+    ]
+    directions = [sympy.sympify(direction) for direction in directions]
+    J = Ideal(qring, directions)
+    Qp = Field("padic", 2 ** 31 - 1, 11)
+    point = J.point_on_variety(Qp, directions=directions, valuations=(2, 2))
+    assert all([Polynomial(direction.expand(), Field("rational", 0, 0))(point, Qp).coeffs[0].n == 2 for direction in directions])
+    # This is (s_123-s_234), which lives in the radical of J but not in J.
+    member_of_radical = '(a1*c1 + a2*c2 + a3*c3)*(b1*d1 + b2*d2 + b3*d3) - (-a1*d1 - a2*d2 - a3*d3)*(-b1*c1 - b2*c2 - b3*c3) - (a2*c2 + a3*c3 + a4*c4)*(b2*d2 + b3*d3 + b4*d4) + (-a2*d2 - a3*d3 - a4*d4)*(-b2*c2 - b3*c3 - b4*c4)'  # noqa
+    member_of_radical = sympy.sympify(member_of_radical).expand()
+    assert Polynomial(member_of_radical.expand(), Field("rational", 0, 0))(point, Qp).coeffs[0].n == 1