diff --git a/mathics/builtin/__init__.py b/mathics/builtin/__init__.py
index d8cbc3ed4..f81819c23 100755
--- a/mathics/builtin/__init__.py
+++ b/mathics/builtin/__init__.py
@@ -50,6 +50,8 @@
 
 
 def add_builtins(new_builtins):
+    from mathics.core.convert.sympy import mathics_to_sympy, sympy_to_mathics
+
     for var_name, builtin in new_builtins:
         name = builtin.get_name()
         if hasattr(builtin, "python_equivalent"):
@@ -255,8 +257,6 @@ def name_is_builtin_symbol(module, name: str) -> Optional[type]:
                 _builtins_list.append((instance.get_name(), instance))
                 builtins_by_module[module.__name__].append(instance)
 
-mathics_to_sympy = {}  # here we have: name -> sympy object
-sympy_to_mathics = {}
 
 new_builtins = _builtins_list
 
diff --git a/mathics/builtin/arithmetic.py b/mathics/builtin/arithmetic.py
index 5a85bbe7e..d54f8f7bd 100644
--- a/mathics/builtin/arithmetic.py
+++ b/mathics/builtin/arithmetic.py
@@ -18,6 +18,7 @@
     IterationFunction,
     Predefined,
     SympyFunction,
+    SympyObject,
     Test,
 )
 from mathics.builtin.inference import evaluate_predicate, get_assumptions_list
@@ -788,7 +789,7 @@ def eval_Element_alternatives(
         return Element(Expression(elems.head, *unknown), domain)
 
 
-class I_(Predefined):
+class I_(Predefined, SympyObject):
     """
     <url>:Imaginary unit:https://en.wikipedia.org/wiki/Imaginary_unit</url> \
     (<url>:WMA:https://reference.wolfram.com/language/ref/I.html</url>)
@@ -805,9 +806,17 @@ class I_(Predefined):
     """
 
     name = "I"
+    sympy_name = "I"
+    sympy_obj = sympy.I
     summary_text = "imaginary unit"
     python_equivalent = 1j
 
+    def is_constant(self) -> bool:
+        return True
+
+    def to_sympy(self, symb, **kwargs):
+        return self.sympy_obj
+
     def evaluate(self, evaluation: Evaluation):
         return Complex(Integer0, Integer1)
 
diff --git a/mathics/builtin/base.py b/mathics/builtin/base.py
index 8eaccab79..ac3f1e1b7 100644
--- a/mathics/builtin/base.py
+++ b/mathics/builtin/base.py
@@ -19,10 +19,10 @@
     String,
 )
 from mathics.core.attributes import A_HOLD_ALL, A_NO_ATTRIBUTES, A_PROTECTED
-from mathics.core.convert.expression import to_expression, to_numeric_sympy_args
+from mathics.core.convert.expression import to_expression
 from mathics.core.convert.op import ascii_operator_to_symbol
 from mathics.core.convert.python import from_bool
-from mathics.core.convert.sympy import from_sympy
+from mathics.core.convert.sympy import from_sympy, to_numeric_sympy_args
 from mathics.core.definitions import Definition
 from mathics.core.evaluation import Evaluation
 from mathics.core.exceptions import MessageException
diff --git a/mathics/builtin/specialfns/elliptic.py b/mathics/builtin/specialfns/elliptic.py
index 0b9009462..928b4bd26 100644
--- a/mathics/builtin/specialfns/elliptic.py
+++ b/mathics/builtin/specialfns/elliptic.py
@@ -15,8 +15,7 @@
 from mathics.builtin.base import SympyFunction
 from mathics.core.atoms import Integer
 from mathics.core.attributes import A_LISTABLE, A_NUMERIC_FUNCTION, A_PROTECTED
-from mathics.core.convert.expression import to_numeric_sympy_args
-from mathics.core.convert.sympy import from_sympy
+from mathics.core.convert.sympy import from_sympy, to_numeric_sympy_args
 from mathics.eval.numerify import numerify
 
 
diff --git a/mathics/core/convert/expression.py b/mathics/core/convert/expression.py
index fc96fd537..d4360a3f2 100644
--- a/mathics/core/convert/expression.py
+++ b/mathics/core/convert/expression.py
@@ -99,23 +99,6 @@ def to_numeric_args(mathics_args: Type[BaseElement], evaluation) -> list:
     )
 
 
-def to_numeric_sympy_args(mathics_args: Type[BaseElement], evaluation) -> list:
-    """
-    Convert Mathics arguments, such as the arguments in an evaluation
-    method a Python list that is sutiable for feeding as arguments
-    into SymPy.
-
-    We make use of fast conversions for literals.
-    """
-    if mathics_args.is_literal:
-        sympy_args = [mathics_args.value]
-    else:
-        args = numerify(mathics_args, evaluation).get_sequence()
-        sympy_args = [a.to_sympy() for a in args]
-
-    return sympy_args
-
-
 expression_constructor_map = {
     SymbolList: lambda head, *args, **kwargs: ListExpression(*args, **kwargs)
 }
diff --git a/mathics/core/convert/sympy.py b/mathics/core/convert/sympy.py
index 980da1d69..843d679f3 100644
--- a/mathics/core/convert/sympy.py
+++ b/mathics/core/convert/sympy.py
@@ -5,17 +5,45 @@
 Conversion to SymPy is handled directly in BaseElement descendants.
 """
 
-from typing import Optional
+from typing import Optional, Type, Union
 
 import sympy
+from sympy import Symbol as Sympy_Symbol, false as SympyFalse, true as SympyTrue
+
+# Import the singleton class
+from sympy.core.numbers import S
 
 BasicSympy = sympy.Expr
 
 
+from mathics.core.atoms import (
+    MATHICS3_COMPLEX_I,
+    Complex,
+    Integer,
+    Integer0,
+    Integer1,
+    IntegerM1,
+    MachineReal,
+    Rational,
+    RationalOneHalf,
+    Real,
+    String,
+)
+from mathics.core.convert.expression import to_expression, to_mathics_list
 from mathics.core.convert.matrix import matrix_data
+from mathics.core.element import BaseElement
+from mathics.core.expression import Expression
+from mathics.core.expression_predefined import (
+    MATHICS3_COMPLEX_INFINITY,
+    MATHICS3_INFINITY,
+    MATHICS3_NEG_INFINITY,
+)
+from mathics.core.list import ListExpression
+from mathics.core.number import FP_MANTISA_BINARY_DIGITS
 from mathics.core.symbols import (
     Symbol,
     SymbolFalse,
+    SymbolNull,
     SymbolPlus,
     SymbolPower,
     SymbolTimes,
@@ -25,16 +53,20 @@
 )
 from mathics.core.systemsymbols import (
     SymbolC,
+    SymbolCatalan,
+    SymbolE,
     SymbolEqual,
+    SymbolEulerGamma,
     SymbolFunction,
+    SymbolGoldenRatio,
     SymbolGreater,
     SymbolGreaterEqual,
     SymbolIndeterminate,
-    SymbolInfinity,
     SymbolLess,
     SymbolLessEqual,
     SymbolMatrixPower,
     SymbolO,
+    SymbolPi,
     SymbolPiecewise,
     SymbolSlot,
     SymbolUnequal,
@@ -45,6 +77,36 @@
 SymbolRootSum = Symbol("RootSum")
 
 
+mathics_to_sympy = {}  # here we have: name -> sympy object
+sympy_to_mathics = {}
+
+
+sympy_singleton_to_mathics = {
+    None: SymbolNull,
+    S.Catalan: SymbolCatalan,
+    S.ComplexInfinity: MATHICS3_COMPLEX_INFINITY,
+    S.EulerGamma: SymbolEulerGamma,
+    S.Exp1: SymbolE,
+    S.GoldenRatio: SymbolGoldenRatio,
+    S.Half: RationalOneHalf,
+    S.ImaginaryUnit: MATHICS3_COMPLEX_I,
+    S.Infinity: MATHICS3_INFINITY,
+    S.NaN: SymbolIndeterminate,
+    S.NegativeInfinity: MATHICS3_NEG_INFINITY,
+    S.NegativeOne: IntegerM1,
+    S.One: Integer1,
+    S.Pi: SymbolPi,
+    S.Zero: Integer0,
+    SympyFalse: SymbolFalse,
+    SympyTrue: SymbolTrue,
+}
+
+
+mathics_to_sympy_singleton = {
+    key: val for val, key in sympy_singleton_to_mathics.items()
+}
+
+
 def is_Cn_expr(name) -> bool:
     if name.startswith(sympy_symbol_prefix) or name.startswith(sympy_slot_prefix):
         return False
@@ -75,7 +137,6 @@ class SympyExpression(BasicSympy):
     def __new__(cls, *exprs):
         # sympy simplify may also recreate the object if simplification occurred
         # in the elements
-        from mathics.core.expression import Expression
 
         if all(isinstance(expr, BasicSympy) for expr in exprs):
             # called with SymPy arguments
@@ -153,23 +214,125 @@ def eval(cls, n):
                 pass
 
 
-def from_sympy(expr):
-    from mathics.builtin import sympy_to_mathics
-    from mathics.core.atoms import (
-        Complex,
-        Integer,
-        Integer0,
-        Integer1,
-        MachineReal,
-        Rational,
-        Real,
-        String,
-    )
-    from mathics.core.convert.expression import to_expression, to_mathics_list
-    from mathics.core.expression import Expression
-    from mathics.core.list import ListExpression
-    from mathics.core.number import FP_MANTISA_BINARY_DIGITS
-    from mathics.core.symbols import Symbol, SymbolNull
+def expression_to_sympy(expr: Expression, **kwargs):
+    """
+    Convert `expr` to its sympy form.
+    """
+
+    if "convert_all_global_functions" in kwargs:
+        if len(expr.elements) > 0 and kwargs["convert_all_global_functions"]:
+            if expr.get_head_name().startswith("Global`"):
+                return expr._as_sympy_function(**kwargs)
+
+    if "converted_functions" in kwargs:
+        functions = kwargs["converted_functions"]
+        if len(expr._elements) > 0 and expr.get_head_name() in functions:
+            sym_args = [element.to_sympy() for element in expr._elements]
+            if None in sym_args:
+                return None
+            func = sympy.Function(str(sympy_symbol_prefix + expr.get_head_name()))(
+                *sym_args
+            )
+            return func
+
+    lookup_name = expr.get_lookup_name()
+    builtin = mathics_to_sympy.get(lookup_name)
+    if builtin is not None:
+        sympy_expr = builtin.to_sympy(expr, **kwargs)
+        if sympy_expr is not None:
+            return sympy_expr
+    return SympyExpression(expr)
+
+
+def symbol_to_sympy(symbol: Symbol, **kwargs) -> Sympy_Symbol:
+    """
+    Convert `symbol` to its sympy form.
+    """
+
+    result = mathics_to_sympy_singleton.get(symbol, None)
+    if result is not None:
+        return result
+
+    if symbol.sympy_dummy is not None:
+        return symbol.sympy_dummy
+
+    builtin = mathics_to_sympy.get(symbol.name)
+    if builtin is None or not builtin.sympy_name or not builtin.is_constant():  # nopep8
+        return Sympy_Symbol(sympy_symbol_prefix + symbol.name)
+    return builtin.to_sympy(symbol, **kwargs)
+
+
+def to_numeric_sympy_args(mathics_args: Type[BaseElement], evaluation) -> list:
+    """
+    Convert Mathics arguments, such as the arguments in an evaluation
+    method a Python list that is sutiable for feeding as arguments
+    into SymPy.
+
+    We make use of fast conversions for literals.
+    """
+    from mathics.eval.numerify import numerify
+
+    if mathics_args.is_literal:
+        sympy_args = [mathics_args.value]
+    else:
+        args = numerify(mathics_args, evaluation).get_sequence()
+        sympy_args = [a.to_sympy() for a in args]
+
+    return sympy_args
+
+
+def from_sympy_matrix(
+    expr: Union[sympy.Matrix, sympy.ImmutableMatrix]
+) -> ListExpression:
+    """
+    Convert `expr` of the type sympy.Matrix or sympy.ImmutableMatrix to
+    a Mathics list.
+    """
+    if len(expr.shape) == 2 and (expr.shape[1] == 1):
+        # This is a vector (only one column)
+        # Transpose and select first row to get result equivalent to Mathematica
+        return to_mathics_list(*expr.T.tolist()[0], elements_conversion_fn=from_sympy)
+    else:
+        return to_mathics_list(*expr.tolist(), elements_conversion_fn=from_sympy)
+
+
+"""
+sympy_conversion_by_type = {
+    complex: lambda expr: Complex(Real(expr.real), Real(expr.imag)),
+    int: lambda x: Integer(x),
+    float: lambda x: Real(x),
+    tuple: lambda expr: to_mathics_list(*expr, elements_conversion_fn=from_sympy),
+    list: lambda expr: to_mathics_list(*expr, elements_conversion_fn=from_sympy),
+    str: lambda x: String(x),
+    sympy.Matrix :from_sympy_matrix,
+    sympy.ImmutableMatrix :from_sympy_matrix,
+    sympy.MatPow: lambda expr: Expression(
+            SymbolMatrixPower, from_sympy(expr.base), from_sympy(expr.exp)
+        ),
+    SympyExpression: lambda expr: expr.expr,
+    SympyPrime: lambda expr: Expression(SymbolPrime, from_sympy(expr.args[0])),
+    sympy.RootSum: lambda expr: Expression(SymbolRootSum, from_sympy(expr.poly), from_sympy(expr.fun)),
+    sympy.Tuple: lambda expr: to_mathics_list(*expr, elements_conversion_fn=from_sympy),
+}
+
+"""
+
+# def new_from_sympy(expr)->BaseElement:
+#    """
+#    converts a SymPy object to a Mathics element.
+#    """
+#    try:
+#        return sympy_singleton_to_mathics[expr]
+#    except (KeyError, TypeError):
+#        pass
+#
+#    return sympy_conversion_by_type.get(type(expr), old_from_sympy)(expr)
+
+
+def old_from_sympy(expr) -> BaseElement:
+    """
+    converts a SymPy object to a Mathics element.
+    """
 
     if isinstance(expr, (tuple, list)):
         return to_mathics_list(*expr, elements_conversion_fn=from_sympy)
@@ -184,14 +347,7 @@ def from_sympy(expr):
     if expr is None:
         return SymbolNull
     if isinstance(expr, sympy.Matrix) or isinstance(expr, sympy.ImmutableMatrix):
-        if len(expr.shape) == 2 and (expr.shape[1] == 1):
-            # This is a vector (only one column)
-            # Transpose and select first row to get result equivalent to Mathematica
-            return to_mathics_list(
-                *expr.T.tolist()[0], elements_conversion_fn=from_sympy
-            )
-        else:
-            return to_mathics_list(*expr.tolist(), elements_conversion_fn=from_sympy)
+        return from_sympy_matrix(expr)
     if isinstance(expr, sympy.MatPow):
         return Expression(
             SymbolMatrixPower, from_sympy(expr.base), from_sympy(expr.exp)
@@ -218,23 +374,23 @@ def from_sympy(expr):
             if builtin is not None:
                 name = builtin.get_name()
             return Symbol(name)
-        elif isinstance(
-            expr, (sympy.core.numbers.Infinity, sympy.core.numbers.ComplexInfinity)
-        ):
-            return Symbol(expr.__class__.__name__)
+        elif isinstance(expr, sympy.core.numbers.Infinity):
+            return MATHICS3_INFINITY
+        elif isinstance(expr, sympy.core.numbers.ComplexInfinity):
+            return MATHICS3_COMPLEX_INFINITY
         elif isinstance(expr, sympy.core.numbers.NegativeInfinity):
-            return Expression(SymbolTimes, Integer(-1), SymbolInfinity)
+            return MATHICS3_NEG_INFINITY
         elif isinstance(expr, sympy.core.numbers.ImaginaryUnit):
-            return Complex(Integer0, Integer1)
+            return MATHICS3_COMPLEX_I
         elif isinstance(expr, sympy.Integer):
             return Integer(int(expr))
         elif isinstance(expr, sympy.Rational):
             numerator, denominator = map(int, expr.as_numer_denom())
             if denominator == 0:
                 if numerator > 0:
-                    return SymbolInfinity
+                    return MATHICS3_INFINITY
                 elif numerator < 0:
-                    return Expression(SymbolTimes, Integer(-1), SymbolInfinity)
+                    return MATHICS3_NEG_INFINITY
                 else:
                     assert numerator == 0
                     return SymbolIndeterminate
@@ -408,3 +564,6 @@ def from_sympy(expr):
                 expr, str(expr.__class__)
             )
         )
+
+
+from_sympy = old_from_sympy
diff --git a/mathics/core/expression.py b/mathics/core/expression.py
index fec9820b9..74b4c6b9e 100644
--- a/mathics/core/expression.py
+++ b/mathics/core/expression.py
@@ -26,7 +26,6 @@
     attribute_string_to_number,
 )
 from mathics.core.convert.python import from_python
-from mathics.core.convert.sympy import SympyExpression, sympy_symbol_prefix
 from mathics.core.element import ElementsProperties, EvalMixin, ensure_context
 from mathics.core.evaluation import Evaluation
 from mathics.core.interrupt import ReturnInterrupt
@@ -258,6 +257,8 @@ def __str__(self) -> str:
         )
 
     def _as_sympy_function(self, **kwargs) -> sympy.Function:
+        from mathics.core.convert.sympy import sympy_symbol_prefix
+
         sym_args = [element.to_sympy(**kwargs) for element in self._elements]
 
         if None in sym_args:
@@ -1461,32 +1462,9 @@ def to_python(self, *args, **kwargs):
         return self
 
     def to_sympy(self, **kwargs):
-        from mathics.builtin import mathics_to_sympy
-
-        if "convert_all_global_functions" in kwargs:
-            if len(self.elements) > 0 and kwargs["convert_all_global_functions"]:
-                if self.get_head_name().startswith("Global`"):
-                    return self._as_sympy_function(**kwargs)
-
-        if "converted_functions" in kwargs:
-            functions = kwargs["converted_functions"]
-            if len(self._elements) > 0 and self.get_head_name() in functions:
-                sym_args = [element.to_sympy() for element in self._elements]
-                if None in sym_args:
-                    return None
-                func = sympy.Function(str(sympy_symbol_prefix + self.get_head_name()))(
-                    *sym_args
-                )
-                return func
-
-        lookup_name = self.get_lookup_name()
-        builtin = mathics_to_sympy.get(lookup_name)
-        if builtin is not None:
-            sympy_expr = builtin.to_sympy(self, **kwargs)
-            if sympy_expr is not None:
-                return sympy_expr
+        from mathics.core.convert.sympy import expression_to_sympy
 
-        return SympyExpression(self)
+        return expression_to_sympy(self, **kwargs)
 
     def process_style_box(self, options):
         if self.has_form("StyleBox", 1, None):
diff --git a/mathics/core/symbols.py b/mathics/core/symbols.py
index cdfd6bf9d..d49a4e031 100644
--- a/mathics/core/symbols.py
+++ b/mathics/core/symbols.py
@@ -4,8 +4,6 @@
 import time
 from typing import Any, FrozenSet, List, Optional, Tuple
 
-import sympy
-
 from mathics.core.element import (
     BaseElement,
     EvalMixin,
@@ -624,19 +622,9 @@ def to_python(self, *args, python_form: bool = False, **kwargs):
         return self.name
 
     def to_sympy(self, **kwargs):
-        from mathics.builtin import mathics_to_sympy
-
-        if self.sympy_dummy is not None:
-            return self.sympy_dummy
-
-        builtin = mathics_to_sympy.get(self.name)
-        if (
-            builtin is None
-            or not builtin.sympy_name
-            or not builtin.is_constant()  # nopep8
-        ):
-            return sympy.Symbol(sympy_symbol_prefix + self.name)
-        return builtin.to_sympy(self, **kwargs)
+        from mathics.core.convert.sympy import symbol_to_sympy
+
+        return symbol_to_sympy(self, **kwargs)
 
 
 class SymbolConstant(Symbol):
diff --git a/mathics/doc/latex/Makefile b/mathics/doc/latex/Makefile
index a66dc37bf..9f4369ad3 100644
--- a/mathics/doc/latex/Makefile
+++ b/mathics/doc/latex/Makefile
@@ -10,7 +10,7 @@ BASH ?= /bin/bash
 DOCTEST_LATEX_DATA_PCL ?= $(HOME)/.local/var/mathics/doctest_latex_data.pcl
 
 # Variable indicating Mathics3 Modules you have available on your system, in latex2doc option format
-MATHICS3_MODULE_OPTION ?= --load-module pymathics.graph,pymathics.natlang
+MATHICS3_MODULE_OPTION ?= # --load-module pymathics.graph,pymathics.natlang
 
 #: Default target: Make everything
 all doc texdoc: mathics.pdf
diff --git a/mathics/docpipeline.py b/mathics/docpipeline.py
index 81c32f69d..0696f4e59 100644
--- a/mathics/docpipeline.py
+++ b/mathics/docpipeline.py
@@ -449,6 +449,12 @@ def save_doctest_data(output_data: Dict[tuple, dict]):
         should_be_readable=False, create_parent=True
     )
     print(f"Writing internal document data to {doctest_latex_data_path}")
+    i = 0
+    for key in output_data:
+        i = i + 1
+        print(key, output_data[key])
+        if i > 9:
+            break
     with open(doctest_latex_data_path, "wb") as output_file:
         pickle.dump(output_data, output_file, 4)
 
diff --git a/mathics/eval/nevaluator.py b/mathics/eval/nevaluator.py
index 17e4f8d39..79d872507 100644
--- a/mathics/eval/nevaluator.py
+++ b/mathics/eval/nevaluator.py
@@ -15,7 +15,6 @@
 
 from mathics.core.atoms import Number
 from mathics.core.attributes import A_N_HOLD_ALL, A_N_HOLD_FIRST, A_N_HOLD_REST
-from mathics.core.convert.sympy import from_sympy
 from mathics.core.element import BaseElement
 from mathics.core.evaluation import Evaluation
 from mathics.core.expression import Expression
@@ -50,6 +49,8 @@ def eval_NValues(
     stored in ``evaluation.definitions``.
     If ``prec`` can not be evaluated as a number, returns None, otherwise, returns an expression.
     """
+    from mathics.core.convert.sympy import from_sympy
+
     # The first step is to determine the precision goal
     try:
         # Here ``get_precision`` is called with ``show_messages``
diff --git a/test/core/convert/mpmath.py b/test/core/convert/test_mpmath.py
similarity index 100%
rename from test/core/convert/mpmath.py
rename to test/core/convert/test_mpmath.py
diff --git a/test/core/convert/test_sympy.py b/test/core/convert/test_sympy.py
new file mode 100644
index 000000000..319274366
--- /dev/null
+++ b/test/core/convert/test_sympy.py
@@ -0,0 +1,144 @@
+import test.helper
+
+import pytest
+from sympy import Float as SympyFloat
+
+from mathics.core.atoms import (
+    MATHICS3_COMPLEX_I,
+    Complex,
+    Integer,
+    Integer0,
+    Integer1,
+    Integer2,
+    Integer3,
+    IntegerM1,
+    MachineReal,
+    PrecisionReal,
+    Rational,
+    RationalOneHalf,
+    Real,
+    String,
+)
+from mathics.core.convert.sympy import from_sympy, sympy_singleton_to_mathics
+from mathics.core.expression import Expression
+from mathics.core.expression_predefined import (
+    MATHICS3_COMPLEX_INFINITY,
+    MATHICS3_INFINITY,
+)
+from mathics.core.symbols import (
+    Symbol,
+    SymbolFalse,
+    SymbolNull,
+    SymbolPlus,
+    SymbolPower,
+    SymbolTimes,
+    SymbolTrue,
+)
+from mathics.core.systemsymbols import (
+    SymbolAnd,
+    SymbolComplexInfinity,
+    SymbolDirectedInfinity,
+    SymbolE,
+    SymbolExp,
+    SymbolI,
+    SymbolIndeterminate,
+    SymbolInfinity,
+    SymbolPi,
+    SymbolSin,
+)
+
+Symbol_a = Symbol("Global`a")
+Symbol_b = Symbol("Global`b")
+Symbol_x = Symbol("Global`x")
+Symbol_y = Symbol("Global`y")
+Symbol_F = Symbol("Global`F")
+Symbol_G = Symbol("Global`G")
+
+
+@pytest.mark.parametrize(
+    ("expr",),
+    [
+        (Symbol_x,),
+        (Expression(Symbol_F, Symbol_x),),
+        (SymbolPi,),
+        (SymbolTrue,),
+        (SymbolFalse,),
+        (Integer1,),
+        (Integer(37),),
+        (Rational(1, 5),),
+        (Real(1.2),),
+        (Real(SympyFloat(1.2, 10)),),
+        (Complex(Real(2.0), Real(3.0)),),
+        (Expression(Symbol_F, Symbol_x, SymbolPi),),
+        (Expression(Symbol_G, Expression(Symbol_F, Symbol_x, SymbolPi)),),
+        (Expression(SymbolPlus, Integer3, Symbol_x, Symbol_y),),
+    ],
+)
+def test_from_to_sympy_invariant(expr):
+    """
+    Check if the conversion back and forward is consistent.
+    """
+    result_sympy = expr.to_sympy()
+    back_to_mathics = from_sympy(result_sympy)
+    print([expr, result_sympy, back_to_mathics])
+    assert expr.sameQ(back_to_mathics)
+
+
+@pytest.mark.parametrize(
+    ("expr", "result", "msg"),
+    [
+        (
+            Expression(SymbolExp, Expression(SymbolTimes, SymbolI, SymbolPi)),
+            IntegerM1,
+            None,
+        ),
+        (
+            Expression(
+                SymbolPower, SymbolE, Expression(SymbolTimes, SymbolI, SymbolPi)
+            ),
+            IntegerM1,
+            None,
+        ),
+        (Expression(SymbolSin, SymbolPi), Integer0, None),
+        (Expression(SymbolPlus, Integer1, Integer2), Integer3, None),
+        (String("Hola"), SymbolNull, None),
+        (Rational(1, 0), MATHICS3_COMPLEX_INFINITY, None),
+        (MATHICS3_COMPLEX_I, MATHICS3_COMPLEX_I, None),
+        (
+            SymbolI,
+            MATHICS3_COMPLEX_I,
+            (
+                "System`I evaluates to Complex[0,1] in the back and forward conversion. "
+                "This prevents an infinite recursion in evaluation"
+            ),
+        ),
+        # (Integer3**Rational(-1, 2), Rational(Integer1, Integer3)* (Integer3 ** (RationalOneHalf)), None ),
+    ],
+)
+def test_from_to_sympy_change(expr, result, msg):
+    """
+    Check if the conversion back and forward produces
+    the expected evaluation.
+    """
+    print([expr, result])
+    if msg:
+        assert result.sameQ(from_sympy(expr.to_sympy())), msg
+    else:
+        assert result.sameQ(from_sympy(expr.to_sympy()))
+
+
+def test_convert_sympy_singletons():
+    """
+    Check conversions between singleton symbols in
+    SymPy and Mathics Symbols.
+    """
+    for key, val in sympy_singleton_to_mathics.items():
+        print("equivalence", key, "<->", val)
+        if key is not None:
+            res = from_sympy(key)
+            print("  ->  ", res)
+            assert from_sympy(key).sameQ(val)
+
+            res = val.to_sympy()
+            print(res, "  <-  ")
+            assert res is key