diff --git a/SYMBOLS_MANIFEST.txt b/SYMBOLS_MANIFEST.txt
index d3ecf4204..43c6e2162 100644
--- a/SYMBOLS_MANIFEST.txt
+++ b/SYMBOLS_MANIFEST.txt
@@ -574,6 +574,7 @@ System`LessEqual
 System`LetterCharacter
 System`LetterNumber
 System`LetterQ
+System`LeviCivitaTensor
 System`Level
 System`LevelQ
 System`LightBlue
diff --git a/mathics/builtin/distance/clusters.py b/mathics/builtin/distance/clusters.py
index ebd904540..382fce982 100644
--- a/mathics/builtin/distance/clusters.py
+++ b/mathics/builtin/distance/clusters.py
@@ -35,6 +35,7 @@
 )
 from mathics.eval.nevaluator import eval_N
 from mathics.eval.parts import walk_levels
+from mathics.eval.tensors import get_default_distance
 
 
 class _LazyDistances(LazyDistances):
@@ -139,8 +140,6 @@ def _cluster(self, p, k, mode, evaluation, options, expr):
             options, "DistanceFunction", evaluation
         )
         if distance_function_string == "Automatic":
-            from mathics.builtin.tensors import get_default_distance
-
             distance_function = get_default_distance(dist_p)
             if distance_function is None:
                 name_of_builtin = strip_context(self.get_name())
@@ -462,8 +461,6 @@ def eval(
             options, "DistanceFunction", evaluation
         )
         if distance_function_string == "Automatic":
-            from mathics.builtin.tensors import get_default_distance
-
             distance_function = get_default_distance(dist_p)
             if distance_function is None:
                 evaluation.message(
diff --git a/mathics/builtin/tensors.py b/mathics/builtin/tensors.py
index 831a4de5a..bad1a83f4 100644
--- a/mathics/builtin/tensors.py
+++ b/mathics/builtin/tensors.py
@@ -19,70 +19,17 @@
 """
 
 
-from sympy.combinatorics import Permutation
-from sympy.utilities.iterables import permutations
-
-from mathics.core.atoms import Integer, String
+from mathics.core.atoms import Integer
 from mathics.core.attributes import A_FLAT, A_ONE_IDENTITY, A_PROTECTED
 from mathics.core.builtin import BinaryOperator, Builtin
-from mathics.core.convert.python import from_python
 from mathics.core.evaluation import Evaluation
-from mathics.core.expression import Expression
 from mathics.core.list import ListExpression
-from mathics.core.symbols import Atom, Symbol, SymbolFalse, SymbolTrue
-from mathics.core.systemsymbols import SymbolRule, SymbolSparseArray
-from mathics.eval.parts import get_part
-
-
-def get_default_distance(p):
-    if all(q.is_numeric() for q in p):
-        return Symbol("SquaredEuclideanDistance")
-    elif all(q.get_head_name() == "System`List" for q in p):
-        dimensions = [get_dimensions(q) for q in p]
-        if len(dimensions) < 1:
-            return None
-        d0 = dimensions[0]
-        if not all(d == d0 for d in dimensions[1:]):
-            return None
-        if len(dimensions[0]) == 1:  # vectors?
-
-            def is_boolean(x):
-                return x.get_head_name() == "System`Symbol" and x in (
-                    SymbolTrue,
-                    SymbolFalse,
-                )
-
-            if all(all(is_boolean(e) for e in q.elements) for q in p):
-                return Symbol("JaccardDissimilarity")
-        return Symbol("SquaredEuclideanDistance")
-    elif all(isinstance(q, String) for q in p):
-        return Symbol("EditDistance")
-    else:
-        from mathics.builtin.colors.color_directives import expression_to_color
-
-        if all(expression_to_color(q) is not None for q in p):
-            return Symbol("ColorDistance")
-
-        return None
-
-
-def get_dimensions(expr, head=None):
-    if isinstance(expr, Atom):
-        return []
-    else:
-        if head is not None and not expr.head.sameQ(head):
-            return []
-        sub_dim = None
-        sub = []
-        for element in expr.elements:
-            sub = get_dimensions(element, expr.head)
-            if sub_dim is None:
-                sub_dim = sub
-            else:
-                if sub_dim != sub:
-                    sub = []
-                    break
-        return [len(expr.elements)] + sub
+from mathics.eval.tensors import (
+    eval_Inner,
+    eval_LeviCivitaTensor,
+    eval_Outer,
+    get_dimensions,
+)
 
 
 class ArrayDepth(Builtin):
@@ -220,46 +167,7 @@ class Inner(Builtin):
     def eval(self, f, list1, list2, g, evaluation: Evaluation):
         "Inner[f_, list1_, list2_, g_]"
 
-        m = get_dimensions(list1)
-        n = get_dimensions(list2)
-
-        if not m or not n:
-            evaluation.message("Inner", "normal")
-            return
-        if list1.get_head() != list2.get_head():
-            evaluation.message("Inner", "heads", list1.get_head(), list2.get_head())
-            return
-        if m[-1] != n[0]:
-            evaluation.message("Inner", "incom", m[-1], len(m), list1, n[0], list2)
-            return
-
-        head = list1.get_head()
-        inner_dim = n[0]
-
-        def rec(i_cur, j_cur, i_rest, j_rest):
-            evaluation.check_stopped()
-            if i_rest:
-                elements = []
-                for i in range(1, i_rest[0] + 1):
-                    elements.append(rec(i_cur + [i], j_cur, i_rest[1:], j_rest))
-                return Expression(head, *elements)
-            elif j_rest:
-                elements = []
-                for j in range(1, j_rest[0] + 1):
-                    elements.append(rec(i_cur, j_cur + [j], i_rest, j_rest[1:]))
-                return Expression(head, *elements)
-            else:
-
-                def summand(i):
-                    part1 = get_part(list1, i_cur + [i])
-                    part2 = get_part(list2, [i] + j_cur)
-                    return Expression(f, part1, part2)
-
-                part = Expression(g, *[summand(i) for i in range(1, inner_dim + 1)])
-                # cur_expr.elements.append(part)
-                return part
-
-        return rec([], [], m[:-1], n[1:])
+        return eval_Inner(f, list1, list2, g, evaluation)
 
 
 class Outer(Builtin):
@@ -278,10 +186,21 @@ class Outer(Builtin):
     Outer product of two matrices:
     >> Outer[Times, {{a, b}, {c, d}}, {{1, 2}, {3, 4}}]
      = {{{{a, 2 a}, {3 a, 4 a}}, {{b, 2 b}, {3 b, 4 b}}}, {{{c, 2 c}, {3 c, 4 c}}, {{d, 2 d}, {3 d, 4 d}}}}
+    
+    Outer product of two sparse arrays:
+    >> Outer[Times, SparseArray[{{1, 2} -> a, {2, 1} -> b}], SparseArray[{{1, 2} -> c, {2, 1} -> d}]]
+     = SparseArray[Automatic, {2, 2, 2, 2}, 0, {{1, 2, 1, 2} -> a c, {1, 2, 2, 1} -> a d, {2, 1, 1, 2} -> b c, {2, 1, 2, 1} -> b d}]
 
     'Outer' of multiple lists:
     >> Outer[f, {a, b}, {x, y, z}, {1, 2}]
      = {{{f[a, x, 1], f[a, x, 2]}, {f[a, y, 1], f[a, y, 2]}, {f[a, z, 1], f[a, z, 2]}}, {{f[b, x, 1], f[b, x, 2]}, {f[b, y, 1], f[b, y, 2]}, {f[b, z, 1], f[b, z, 2]}}}
+    
+    'Outer' converts input sparse arrays to lists if f=!=Times, or if the input is a mixture of sparse arrays and lists:
+    >> Outer[f, SparseArray[{{1, 2} -> a, {2, 1} -> b}], SparseArray[{{1, 2} -> c, {2, 1} -> d}]]
+     = {{{{f[0, 0], f[0, c]}, {f[0, d], f[0, 0]}}, {{f[a, 0], f[a, c]}, {f[a, d], f[a, 0]}}}, {{{f[b, 0], f[b, c]}, {f[b, d], f[b, 0]}}, {{f[0, 0], f[0, c]}, {f[0, d], f[0, 0]}}}}
+    
+    >> Outer[Times, SparseArray[{{1, 2} -> a, {2, 1} -> b}], {c, d}]
+     = {{{0, 0}, {a c, a d}}, {{b c, b d}, {0, 0}}}
 
     Arrays can be ragged:
     >> Outer[Times, {{1, 2}}, {{a, b}, {c, d, e}}]
@@ -304,32 +223,7 @@ class Outer(Builtin):
     def eval(self, f, lists, evaluation: Evaluation):
         "Outer[f_, lists__]"
 
-        lists = lists.get_sequence()
-        head = None
-        for list in lists:
-            if isinstance(list, Atom):
-                evaluation.message("Outer", "normal")
-                return
-            if head is None:
-                head = list.head
-            elif not list.head.sameQ(head):
-                evaluation.message("Outer", "heads", head, list.head)
-                return
-
-        def rec(item, rest_lists, current):
-            evaluation.check_stopped()
-            if isinstance(item, Atom) or not item.head.sameQ(head):
-                if rest_lists:
-                    return rec(rest_lists[0], rest_lists[1:], current + [item])
-                else:
-                    return Expression(f, *(current + [item]))
-            else:
-                elements = []
-                for element in item.elements:
-                    elements.append(rec(element, rest_lists, current))
-                return Expression(head, *elements)
-
-        return rec(lists[0], lists[1:], [])
+        return eval_Outer(f, lists, evaluation)
 
 
 class RotationTransform(Builtin):
@@ -548,17 +442,4 @@ class LeviCivitaTensor(Builtin):
     def eval(self, d, type, evaluation: Evaluation):
         "LeviCivitaTensor[d_Integer, type_]"
 
-        if isinstance(d, Integer) and type == SymbolSparseArray:
-            d = d.get_int_value()
-            perms = list(permutations([i for i in range(1, d + 1)]))
-            rules = [
-                Expression(
-                    SymbolRule,
-                    from_python(p),
-                    from_python(Permutation.from_sequence(p).signature()),
-                )
-                for p in perms
-            ]
-            return Expression(
-                SymbolSparseArray, from_python(rules), from_python([d] * d)
-            )
+        return eval_LeviCivitaTensor(d, type)
diff --git a/mathics/eval/tensors.py b/mathics/eval/tensors.py
new file mode 100644
index 000000000..596c792b0
--- /dev/null
+++ b/mathics/eval/tensors.py
@@ -0,0 +1,244 @@
+import itertools
+
+from sympy.combinatorics import Permutation
+from sympy.utilities.iterables import permutations
+
+from mathics.core.atoms import Integer, Integer0, Integer1, String
+from mathics.core.convert.python import from_python
+from mathics.core.evaluation import Evaluation
+from mathics.core.expression import Expression
+from mathics.core.list import ListExpression
+from mathics.core.symbols import (
+    Atom,
+    Symbol,
+    SymbolFalse,
+    SymbolList,
+    SymbolTimes,
+    SymbolTrue,
+)
+from mathics.core.systemsymbols import (
+    SymbolAutomatic,
+    SymbolNormal,
+    SymbolRule,
+    SymbolSparseArray,
+)
+from mathics.eval.parts import get_part
+
+
+def get_default_distance(p):
+    if all(q.is_numeric() for q in p):
+        return Symbol("SquaredEuclideanDistance")
+    elif all(q.get_head_name() == "System`List" for q in p):
+        dimensions = [get_dimensions(q) for q in p]
+        if len(dimensions) < 1:
+            return None
+        d0 = dimensions[0]
+        if not all(d == d0 for d in dimensions[1:]):
+            return None
+        if len(dimensions[0]) == 1:  # vectors?
+
+            def is_boolean(x):
+                return x.get_head_name() == "System`Symbol" and x in (
+                    SymbolTrue,
+                    SymbolFalse,
+                )
+
+            if all(all(is_boolean(e) for e in q.elements) for q in p):
+                return Symbol("JaccardDissimilarity")
+        return Symbol("SquaredEuclideanDistance")
+    elif all(isinstance(q, String) for q in p):
+        return Symbol("EditDistance")
+    else:
+        from mathics.builtin.colors.color_directives import expression_to_color
+
+        if all(expression_to_color(q) is not None for q in p):
+            return Symbol("ColorDistance")
+
+        return None
+
+
+def get_dimensions(expr, head=None):
+    if isinstance(expr, Atom):
+        return []
+    else:
+        if head is not None and not expr.head.sameQ(head):
+            return []
+        sub_dim = None
+        sub = []
+        for element in expr.elements:
+            sub = get_dimensions(element, expr.head)
+            if sub_dim is None:
+                sub_dim = sub
+            else:
+                if sub_dim != sub:
+                    sub = []
+                    break
+        return [len(expr.elements)] + sub
+
+
+def eval_Inner(f, list1, list2, g, evaluation: Evaluation):
+    "Evaluates recursively the inner product of list1 and list2"
+
+    m = get_dimensions(list1)
+    n = get_dimensions(list2)
+
+    if not m or not n:
+        evaluation.message("Inner", "normal")
+        return
+    if list1.get_head() != list2.get_head():
+        evaluation.message("Inner", "heads", list1.get_head(), list2.get_head())
+        return
+    if m[-1] != n[0]:
+        evaluation.message("Inner", "incom", m[-1], len(m), list1, n[0], list2)
+        return
+
+    head = list1.get_head()
+    inner_dim = n[0]
+
+    def rec(i_cur, j_cur, i_rest, j_rest):
+        evaluation.check_stopped()
+        if i_rest:
+            elements = []
+            for i in range(1, i_rest[0] + 1):
+                elements.append(rec(i_cur + [i], j_cur, i_rest[1:], j_rest))
+            return Expression(head, *elements)
+        elif j_rest:
+            elements = []
+            for j in range(1, j_rest[0] + 1):
+                elements.append(rec(i_cur, j_cur + [j], i_rest, j_rest[1:]))
+            return Expression(head, *elements)
+        else:
+
+            def summand(i):
+                part1 = get_part(list1, i_cur + [i])
+                part2 = get_part(list2, [i] + j_cur)
+                return Expression(f, part1, part2)
+
+            part = Expression(g, *[summand(i) for i in range(1, inner_dim + 1)])
+            # cur_expr.elements.append(part)
+            return part
+
+    return rec([], [], m[:-1], n[1:])
+
+
+def eval_Outer(f, lists, evaluation: Evaluation):
+    "Evaluates recursively the outer product of lists"
+
+    # If f=!=Times, or lists contain both SparseArray and List, then convert all SparseArrays to Lists
+    lists = lists.get_sequence()
+    head = None
+    sparse_to_list = f != SymbolTimes
+    contain_sparse = False
+    contain_list = False
+    for _list in lists:
+        if _list.head.sameQ(SymbolSparseArray):
+            contain_sparse = True
+        if _list.head.sameQ(SymbolList):
+            contain_list = True
+        sparse_to_list = sparse_to_list or (contain_sparse and contain_list)
+        if sparse_to_list:
+            break
+    if sparse_to_list:
+        new_lists = []
+    for _list in lists:
+        if isinstance(_list, Atom):
+            evaluation.message("Outer", "normal")
+            return
+        if sparse_to_list:
+            if _list.head.sameQ(SymbolSparseArray):
+                _list = Expression(SymbolNormal, _list).evaluate(evaluation)
+            new_lists.append(_list)
+        if head is None:
+            head = _list.head
+        elif not _list.head.sameQ(head):
+            evaluation.message("Outer", "heads", head, _list.head)
+            return
+    if sparse_to_list:
+        lists = new_lists
+
+    def rec(item, rest_lists, current):
+        evaluation.check_stopped()
+        if isinstance(item, Atom) or not item.head.sameQ(head):
+            if rest_lists:
+                return rec(rest_lists[0], rest_lists[1:], current + [item])
+            else:
+                return Expression(f, *(current + [item]))
+        else:
+            elements = []
+            for element in item.elements:
+                elements.append(rec(element, rest_lists, current))
+            return Expression(head, *elements)
+
+    def rec_sparse(item, rest_lists, current):
+        evaluation.check_stopped()
+        if isinstance(item, tuple):  # (rules)
+            elements = []
+            for element in item:
+                elements.extend(rec_sparse(element, rest_lists, current))
+            return tuple(elements)
+        else:  # rule
+            _pos, _val = item.elements
+            if rest_lists:
+                return rec_sparse(
+                    rest_lists[0],
+                    rest_lists[1:],
+                    (current[0] + _pos.elements, current[1] * _val),
+                )
+            else:
+                return (
+                    Expression(
+                        SymbolRule,
+                        ListExpression(*(current[0] + _pos.elements)),
+                        current[1] * _val,
+                    ),
+                )
+
+    # head != SparseArray
+    if not head.sameQ(SymbolSparseArray):
+        return rec(lists[0], lists[1:], [])
+
+    # head == SparseArray
+    dims = []
+    val = Integer1
+    data = []  # data = [(rules), ...]
+    for _list in lists:
+        _dims, _val, _rules = _list.elements[1:]
+        dims.extend(_dims)
+        val *= _val
+        if _val == Integer0:  # _val==0, append (_rules)
+            data.append(_rules.elements)
+        else:  # _val!=0, append (_rules, other pos->_val)
+            other_pos = []
+            for pos in itertools.product(*(range(1, d.value + 1) for d in _dims)):
+                other_pos.append(ListExpression(*(Integer(i) for i in pos)))
+            rules_pos = set(rule.elements[0] for rule in _rules.elements)
+            other_pos = set(other_pos) - rules_pos
+            other_rules = []
+            for pos in other_pos:
+                other_rules.append(Expression(SymbolRule, pos, _val))
+            data.append(_rules.elements + tuple(other_rules))
+    dims = ListExpression(*dims)
+    return Expression(
+        SymbolSparseArray,
+        SymbolAutomatic,
+        dims,
+        val,
+        ListExpression(*rec_sparse(data[0], data[1:], ((), Integer1))),
+    )
+
+
+def eval_LeviCivitaTensor(d, type):
+    "Evaluates Levi-Civita tensor of rank d"
+
+    if isinstance(d, Integer) and type == SymbolSparseArray:
+        d = d.get_int_value()
+        perms = list(permutations(list(range(1, d + 1))))
+        rules = [
+            Expression(
+                SymbolRule,
+                from_python(p),
+                from_python(Permutation.from_sequence(p).signature()),
+            )
+            for p in perms
+        ]
+        return Expression(SymbolSparseArray, from_python(rules), from_python([d] * d))