diff --git a/mathics/builtin/arithfns/basic.py b/mathics/builtin/arithfns/basic.py
index 07ce6c364..439a0fdac 100644
--- a/mathics/builtin/arithfns/basic.py
+++ b/mathics/builtin/arithfns/basic.py
@@ -31,7 +31,7 @@
     A_READ_PROTECTED,
 )
 from mathics.core.convert.expression import to_expression
-from mathics.core.convert.sympy import from_sympy, to_sympy
+from mathics.core.convert.sympy import eval_sympy, from_sympy, to_sympy
 from mathics.core.expression import Expression
 from mathics.core.list import ListExpression
 from mathics.core.symbols import (
@@ -347,7 +347,7 @@ def negate(item):  # -> Expression (see FIXME below)
                 else:
                     return Expression(SymbolTimes, IntegerM1, *item.elements)
             elif isinstance(item, Number):
-                return from_sympy(to_sympy(-item))
+                return eval_sympy(-item)
             else:
                 return Expression(SymbolTimes, IntegerM1, item)
 
@@ -384,6 +384,7 @@ def is_negative(value) -> bool:
     def eval(self, items, evaluation):
         "Plus[items___]"
         items_tuple = numerify(items, evaluation).get_sequence()
+
         return eval_Plus(*items_tuple)
 
 
@@ -838,5 +839,6 @@ def format_outputform(self, items, evaluation):
 
     def eval(self, items, evaluation):
         "Times[items___]"
-        items = numerify(items, evaluation).get_sequence()
-        return eval_Times(*items)
+        items_tuple = numerify(items, evaluation).get_sequence()
+
+        return eval_Times(*items_tuple)
diff --git a/mathics/builtin/arithmetic.py b/mathics/builtin/arithmetic.py
index 6d73aa2b4..ab5783ba6 100644
--- a/mathics/builtin/arithmetic.py
+++ b/mathics/builtin/arithmetic.py
@@ -50,6 +50,7 @@
     from_sympy,
     sympy_symbol_prefix,
     to_sympy,
+    to_sympy_with_kwargs,
 )
 from mathics.core.element import BaseElement
 from mathics.core.evaluation import Evaluation
@@ -970,16 +971,16 @@ def to_sympy(self, expr, **kwargs):
                 elif cond is SymbolFalse:
                     sympy_cond = False
             if sympy_cond is None:
-                sympy_cond = to_sympy(cond, **kwargs)
+                sympy_cond = to_sympy_with_kwargs(cond, **kwargs)
                 if not (sympy_cond.is_Relational or sympy_cond.is_Boolean):
                     return
 
-            sympy_cases.append((to_sympy(then, **kwargs), sympy_cond))
+            sympy_cases.append((to_sympy_with_kwargs(then, **kwargs), sympy_cond))
 
         if len(elements) == 2:  # default case
-            sympy_cases.append((to_sympy(elements[1], **kwargs), True))
+            sympy_cases.append((to_sympy_with_kwargs(elements[1], **kwargs), True))
         else:
-            sympy_cases.append((to_sympy(Integer0, **kwargs), True))
+            sympy_cases.append((to_sympy_with_kwargs(Integer0, **kwargs), True))
 
         return sympy.Piecewise(*sympy_cases)
 
@@ -1061,10 +1062,10 @@ def to_sympy(self, expr, **kwargs):
             try:
                 e_kwargs = kwargs.copy()
                 e_kwargs["convert_all_global_functions"] = True
-                e = to_sympy(expr.elements[0], **e_kwargs)
-                i = to_sympy(index.elements[0], **kwargs)
-                start = to_sympy(index.elements[1], **kwargs)
-                stop = to_sympy(index.elements[2], **kwargs)
+                e = to_sympy_with_kwargs(expr.elements[0], **e_kwargs)
+                i = to_sympy_with_kwargs(index.elements[0], **kwargs)
+                start = to_sympy_with_kwargs(index.elements[1], **kwargs)
+                stop = to_sympy_with_kwargs(index.elements[2], **kwargs)
 
                 return sympy.product(e, (i, start, stop))
             except ZeroDivisionError:
@@ -1510,7 +1511,7 @@ def to_sympy(self, expr, **kwargs) -> Optional[SympyExpression]:
             index = expr.elements[1]
             arg_kwargs = kwargs.copy()
             arg_kwargs["convert_all_global_functions"] = True
-            f_sympy = to_sympy(expr.elements[0], **arg_kwargs)
+            f_sympy = to_sympy_with_kwargs(expr.elements[0], **arg_kwargs)
             if f_sympy is None:
                 return
 
@@ -1518,7 +1519,7 @@ def to_sympy(self, expr, **kwargs) -> Optional[SympyExpression]:
 
             # Handle summation parameters: variable, min, max
             var_min_max = index.elements[:3]
-            bounds = [to_sympy(expr, **kwargs) for expr in var_min_max]
+            bounds = [to_sympy_with_kwargs(expr, **kwargs) for expr in var_min_max]
 
             if evaluation:
                 # Min and max might be Mathics expressions. If so, evaluate them.
@@ -1526,7 +1527,7 @@ def to_sympy(self, expr, **kwargs) -> Optional[SympyExpression]:
                     min_max_expr = var_min_max[i]
                     if not isinstance(expr, Symbol):
                         min_max_expr_eval = min_max_expr.evaluate(evaluation)
-                        value = to_sympy(min_max_expr_eval, **kwargs)
+                        value = to_sympy_with_kwargs(min_max_expr_eval, **kwargs)
                         bounds[i] = value
 
             # FIXME: The below tests on SympyExpression, but really the
diff --git a/mathics/builtin/base.py b/mathics/builtin/base.py
index e02f9a3a3..f445a19d8 100644
--- a/mathics/builtin/base.py
+++ b/mathics/builtin/base.py
@@ -22,7 +22,12 @@
 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, to_numeric_sympy_args, to_sympy
+from mathics.core.convert.sympy import (
+    from_sympy,
+    to_numeric_sympy_args,
+    to_sympy,
+    to_sympy_with_kwargs,
+)
 from mathics.core.definitions import Definition
 from mathics.core.evaluation import Evaluation
 from mathics.core.exceptions import MessageException
@@ -691,7 +696,7 @@ def eval_iter(self, expr, i, imin, imax, di, evaluation):
             whole_expr = to_expression(
                 self.get_name(), expr, ListExpression(i, imin, imax)
             )
-            sympy_expr = to_sympy(whole_expr, evaluation=evaluation)
+            sympy_expr = to_sympy_with_kwargs(whole_expr, evaluation=evaluation)
             if sympy_expr is None:
                 return None
 
@@ -964,7 +969,9 @@ def to_sympy(self, expr, **kwargs):
         try:
             if self.sympy_name:
                 elements = self.prepare_sympy(expr.elements)
-                sympy_args = [to_sympy(element, **kwargs) for element in elements]
+                sympy_args = [
+                    to_sympy_with_kwargs(element, **kwargs) for element in elements
+                ]
                 if None in sympy_args:
                     return None
                 sympy_function = self.get_sympy_function(elements)
diff --git a/mathics/builtin/compilation.py b/mathics/builtin/compilation.py
index 7b33f5437..7f3032c9b 100644
--- a/mathics/builtin/compilation.py
+++ b/mathics/builtin/compilation.py
@@ -176,7 +176,12 @@ def to_sympy(self, *args, **kwargs):
         raise NotImplementedError
 
     def __hash__(self):
-        return hash(("CompiledCode", ctypes.addressof(self.cfunc)))  # XXX hack
+        try:
+            address = ctypes.addressof(self.cfunc)
+        except TypeError:
+            address = id(self.cfunc)
+
+        return hash(("CompiledCode", address))  # XXX hack
 
     def atom_to_boxes(self, f, evaluation: Evaluation):
         return CompiledCodeBox(String(self.__str__()), evaluation=evaluation)
diff --git a/mathics/builtin/numbers/calculus.py b/mathics/builtin/numbers/calculus.py
index 2737ae4b2..8b73e6328 100644
--- a/mathics/builtin/numbers/calculus.py
+++ b/mathics/builtin/numbers/calculus.py
@@ -57,6 +57,7 @@
     from_sympy,
     sympy_symbol_prefix,
     to_sympy,
+    to_sympy_with_kwargs,
 )
 from mathics.core.evaluation import Evaluation
 from mathics.core.expression import Expression
@@ -562,7 +563,7 @@ class DiscreteLimit(Builtin):
     def eval(self, f, n, n0, evaluation: Evaluation, options: dict = {}):
         "DiscreteLimit[f_, n_->n0_, OptionsPattern[DiscreteLimit]]"
 
-        f = to_sympy(f, convert_all_global_functions=True)
+        f = to_sympy_with_kwargs(f, convert_all_global_functions=True)
         n = to_sympy(n)
         n0 = to_sympy(n0)
 
diff --git a/mathics/builtin/numbers/diffeqns.py b/mathics/builtin/numbers/diffeqns.py
index 94860a942..b56a8c026 100644
--- a/mathics/builtin/numbers/diffeqns.py
+++ b/mathics/builtin/numbers/diffeqns.py
@@ -7,7 +7,7 @@
 import sympy
 
 from mathics.builtin.base import Builtin
-from mathics.core.convert.sympy import from_sympy, to_sympy
+from mathics.core.convert.sympy import from_sympy, to_sympy_with_kwargs
 from mathics.core.evaluation import Evaluation
 from mathics.core.expression import Expression
 from mathics.core.list import ListExpression
@@ -153,8 +153,8 @@ def eval(self, eqn, y, x, evaluation: Evaluation):
         f_name = func.get_head_name()
 
         conversion_args = {"converted_functions": set([f_name])}
-        sym_func = to_sympy(func, **conversion_args)
-        sym_eq = to_sympy(eqn, **conversion_args)
+        sym_func = to_sympy_with_kwargs(func, **conversion_args)
+        sym_eq = to_sympy_with_kwargs(eqn, **conversion_args)
 
         # XXX when sympy adds support for higher-order PDE we will have to
         # change this to a tuple of solvefuns
diff --git a/mathics/builtin/numbers/numbertheory.py b/mathics/builtin/numbers/numbertheory.py
index 169d31e07..9b0710f25 100644
--- a/mathics/builtin/numbers/numbertheory.py
+++ b/mathics/builtin/numbers/numbertheory.py
@@ -18,7 +18,12 @@
 )
 from mathics.core.convert.expression import to_mathics_list
 from mathics.core.convert.python import from_bool, from_python
-from mathics.core.convert.sympy import SympyPrime, from_sympy, to_sympy
+from mathics.core.convert.sympy import (
+    SympyPrime,
+    from_sympy,
+    to_sympy,
+    to_sympy_with_kwargs,
+)
 from mathics.core.evaluation import Evaluation
 from mathics.core.expression import Expression
 from mathics.core.list import ListExpression
@@ -612,7 +617,7 @@ def eval(self, n, evaluation: Evaluation):
 
     def to_sympy(self, expr, **kwargs):
         if expr.has_form("Prime", 1):
-            return SympyPrime(to_sympy(expr.elements[0], **kwargs))
+            return SympyPrime(to_sympy_with_kwargs(expr.elements[0], **kwargs))
 
 
 class PrimePi(SympyFunction):
diff --git a/mathics/builtin/optimization.py b/mathics/builtin/optimization.py
index d9454e851..170841b9d 100644
--- a/mathics/builtin/optimization.py
+++ b/mathics/builtin/optimization.py
@@ -22,7 +22,7 @@
 from mathics.core.atoms import IntegerM1
 from mathics.core.attributes import A_CONSTANT, A_PROTECTED, A_READ_PROTECTED
 from mathics.core.convert.python import from_python
-from mathics.core.convert.sympy import from_sympy, to_sympy
+from mathics.core.convert.sympy import eval_sympy, from_sympy, to_sympy
 from mathics.core.evaluation import Evaluation
 from mathics.core.expression import Expression
 from mathics.core.list import ListExpression
@@ -79,7 +79,7 @@ def eval_constraints(self, f, vars, evaluation: Evaluation):
         "Maximize[f_List, vars_]"
 
         constraints = [function for function in f.elements]
-        constraints[0] = from_sympy(to_sympy(constraints[0]) * IntegerM1)
+        constraints[0] = eval_sympy(constraints[0] * IntegerM1)
 
         dual_solutions = (
             Expression(SymbolMinimize, constraints, vars).evaluate(evaluation).elements
diff --git a/mathics/builtin/recurrence.py b/mathics/builtin/recurrence.py
index 53d07cd02..01c8308d5 100644
--- a/mathics/builtin/recurrence.py
+++ b/mathics/builtin/recurrence.py
@@ -15,7 +15,12 @@
 from mathics.builtin.base import Builtin
 from mathics.core.atoms import IntegerM1
 from mathics.core.attributes import A_CONSTANT
-from mathics.core.convert.sympy import from_sympy, sympy_symbol_prefix, to_sympy
+from mathics.core.convert.sympy import (
+    from_sympy,
+    sympy_symbol_prefix,
+    to_sympy,
+    to_sympy_with_kwargs,
+)
 from mathics.core.evaluation import Evaluation
 from mathics.core.expression import Expression
 from mathics.core.list import ListExpression
@@ -138,7 +143,9 @@ def is_relation(eqn):
             SymbolPlus, left, Expression(SymbolTimes, IntegerM1, right)
         ).evaluate(evaluation)
 
-        sym_eq = to_sympy(relation, converted_functions={func.get_head_name()})
+        sym_eq = to_sympy_with_kwargs(
+            relation, converted_functions={func.get_head_name()}
+        )
         if sym_eq is None:
             return
         sym_n = sympy.core.symbols(str(sympy_symbol_prefix + n.name))
diff --git a/mathics/builtin/testing_expressions/equality_inequality.py b/mathics/builtin/testing_expressions/equality_inequality.py
index e99456001..72d5463e5 100644
--- a/mathics/builtin/testing_expressions/equality_inequality.py
+++ b/mathics/builtin/testing_expressions/equality_inequality.py
@@ -18,7 +18,7 @@
     A_PROTECTED,
 )
 from mathics.core.convert.expression import to_expression, to_numeric_args
-from mathics.core.convert.sympy import to_sympy
+from mathics.core.convert.sympy import to_sympy, to_sympy_with_kwargs
 from mathics.core.expression import Expression
 from mathics.core.expression_predefined import (
     MATHICS3_COMPLEX_INFINITY,
@@ -155,8 +155,8 @@ def infty_equal(self, lhs, rhs, max_extra_prec=None) -> Optional[bool]:
 
     def sympy_equal(self, lhs, rhs, max_extra_prec=None) -> Optional[bool]:
         try:
-            lhs_sympy = to_sympy(lhs, evaluate=True, prec=COMPARE_PREC)
-            rhs_sympy = to_sympy(rhs, evaluate=True, prec=COMPARE_PREC)
+            lhs_sympy = to_sympy_with_kwargs(lhs, evaluate=True, prec=COMPARE_PREC)
+            rhs_sympy = to_sympy_with_kwargs(rhs, evaluate=True, prec=COMPARE_PREC)
         except NotImplementedError:
             return None
 
diff --git a/mathics/core/convert/sympy.py b/mathics/core/convert/sympy.py
index f9aeee632..831ecdff0 100644
--- a/mathics/core/convert/sympy.py
+++ b/mathics/core/convert/sympy.py
@@ -4,8 +4,8 @@
 Converts expressions from SymPy to Mathics expressions.
 Conversion to SymPy is handled directly in BaseElement descendants.
 """
-
-from typing import Optional, Type, Union
+from functools import lru_cache
+from typing import Optional, Tuple, Type, Union
 
 import sympy
 from sympy import Symbol as Sympy_Symbol, false as SympyFalse, true as SympyTrue
@@ -31,7 +31,7 @@
 )
 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.element import BaseElement, ElementsProperties
 from mathics.core.expression import Expression
 from mathics.core.expression_predefined import (
     MATHICS3_COMPLEX_INFINITY,
@@ -118,16 +118,31 @@ def is_Cn_expr(name) -> bool:
     return False
 
 
-def to_sympy(expr, **kwargs):
+def eval_sympy(expr):
+    sp = to_sympy(expr)
+    return None if sp is None else from_sympy(sp)
+
+
+def to_sympy(expr):
+    #    print("  to_sympy", expr)
+    if hasattr(expr, "to_sympy"):
+        return expr.to_sympy()
+    if isinstance(expr, Symbol):
+        return symbol_to_sympy(expr)
+    if isinstance(expr, Expression):
+        return expression_to_sympy(expr)
+
+
+def to_sympy_with_kwargs(expr, **kwargs):
     if hasattr(expr, "to_sympy"):
         return expr.to_sympy(**kwargs)
     if isinstance(expr, Symbol):
-        return symbol_to_sympy(expr, **kwargs)
+        return symbol_to_sympy_with_kwargs(expr, **kwargs)
     if isinstance(expr, Expression):
-        return expression_to_sympy(expr, **kwargs)
+        return expression_to_sympy_with_kwargs(expr, **kwargs)
 
 
-def to_sympy_matrix(data, **kwargs) -> Optional[sympy.MutableDenseMatrix]:
+def to_sympy_matrix(data) -> Optional[sympy.MutableDenseMatrix]:
     """Convert a Mathics matrix to one that can be used by Sympy.
     None is returned if we can't convert to a Sympy matrix.
     """
@@ -139,6 +154,33 @@ def to_sympy_matrix(data, **kwargs) -> Optional[sympy.MutableDenseMatrix]:
         return None
 
 
+class ExpressionWithToSympy(Expression):
+    """
+    This is an Expression that stores a Sympy expression from where
+    it was converted.
+
+    """
+
+    def __init__(
+        self,
+        head: BaseElement,
+        *elements: Tuple[BaseElement],
+        elements_properties: Optional[ElementsProperties] = None,
+        literal_values: Optional[tuple] = None,
+        sympy_expr=None,
+    ):
+        super().__init__(
+            head,
+            *elements,
+            elements_properties=elements_properties,
+            literal_values=literal_values,
+        )
+        self._sympy_form = sympy_expr
+
+    def to_sympy(self, **kwargs):
+        return self._sympy_form
+
+
 class SympyExpression(BasicSympy):
     is_Function = True
     nargs = None
@@ -223,13 +265,33 @@ def eval(cls, n):
                 pass
 
 
-def expression_to_sympy(expr: Expression, **kwargs):
+@lru_cache(maxsize=512)
+def expression_to_sympy(expr: Expression):
+    """
+    Convert `expr` to its sympy form.
+    """
+    lookup_name = expr.get_lookup_name()
+    # if lookup_name == "System`Pattern":
+    #    return None
+    builtin = mathics_to_sympy.get(lookup_name)
+    if builtin is not None:
+        sympy_expr = builtin.to_sympy(expr)
+        if sympy_expr is not None:
+            return sympy_expr
+    return SympyExpression(expr)
+
+
+def expression_to_sympy_with_kwargs(expr: Expression, **kwargs):
     """
     Convert `expr` to its sympy form.
     """
+    if len(kwargs) == 0:
+        return expression_to_sympy(expr)
 
     def as_sympy_function(expr, **kwargs) -> Optional[sympy.Function]:
-        sym_args = [to_sympy(element, **kwargs) for element in expr._elements]
+        sym_args = [
+            to_sympy_with_kwargs(element, **kwargs) for element in expr._elements
+        ]
 
         if None in sym_args:
             return None
@@ -262,7 +324,7 @@ def as_sympy_function(expr, **kwargs) -> Optional[sympy.Function]:
     return SympyExpression(expr)
 
 
-def symbol_to_sympy(symbol: Symbol, **kwargs) -> Sympy_Symbol:
+def symbol_to_sympy(symbol: Symbol) -> Sympy_Symbol:
     """
     Convert `symbol` to its sympy form.
     """
@@ -271,8 +333,25 @@ def symbol_to_sympy(symbol: Symbol, **kwargs) -> Sympy_Symbol:
     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)
+
+
+def symbol_to_sympy_with_kwargs(symbol: Symbol, **kwargs) -> Sympy_Symbol:
+    """
+    Convert `symbol` to its sympy form.
+    """
+    if len(kwargs) == 0:
+        return symbol_to_sympy(symbol)
+
+    result = mathics_to_sympy_singleton.get(symbol, None)
+    if result is not None:
+        return result
+
+    if symbol.name in kwargs.get("dummies", {}):
+        return None
 
     builtin = mathics_to_sympy.get(symbol.name)
     if builtin is None or not builtin.sympy_name or not builtin.is_constant():  # nopep8
@@ -299,6 +378,54 @@ def to_numeric_sympy_args(mathics_args: Type[BaseElement], evaluation) -> list:
     return sympy_args
 
 
+def from_sympy_O(expr):
+    if expr.args[0].func == sympy.core.power.Pow:
+        [var, power] = [from_sympy(arg) for arg in expr.args[0].args]
+        o = Expression(SymbolO, var)
+        return ExpressionWithToSympy(SymbolPower, o, power, sympy_expr=expr)
+    else:
+        return ExpressionWithToSympy(SymbolO, from_sympy(expr.args[0]), sympy_expr=expr)
+
+
+def from_sympy_add(expr):
+    if expr.is_number and all([x.is_Number for x in expr.as_real_imag()]):
+        # Hack to convert <Integer> * I to Complex[0, <Integer>]
+        try:
+            return Complex(*[from_sympy(arg) for arg in expr.as_real_imag()])
+        except ValueError:
+            # The exception happens if one of the components is infinity
+            pass
+
+    terms = sorted([from_sympy(arg) for arg in expr.args])
+    # TODO: Plug code the to reinforce the WL canonical form
+
+    return ExpressionWithToSympy(
+        SymbolPlus, *sorted([from_sympy(arg) for arg in expr.args]), sympy_expr=expr
+    )
+
+
+def from_sympy_c_roof_of(expr):
+    try:
+        e, i = expr.args
+    except ValueError:
+        return SymbolNull
+
+    try:
+        e = sympy.PurePoly(e)
+    except Exception:
+        pass
+
+    return Expression(SymbolRoot, from_sympy(e), Integer(i + 1))
+
+
+def from_sympy_lambda(expr):
+    vars = [
+        sympy.Symbol("%s%d" % (sympy_slot_prefix, index + 1))
+        for index in range(len(expr.variables))
+    ]
+    return Expression(SymbolFunction, from_sympy(expr(*vars)))
+
+
 def from_sympy_matrix(
     expr: Union[sympy.Matrix, sympy.ImmutableMatrix]
 ) -> ListExpression:
@@ -314,197 +441,12 @@ def from_sympy_matrix(
         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:
+@lru_cache(maxsize=512)
+def from_sympy_expression(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)
-    if isinstance(expr, int):
-        return Integer(expr)
-    if isinstance(expr, float):
-        return Real(expr)
-    if isinstance(expr, complex):
-        return Complex(Real(expr.real), Real(expr.imag))
-    if isinstance(expr, str):
-        return String(expr)
-    if expr is None:
-        return SymbolNull
-    if isinstance(expr, sympy.Matrix) or isinstance(expr, sympy.ImmutableMatrix):
-        return from_sympy_matrix(expr)
-    if isinstance(expr, sympy.MatPow):
-        return Expression(
-            SymbolMatrixPower, from_sympy(expr.base), from_sympy(expr.exp)
-        )
-    if expr.is_Atom:
-        name = None
-        if expr.is_Symbol:
-            name = str(expr)
-            if isinstance(expr, sympy.Dummy):
-                name = name + ("__Dummy_%d" % expr.dummy_index)
-                # Probably, this should be the value attribute
-                return Symbol(name, sympy_dummy=expr)
-            if is_Cn_expr(name):
-                return Expression(SymbolC, Integer(int(name[1:])))
-            if name.startswith(sympy_symbol_prefix):
-                name = name[len(sympy_symbol_prefix) :]
-            if name.startswith(sympy_slot_prefix):
-                index = name[len(sympy_slot_prefix) :]
-                return Expression(SymbolSlot, Integer(int(index)))
-        elif expr.is_NumberSymbol:
-            name = str(expr)
-        if name is not None:
-            builtin = sympy_to_mathics.get(name)
-            if builtin is not None:
-                name = builtin.get_name()
-            return Symbol(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 MATHICS3_NEG_INFINITY
-        elif isinstance(expr, sympy.core.numbers.ImaginaryUnit):
-            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 MATHICS3_INFINITY
-                elif numerator < 0:
-                    return MATHICS3_NEG_INFINITY
-                else:
-                    assert numerator == 0
-                    return SymbolIndeterminate
-            return Rational(numerator, denominator)
-        elif isinstance(expr, sympy.Float):
-            if expr._prec == FP_MANTISA_BINARY_DIGITS:
-                return MachineReal(float(expr))
-            return Real(expr)
-        elif isinstance(expr, sympy.core.numbers.NaN):
-            return SymbolIndeterminate
-        elif isinstance(expr, sympy.core.function.FunctionClass):
-            return Symbol(str(expr))
-        elif expr is sympy.true:
-            return SymbolTrue
-        elif expr is sympy.false:
-            return SymbolFalse
-
-    if expr.is_number and all([x.is_Number for x in expr.as_real_imag()]):
-        # Hack to convert <Integer> * I to Complex[0, <Integer>]
-        try:
-            return Complex(*[from_sympy(arg) for arg in expr.as_real_imag()])
-        except ValueError:
-            # The exception happens if one of the components is infinity
-            pass
-    if expr.is_Add:
-        return to_expression(
-            SymbolPlus, *sorted([from_sympy(arg) for arg in expr.args])
-        )
-    elif expr.is_Mul:
-        return to_expression(
-            SymbolTimes, *sorted([from_sympy(arg) for arg in expr.args])
-        )
-    elif expr.is_Pow:
-        return to_expression(SymbolPower, *[from_sympy(arg) for arg in expr.args])
-    elif expr.is_Equality:
-        return to_expression(SymbolEqual, *[from_sympy(arg) for arg in expr.args])
-
-    elif isinstance(expr, SympyExpression):
-        return expr.expr
-
-    elif isinstance(expr, sympy.Piecewise):
-        args = expr.args
-        return Expression(
-            SymbolPiecewise,
-            ListExpression(
-                *[
-                    to_mathics_list(from_sympy(case), from_sympy(cond))
-                    for case, cond in args
-                ]
-            ),
-        )
-
-    elif isinstance(expr, SympyPrime):
-        return Expression(SymbolPrime, from_sympy(expr.args[0]))
-    elif isinstance(expr, sympy.RootSum):
-        return Expression(SymbolRootSum, from_sympy(expr.poly), from_sympy(expr.fun))
-    elif isinstance(expr, sympy.PurePoly):
-        coeffs = expr.coeffs()
-        monoms = expr.monoms()
-        result = []
-        for coeff, monom in zip(coeffs, monoms):
-            factors = []
-            if coeff != 1:
-                factors.append(from_sympy(coeff))
-            for index, exp in enumerate(monom):
-                if exp != 0:
-                    slot = Expression(SymbolSlot, Integer(index + 1))
-                    if exp == 1:
-                        factors.append(slot)
-                    else:
-                        factors.append(Expression(SymbolPower, slot, from_sympy(exp)))
-            if factors:
-                result.append(Expression(SymbolTimes, *factors))
-            else:
-                result.append(Integer1)
-        return Expression(SymbolFunction, Expression(SymbolPlus, *result))
-    elif isinstance(expr, sympy.CRootOf):
-        try:
-            e, i = expr.args
-        except ValueError:
-            return SymbolNull
-
-        try:
-            e = sympy.PurePoly(e)
-        except Exception:
-            pass
-
-        return Expression(SymbolRoot, from_sympy(e), Integer(i + 1))
-    elif isinstance(expr, sympy.Lambda):
-        vars = [
-            sympy.Symbol("%s%d" % (sympy_slot_prefix, index + 1))
-            for index in range(len(expr.variables))
-        ]
-        return Expression(SymbolFunction, from_sympy(expr(*vars)))
-
-    elif expr.is_Function or isinstance(
+    if expr.is_Function or isinstance(
         expr, (sympy.Integral, sympy.Derivative, sympy.Sum, sympy.Product)
     ):
         if isinstance(expr, sympy.Integral):
@@ -532,7 +474,7 @@ def old_from_sympy(expr) -> BaseElement:
             if is_Cn_expr(name):
                 return Expression(
                     Expression(Symbol("C"), Integer(int(name[1:]))),
-                    *[from_sympy(arg) for arg in expr.args]
+                    *[from_sympy(arg) for arg in expr.args],
                 )
             if name.startswith(sympy_symbol_prefix):
                 name = name[len(sympy_symbol_prefix) :]
@@ -540,42 +482,7 @@ def old_from_sympy(expr) -> BaseElement:
         builtin = sympy_to_mathics.get(name)
         if builtin is not None:
             return builtin.from_sympy(name, args)
-        return Expression(Symbol(name), *args)
-
-    elif isinstance(expr, sympy.Tuple):
-        return to_mathics_list(*expr.args, elements_conversion_fn=from_sympy)
-
-    # elif isinstance(expr, sympy.Sum):
-    #    return Expression('Sum', )
-
-    elif isinstance(expr, sympy.LessThan):
-        return to_expression(
-            SymbolLessEqual, *expr.args, elements_conversion_fn=from_sympy
-        )
-    elif isinstance(expr, sympy.StrictLessThan):
-        return to_expression(SymbolLess, *expr.args, elements_conversion_fn=from_sympy)
-    elif isinstance(expr, sympy.GreaterThan):
-        return to_expression(
-            SymbolGreaterEqual, *expr.args, elements_conversion_fn=from_sympy
-        )
-    elif isinstance(expr, sympy.StrictGreaterThan):
-        return to_expression(
-            SymbolGreater, *expr.args, elements_conversion_fn=from_sympy
-        )
-    elif isinstance(expr, sympy.Unequality):
-        return to_expression(
-            SymbolUnequal, *expr.args, elements_conversion_fn=from_sympy
-        )
-    elif isinstance(expr, sympy.Equality):
-        return to_expression(SymbolEqual, *expr.args, elements_conversion_fn=from_sympy)
-
-    elif isinstance(expr, sympy.O):
-        if expr.args[0].func == sympy.core.power.Pow:
-            [var, power] = [from_sympy(arg) for arg in expr.args[0].args]
-            o = Expression(SymbolO, var)
-            return Expression(SymbolPower, o, power)
-        else:
-            return Expression(SymbolO, from_sympy(expr.args[0]))
+        return ExpressionWithToSympy(Symbol(name), *args, sympy_expr=expr)
     else:
         raise ValueError(
             "Unknown SymPy expression: {} (instance of {})".format(
@@ -584,4 +491,165 @@ def old_from_sympy(expr) -> BaseElement:
         )
 
 
-from_sympy = old_from_sympy
+def from_sympy_float(expr):
+    if expr._prec == FP_MANTISA_BINARY_DIGITS:
+        return MachineReal(float(expr))
+    return Real(expr)
+
+
+def from_sympy_mul(expr):
+    if expr.is_number and all([x.is_Number for x in expr.as_real_imag()]):
+        # Hack to convert <Integer> * I to Complex[0, <Integer>]
+        try:
+            return Complex(*[from_sympy(arg) for arg in expr.as_real_imag()])
+        except ValueError:
+            # The exception happens if one of the components is infinity
+            pass
+
+    factors = sorted([from_sympy(arg) for arg in expr.args])
+    # TODO: Plug code here to reinforce the WL canonical form.
+
+    return ExpressionWithToSympy(SymbolTimes, *factors, sympy_expr=expr)
+
+
+def from_sympy_piecewise(expr):
+    args = expr.args
+    return ExpressionWithToSympy(
+        SymbolPiecewise,
+        ListExpression(
+            *[
+                to_mathics_list(from_sympy(case), from_sympy(cond))
+                for case, cond in args
+            ]
+        ),
+        sympy_expr=expr,
+    )
+
+
+def from_sympy_pure_poly(expr):
+    coeffs = expr.coeffs()
+    monoms = expr.monoms()
+    result = []
+    for coeff, monom in zip(coeffs, monoms):
+        factors = []
+        if coeff != 1:
+            factors.append(from_sympy(coeff))
+        for index, exp in enumerate(monom):
+            if exp != 0:
+                slot = Expression(SymbolSlot, Integer(index + 1))
+                if exp == 1:
+                    factors.append(slot)
+                else:
+                    factors.append(Expression(SymbolPower, slot, from_sympy(exp)))
+        if factors:
+            result.append(Expression(SymbolTimes, *factors))
+        else:
+            result.append(Integer1)
+    return ExpressionWithToSympy(
+        SymbolFunction, Expression(SymbolPlus, *result), sympy_expr=expr
+    )
+
+
+def from_sympy_rational(expr):
+    numerator, denominator = map(int, expr.as_numer_denom())
+    if denominator == 0:
+        if numerator > 0:
+            return MATHICS3_INFINITY
+        elif numerator < 0:
+            return MATHICS3_NEG_INFINITY
+        else:
+            assert numerator == 0
+            return SymbolIndeterminate
+    return Rational(numerator, denominator)
+
+
+def from_sympy_symbol(s):
+    name = str(s)
+    if name.startswith(sympy_symbol_prefix):
+        return Symbol(name[len(sympy_symbol_prefix) :])
+    if name.startswith(sympy_slot_prefix):
+        index = name[len(sympy_slot_prefix) :]
+        return ExpressionWithToSympy(SymbolSlot, Integer(int(index)), sympy_expr=s)
+    if is_Cn_expr(name):
+        return ExpressionWithToSympy(SymbolC, Integer(int(name[1:])), sympy_expr=s)
+    return None
+
+
+def from_sympy_pow(expr):
+    if expr.is_number and all([x.is_Number for x in expr.as_real_imag()]):
+        # Hack to convert <Integer> * I to Complex[0, <Integer>]
+        try:
+            return Complex(*[from_sympy(arg) for arg in expr.as_real_imag()])
+        except ValueError:
+            # The exception happens if one of the components is infinity
+            pass
+    return ExpressionWithToSympy(
+        SymbolPower, *[from_sympy(arg) for arg in expr.args], sympy_expr=expr
+    )
+
+
+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),
+    SympyExpression: lambda expr: expr.expr,
+    SympyPrime: lambda expr: Expression(SymbolPrime, from_sympy(expr.args[0])),
+    sympy.CRootOf: from_sympy_c_roof_of,
+    sympy.Dummy: lambda expr: Symbol(f"SympyDummy`{str(expr)}$${expr.dummy_index}"),
+    sympy.Equality: lambda expr: to_expression(
+        SymbolEqual, *[from_sympy(arg) for arg in expr.args]
+    ),
+    sympy.GreaterThan: lambda expr: to_expression(
+        SymbolGreaterEqual, *expr.args, elements_conversion_fn=from_sympy
+    ),
+    sympy.ImmutableMatrix: from_sympy_matrix,
+    sympy.Lambda: from_sympy_lambda,
+    sympy.LessThan: lambda expr: to_expression(
+        SymbolLessEqual, *expr.args, elements_conversion_fn=from_sympy
+    ),
+    sympy.MatPow: lambda expr: Expression(
+        SymbolMatrixPower, from_sympy(expr.base), from_sympy(expr.exp)
+    ),
+    sympy.Matrix: from_sympy_matrix,
+    sympy.O: from_sympy_O,
+    sympy.Piecewise: from_sympy_piecewise,
+    sympy.PurePoly: from_sympy_pure_poly,
+    sympy.RootSum: lambda expr: Expression(
+        SymbolRootSum, from_sympy(expr.poly), from_sympy(expr.fun)
+    ),
+    sympy.StrictGreaterThan: lambda expr: to_expression(
+        SymbolGreater, *expr.args, elements_conversion_fn=from_sympy
+    ),
+    sympy.StrictLessThan: lambda expr: to_expression(
+        SymbolLess, *expr.args, elements_conversion_fn=from_sympy
+    ),
+    sympy.Tuple: lambda expr: to_mathics_list(*expr, elements_conversion_fn=from_sympy),
+    sympy.Unequality: lambda expr: to_expression(
+        SymbolUnequal, *expr.args, elements_conversion_fn=from_sympy
+    ),
+    sympy.core.add.Add: from_sympy_add,
+    sympy.core.mul.Mul: from_sympy_mul,
+    sympy.core.numbers.Float: from_sympy_float,
+    sympy.core.numbers.Integer: lambda x: Integer(x),
+    sympy.core.numbers.Rational: from_sympy_rational,
+    sympy.core.power.Pow: from_sympy_pow,
+    sympy.core.symbol.Symbol: from_sympy_symbol,
+}
+
+
+def from_sympy(expr) -> BaseElement:
+    """
+    converts a SymPy object to a Mathics element.
+    """
+    #    print("  from_sympy", expr)
+    if expr is None:
+        return SymbolNull
+    try:
+        return sympy_singleton_to_mathics[expr]
+    except (KeyError, TypeError):
+        pass
+
+    return sympy_conversion_by_type.get(type(expr), from_sympy_expression)(expr)
diff --git a/mathics/core/symbols.py b/mathics/core/symbols.py
index 41f842ddd..687cf6f32 100644
--- a/mathics/core/symbols.py
+++ b/mathics/core/symbols.py
@@ -360,7 +360,6 @@ class Symbol(Atom, NumericOperators, EvalMixin):
 
     name: str
     hash: str
-    sympy_dummy: Any
 
     # Dictionary of Symbols defined so far.
     # We use this for object uniqueness.
@@ -372,7 +371,7 @@ class Symbol(Atom, NumericOperators, EvalMixin):
 
     # __new__ instead of __init__ is used here because we want
     # to return the same object for a given "name" value.
-    def __new__(cls, name: str, sympy_dummy=None):
+    def __new__(cls, name: str):
         """
         Allocate an object ensuring that for a given ``name`` and ``cls`` we get back the same object,
         id(object) is the same and its object.__hash__() is the same.
@@ -402,18 +401,6 @@ def __new__(cls, name: str, sympy_dummy=None):
             # For example, this can happen with String constants.
 
             self.hash = hash((cls, name))
-
-            # TODO: revise how we convert sympy.Dummy
-            # symbols.
-            #
-            # In some cases, SymPy returns a sympy.Dummy
-            # object. It is converted to Mathics as a
-            # Symbol. However, we probably should have
-            # a different class for this kind of symbols.
-            # Also, sympy_dummy should be stored as the
-            # value attribute.
-            self.sympy_dummy = sympy_dummy
-
             self._short_name = strip_context(name)
 
         return self
@@ -422,7 +409,7 @@ def __eq__(self, other) -> bool:
         return self is other
 
     def __getnewargs__(self):
-        return (self.name, self.sympy_dummy)
+        return tuple([self.name])
 
     def __hash__(self) -> int:
         """
@@ -674,6 +661,9 @@ def __new__(cls, name, value):
             self.hash = hash((cls, name))
         return self
 
+    def __getnewargs__(self):
+        return tuple([self.name, self._value])
+
     @property
     def is_literal(self) -> bool:
         """
diff --git a/mathics/eval/arithmetic.py b/mathics/eval/arithmetic.py
index 84b3ab86d..125d56fa8 100644
--- a/mathics/eval/arithmetic.py
+++ b/mathics/eval/arithmetic.py
@@ -49,6 +49,9 @@
 RealOne = Real(1.0)
 
 
+use_sympy_for_arithmetic = False
+
+
 # This cache might not be used that much.
 @lru_cache()
 def call_mpmath(
@@ -347,7 +350,15 @@ def eval_mpmath_function(
             return call_mpmath(mpmath_function, tuple(mpmath_args), prec)
 
 
-def eval_Plus(*items: BaseElement) -> BaseElement:
+def eval_Plus_sympy(*items: BaseElement) -> BaseElement:
+    items_tuple = tuple(to_sympy(it) for it in items)
+    if any(it is None for it in items_tuple):
+        return
+    res = sympy.Add(*items_tuple)
+    return from_sympy(res) if res is not None else None
+
+
+def eval_Plus_native(*items: BaseElement) -> BaseElement:
     "evaluate Plus for general elements"
     numbers, items_tuple = segregate_numbers_from_sorted_list(*items)
     elements = []
@@ -559,7 +570,15 @@ def eval_Power_inexact(base: Number, exp: Number) -> BaseElement:
                 return from_mpmath(number, prec)
 
 
-def eval_Times(*items: BaseElement) -> BaseElement:
+def eval_Times_sympy(*items: BaseElement) -> BaseElement:
+    items_tuple = tuple(to_sympy(it) for it in items)
+    if any(it is None for it in items_tuple):
+        return
+    res = sympy.Mul(*items_tuple)
+    return from_sympy(res) if res is not None else None
+
+
+def eval_Times_native(*items: BaseElement) -> BaseElement:
     elements = []
     numbers = []
     # find numbers and simplify Times -> Power
@@ -898,3 +917,11 @@ def to_inexact_value(expr: BaseElement) -> BaseElement:
         except Exception:
             pass
     return None
+
+
+if use_sympy_for_arithmetic:
+    eval_Plus = eval_Plus_sympy
+    eval_Times = eval_Times_sympy
+else:
+    eval_Plus = eval_Plus_native
+    eval_Times = eval_Times_native
diff --git a/mathics/session.py b/mathics/session.py
index ae23eaa6c..d9baf113e 100644
--- a/mathics/session.py
+++ b/mathics/session.py
@@ -14,8 +14,11 @@
 
 from mathics.core.definitions import Definitions, autoload_files
 from mathics.core.evaluation import Evaluation
+from mathics.core.load_builtin import import_and_load_builtins
 from mathics.core.parser import MathicsSingleLineFeeder, parse
 
+import_and_load_builtins()
+
 
 def load_default_settings_files(
     definitions: Definitions, load_cli_settings: bool = True
diff --git a/mathics/timing.py b/mathics/timing.py
index f8ab01ac7..56ef7ad8d 100644
--- a/mathics/timing.py
+++ b/mathics/timing.py
@@ -70,6 +70,7 @@ def show_lru_cache_statistics():
     from mathics.builtin.base import run_sympy
     from mathics.core.atoms import Integer, Rational
     from mathics.core.convert.mpmath import from_mpmath
+    from mathics.core.convert.sympy import expression_to_sympy, from_sympy_expression
     from mathics.eval.arithmetic import call_mpmath
 
     print(f"Integer             {len(Integer._integers)}")
@@ -78,5 +79,6 @@ def show_lru_cache_statistics():
     print(f"log_n_b             {log_n_b.cache_info()}")
     print(f"from_mpmath         {from_mpmath.cache_info()}")
     print(f"get_mpmath_function {_MPMathFunction.get_mpmath_function.cache_info()}")
-
+    print(f"expression_to_sympy            {expression_to_sympy.cache_info()}")
+    print(f"from_sympy_expression          {from_sympy_expression.cache_info()}")
     print(f"run_sympy           {run_sympy.cache_info()}")
diff --git a/test/core/test_sympy_python_convert.py b/test/core/test_sympy_python_convert.py
index 9f20b3e8b..763555104 100644
--- a/test/core/test_sympy_python_convert.py
+++ b/test/core/test_sympy_python_convert.py
@@ -17,7 +17,7 @@
     String,
 )
 from mathics.core.convert.python import from_python
-from mathics.core.convert.sympy import from_sympy, to_sympy
+from mathics.core.convert.sympy import from_sympy, to_sympy, to_sympy_with_kwargs
 from mathics.core.expression import Expression
 from mathics.core.list import ListExpression
 from mathics.core.symbols import Symbol, SymbolPlus
@@ -32,7 +32,7 @@
 
 class SympyConvert(unittest.TestCase):
     def compare_to_sympy(self, mathics_expr, sympy_expr, **kwargs):
-        to_sympy(mathics_expr, **kwargs) == sympy_expr
+        to_sympy_with_kwargs(mathics_expr, **kwargs) == sympy_expr
 
     def compare_to_mathics(self, mathics_expr, sympy_expr, **kwargs):
         mathics_expr == from_sympy(sympy_expr, **kwargs)