Repository: pyccel/pyccel (language: python)

Issue:
Add `cls_base` property to `TypedAstNode`
### Describe the problem 
Containers (e.g. NumPy array, CuPy array, lists) all have an underlying datatype for their elements. As a result the 4 properties in `TypedAstNode` are not sufficient to describe the type.

### Describe the solution
The `cls_base` property would complete the description. Currently this only exists in `Variable` but generalising it should help fix this problem.

### Potential issues
The definition of a `ClassDef` includes its methods which will now refer to this `ClassDef`. This could create circular dependencies. This can be avoided by declaring the `ClassDef` in `pyccel.ast.class_def` but adding the methods later in the file where they are created.

The final `class_def` must be calculated for operators. In general this should not cause too many problems as most containers are not compatible (Cuda Array + NumPy Array does work)

### Additional context
This will affect github.com/pyccel/pyccel-cuda.

cc: @bauom 


Discussion / hints:
Using `ClassDef` directly results in circular imports. I suggest having 2 properties:
- datatype
- class_type

Both these types can be declared in `ast.dtype` but a container will contain other data types

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index af0d561c8b..918b8964d4 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,7 @@ All notable changes to this project will be documented in this file.
 
 -   #1571 : Add support for the function `tuple`.
 -   #1493 : Add preliminary support for importing classes.
+-   \[INTERNALS\] Add `class_type` attribute to `TypedAstNode`.
 
 ### Fixed
 
diff --git a/pyccel/ast/basic.py b/pyccel/ast/basic.py
index 08cc0871ec..fa38daac39 100644
--- a/pyccel/ast/basic.py
+++ b/pyccel/ast/basic.py
@@ -560,6 +560,17 @@ def order(self):
         """
         return self._order # pylint: disable=no-member
 
+    @property
+    def class_type(self):
+        """
+        The type of the object.
+
+        The Python type of the object. In the case of scalars this is equivalent to
+        the datatype. For objects in (homogeneous) containers (e.g. list/ndarray/tuple),
+        this is the type of the container.
+        """
+        return self._class_type # pylint: disable=no-member
+
     @classmethod
     def static_rank(cls):
         """
@@ -618,6 +629,20 @@ class does not have a predetermined order.
         """
         return cls._order # pylint: disable=no-member
 
+    @classmethod
+    def static_class_type(cls):
+        """
+        The type of the object.
+
+        The Python type of the object. In the case of scalars this is equivalent to
+        the datatype. For objects in (homogeneous) containers (e.g. list/ndarray/tuple),
+        this is the type of the container.
+
+        This function is static and will return an AttributeError if the
+        class does not have a predetermined order.
+        """
+        return cls._class_type # pylint: disable=no-member
+
     def copy_attributes(self, x):
         """
         Copy the attributes describing a TypedAstNode into this node.
@@ -631,11 +656,12 @@ def copy_attributes(self, x):
         x : TypedAstNode
             The node from which the attributes should be copied.
         """
-        self._shape     = x.shape
-        self._rank      = x.rank
-        self._dtype     = x.dtype
-        self._precision = x.precision
-        self._order     = x.order
+        self._shape      = x.shape
+        self._rank       = x.rank
+        self._dtype      = x.dtype
+        self._precision  = x.precision
+        self._order      = x.order
+        self._class_type = x.class_type
 
 
 #------------------------------------------------------------------------------
diff --git a/pyccel/ast/bitwise_operators.py b/pyccel/ast/bitwise_operators.py
index 1cbbf3e4a1..2a867f9340 100644
--- a/pyccel/ast/bitwise_operators.py
+++ b/pyccel/ast/bitwise_operators.py
@@ -10,7 +10,7 @@
 """
 from .builtins     import PythonInt
 from .datatypes    import (NativeBool, NativeInteger, NativeFloat,
-                           NativeComplex, NativeString)
+                           NativeComplex, NativeString, NativeGeneric)
 from .internals    import max_precision
 from .operators    import PyccelUnaryOperator, PyccelOperator
 
@@ -70,7 +70,7 @@ def _calculate_dtype(self, arg):
 
         self._args      = (PythonInt(arg) if arg.dtype is NativeBool() else arg,)
         precision = arg.precision
-        return dtype, precision
+        return dtype, precision, dtype
 
     def __repr__(self):
         return f'~{repr(self.args[0])}'
@@ -93,7 +93,7 @@ class PyccelBitOperator(PyccelOperator):
     _shape = None
     _rank  = 0
     _order = None
-    __slots__ = ('_dtype','_precision')
+    __slots__ = ('_dtype','_precision','_class_type')
 
     def __init__(self, arg1, arg2):
         super().__init__(arg1, arg2)
@@ -124,17 +124,20 @@ def _calculate_dtype(self, *args):
         precision : integer
             The precision of the result of the operation.
         """
-        integers  = [a for a in args if a.dtype in (NativeInteger(),NativeBool())]
-        floats    = [a for a in args if a.dtype is NativeFloat()]
-        complexes = [a for a in args if a.dtype is NativeComplex()]
-        strs      = [a for a in args if a.dtype is NativeString()]
+        try:
+            dtype = sum((a.dtype for a in args), start=NativeGeneric())
+            class_type = sum((a.class_type for a in args), start=NativeGeneric())
+        except NotImplementedError:
+            raise TypeError(f'Cannot determine the type of {args}') #pylint: disable=raise-missing-from
 
-        if strs or complexes or floats:
+        if dtype in (NativeString(), NativeComplex(), NativeFloat()):
             raise TypeError(f'unsupported operand type(s): {args}')
-        elif integers:
-            return self._handle_integer_type(integers)
+        elif (dtype in (NativeInteger(), NativeBool())):
+            if class_type is NativeBool():
+                class_type = NativeInteger()
+            return *self._handle_integer_type(args), class_type
         else:
-            raise TypeError(f'cannot determine the type of {args}')
+            raise TypeError(f'Cannot determine the type of {args}')
 
     def _set_shape_rank(self):
         pass
diff --git a/pyccel/ast/builtins.py b/pyccel/ast/builtins.py
index 75d25e0378..804fb48497 100644
--- a/pyccel/ast/builtins.py
+++ b/pyccel/ast/builtins.py
@@ -16,8 +16,10 @@
 
 from .basic     import PyccelAstNode, TypedAstNode
 from .datatypes import (NativeInteger, NativeBool, NativeFloat,
-                        NativeComplex, NativeString, NativeGeneric)
-from .internals import PyccelInternalFunction, max_precision, Slice
+                        NativeComplex, NativeGeneric)
+from .datatypes import NativeHomogeneousTuple, NativeInhomogeneousTuple
+from .datatypes import NativeHomogeneousList
+from .internals import PyccelInternalFunction, Slice, get_final_precision
 from .literals  import LiteralInteger, LiteralFloat, LiteralComplex, Nil
 from .literals  import Literal, LiteralImaginaryUnit, convert_to_literal
 from .literals  import LiteralString
@@ -79,6 +81,7 @@ class PythonComplexProperty(PyccelInternalFunction):
     _rank  = 0
     _shape = None
     _order = None
+    _class_type = NativeFloat()
 
     def __init__(self, arg):
         super().__init__(arg)
@@ -168,6 +171,7 @@ class PythonConjugate(PyccelInternalFunction):
     _rank  = 0
     _shape = None
     _order = None
+    _class_type = NativeComplex()
     name = 'conjugate'
 
     def __new__(cls, arg):
@@ -209,6 +213,7 @@ class PythonBool(PyccelInternalFunction):
     _rank  = 0
     _shape = None
     _order = None
+    _class_type = NativeBool()
 
     def __new__(cls, arg):
         if getattr(arg, 'is_optional', None):
@@ -254,6 +259,7 @@ class PythonComplex(TypedAstNode):
     _rank  = 0
     _shape = None
     _order = None
+    _class_type = NativeComplex()
     _real_cast = PythonReal
     _imag_cast = PythonImag
     _attribute_nodes = ('_real_part', '_imag_part', '_internal_var')
@@ -416,6 +422,7 @@ class PythonFloat(PyccelInternalFunction):
     _rank  = 0
     _shape = None
     _order = None
+    _class_type = NativeFloat()
 
     def __new__(cls, arg):
         if isinstance(arg, LiteralFloat) and arg.precision == cls._precision:
@@ -460,6 +467,7 @@ class PythonInt(PyccelInternalFunction):
     _rank  = 0
     _shape = None
     _order = None
+    _class_type = NativeInteger()
 
     def __new__(cls, arg):
         if isinstance(arg, LiteralInteger):
@@ -482,11 +490,10 @@ def arg(self):
 #==============================================================================
 class PythonTuple(TypedAstNode):
     """
-    Class representing a call to Python's native tuple() function.
+    Class representing a call to Python's native (,) function which creates tuples.
 
-    Class representing a call to Python's native tuple() function
-    which either converts an object to a tuple or initialises a
-    literal tuple.
+    Class representing a call to Python's native (,) function
+    which initialises a literal tuple.
 
     Parameters
     ----------
@@ -494,7 +501,7 @@ class PythonTuple(TypedAstNode):
         The arguments passed to the tuple function.
     """
     __slots__ = ('_args','_inconsistent_shape','_is_homogeneous',
-            '_dtype','_precision','_rank','_shape','_order')
+            '_dtype','_precision','_rank','_shape','_order', '_class_type')
     _iterable        = True
     _attribute_nodes = ('_args',)
 
@@ -512,49 +519,31 @@ def __init__(self, *args):
             self._is_homogeneous = False
             return
         arg0 = args[0]
+        precision = get_final_precision(arg0)
         is_homogeneous = arg0.dtype is not NativeGeneric() and \
                          all(a.dtype is not NativeGeneric() and \
                              arg0.dtype == a.dtype and \
+                             precision == get_final_precision(a) and \
                              arg0.rank  == a.rank  and \
                              arg0.order == a.order for a in args[1:])
         self._inconsistent_shape = not all(arg0.shape==a.shape   for a in args[1:])
         self._is_homogeneous = is_homogeneous
         if is_homogeneous:
-            integers  = [a for a in args if a.dtype is NativeInteger()]
-            floats    = [a for a in args if a.dtype is NativeFloat()]
-            complexes = [a for a in args if a.dtype is NativeComplex()]
-            bools     = [a for a in args if a.dtype is NativeBool()]
-            strs      = [a for a in args if a.dtype is NativeString()]
-            if strs:
-                self._dtype = NativeString()
-                self._precision = 0
-            else:
-                if complexes:
-                    self._dtype     = NativeComplex()
-                    self._precision = max_precision(complexes)
-                elif floats:
-                    self._dtype     = NativeFloat()
-                    self._precision = max_precision(floats)
-                elif integers:
-                    self._dtype     = NativeInteger()
-                    self._precision = max_precision(integers)
-                elif bools:
-                    self._dtype     = NativeBool()
-                    self._precision  = max_precision(bools)
-                else:
-                    raise TypeError(f'cannot determine the type of {self}')
-
-
+            self._dtype = arg0.dtype
+            self._precision = arg0.precision
             inner_shape = [() if a.rank == 0 else a.shape for a in args]
             self._rank = max(a.rank for a in args) + 1
             self._shape = (LiteralInteger(len(args)), ) + inner_shape[0]
             self._rank  = len(self._shape)
 
+            self._class_type = NativeHomogeneousTuple()
+
         else:
             max_rank = max(a.rank for a in args)
-            self._rank      = max_rank + 1
-            self._dtype     = NativeGeneric()
-            self._precision = 0
+            self._rank       = max_rank + 1
+            self._dtype      = NativeInhomogeneousTuple(*[a.dtype for a in args])
+            self._precision  = 0
+            self._class_type = self._dtype
             if self._rank == 1:
                 self._shape     = (LiteralInteger(len(args)), )
             elif any(a.rank != max_rank for a in args):
@@ -617,16 +606,12 @@ def is_homogeneous(self):
 
     @property
     def args(self):
-        """ Arguments of the tuple
         """
-        return self._args
+        Arguments of the tuple.
 
-    @property
-    def allows_negative_indexes(self):
-        """ Indicates whether variables used to
-        index this Variable can be negative
+        The arguments that were used to initialise the tuple.
         """
-        return False
+        return self._args
 
 class PythonTupleFunction(TypedAstNode):
     """
@@ -676,6 +661,7 @@ class PythonLen(PyccelInternalFunction):
     _rank      = 0
     _shape     = None
     _order     = None
+    _class_type = NativeInteger()
 
     def __new__(cls, arg):
         if not getattr(arg, 'is_homogeneous', False):
@@ -699,10 +685,91 @@ def __str__(self):
         return f'len({self.arg})'
 
 #==============================================================================
-class PythonList(PythonTuple):
-    """ Represents a call to Python's native list() function.
+class PythonList(TypedAstNode):
     """
-    __slots__ = ()
+    Class representing a call to Python's `[,]` function.
+
+    Class representing a call to Python's `[,]` function which generates
+    a literal Python list.
+
+    Parameters
+    ----------
+    *args : tuple of TypedAstNodes
+        The arguments passed to the operator.
+
+    See Also
+    --------
+    FunctionalFor
+        The `[]` function when it describes a comprehension.
+    """
+    __slots__ = ('_args','_dtype','_precision','_rank','_shape','_order')
+    _attribute_nodes = ('_args',)
+    _class_type = NativeHomogeneousList()
+
+    def __init__(self, *args):
+        self._args = args
+        super().__init__()
+        if pyccel_stage == 'syntactic':
+            return
+        elif len(args) == 0:
+            self._dtype = NativeGeneric()
+            self._precision = 0
+            self._rank  = 0
+            self._shape = None
+            self._order = None
+            return
+        arg0 = args[0]
+        precision = get_final_precision(arg0)
+        is_homogeneous = arg0.dtype is not NativeGeneric() and \
+                         all(a.dtype is not NativeGeneric() and \
+                             arg0.dtype == a.dtype and \
+                             precision == get_final_precision(a) and \
+                             arg0.rank  == a.rank  and \
+                             arg0.order == a.order for a in args[1:])
+        if is_homogeneous:
+            self._dtype = arg0.dtype
+            self._precision = arg0.precision
+
+            inner_shape = [() if a.rank == 0 else a.shape for a in args]
+            self._rank = max(a.rank for a in args) + 1
+            self._shape = (LiteralInteger(len(args)), ) + inner_shape[0]
+            self._rank  = len(self._shape)
+
+        else:
+            raise TypeError("Can't create an inhomogeneous list")
+
+        self._order = None if self._rank < 2 else 'C'
+
+    def __iter__(self):
+        return self._args.__iter__()
+
+    def __str__(self):
+        args = ', '.join(str(a) for a in self)
+        return f'({args})'
+
+    def __repr__(self):
+        args = ', '.join(str(a) for a in self)
+        return f'PythonList({args})'
+
+    @property
+    def args(self):
+        """
+        Arguments of the list.
+
+        The arguments that were used to initialise the list.
+        """
+        return self._args
+
+    @property
+    def is_homogeneous(self):
+        """
+        Indicates whether the list is homogeneous or inhomogeneous.
+
+        Indicates whether all elements of the list have the same dtype, precision,
+        rank, etc (homogenous) or if these values can vary (inhomogeneous). Lists
+        are always homogeneous.
+        """
+        return True
 
 #==============================================================================
 class PythonMap(PyccelAstNode):
@@ -905,7 +972,7 @@ class PythonAbs(PyccelInternalFunction):
     x : TypedAstNode
         The argument passed to the function.
     """
-    __slots__ = ('_dtype','_precision','_rank','_shape','_order')
+    __slots__ = ('_dtype','_precision','_rank','_shape','_order','_class_type')
     name = 'abs'
     def __init__(self, x):
         self._shape     = x.shape
@@ -913,6 +980,7 @@ def __init__(self, x):
         self._dtype     = NativeInteger() if x.dtype is NativeInteger() else NativeFloat()
         self._precision = -1
         self._order     = x.order
+        self._class_type = x.class_type
         super().__init__(x)
 
     @property
@@ -936,7 +1004,7 @@ class PythonSum(PyccelInternalFunction):
     arg : TypedAstNode
         The argument passed to the function.
     """
-    __slots__ = ('_dtype','_precision')
+    __slots__ = ('_dtype','_precision','_class_type')
     name   = 'sum'
     _rank  = 0
     _shape = None
@@ -947,6 +1015,7 @@ def __init__(self, arg):
             raise TypeError(f'Unknown type of {type(arg)}.' )
         self._dtype = arg.dtype
         self._precision = -1
+        self._class_type = arg.dtype
         super().__init__(arg)
 
     @property
@@ -970,7 +1039,7 @@ class PythonMax(PyccelInternalFunction):
     *x : list, tuple, PythonTuple, PythonList
         The arguments passed to the funciton.
     """
-    __slots__ = ('_dtype','_precision')
+    __slots__ = ('_dtype','_precision','_class_type')
     name   = 'max'
     _rank  = 0
     _shape = None
@@ -990,6 +1059,7 @@ def __init__(self, *x):
                              f"types ({types}). Please cast arguments to the desired dtype")
         self._dtype     = x.dtype
         self._precision = x.precision
+        self._class_type = x.class_type
         super().__init__(x)
 
 
@@ -1005,7 +1075,7 @@ class PythonMin(PyccelInternalFunction):
     *x : list, tuple, PythonTuple, PythonList
         The arguments passed to the funciton.
     """
-    __slots__ = ('_dtype','_precision')
+    __slots__ = ('_dtype','_precision','_class_type')
     name   = 'min'
     _rank  = 0
     _shape = None
@@ -1024,6 +1094,7 @@ def __init__(self, *x):
                               f"types ({types}). Please cast arguments to the desired dtype")
         self._dtype     = x.dtype
         self._precision = x.precision
+        self._class_type = x.class_type
         super().__init__(x)
 
 #==============================================================================
diff --git a/pyccel/ast/c_concepts.py b/pyccel/ast/c_concepts.py
index c02586956a..e2ee6d5370 100644
--- a/pyccel/ast/c_concepts.py
+++ b/pyccel/ast/c_concepts.py
@@ -15,23 +15,40 @@
            'ObjectAddress')
 
 class ObjectAddress(TypedAstNode):
-    """Represents the address of an object.
-    ObjectAddress(Variable('int','a'))                            is  &a
-    ObjectAddress(Variable('int','a', memory_handling='alias'))   is   a
     """
+    Class representing the address of an object.
 
-    __slots__ = ('_obj', '_rank', '_precision', '_dtype', '_shape', '_order')
+    Class representing the address of an object. In most situations it will not be
+    necessary to use this object explicitly. E.g. if you assign a pointer to a
+    target then the pointer will be printed using `AliasAssign`. However for the
+    `_print_AliasAssign` function to print neatly, this class will be used.
+
+    Parameters
+    ----------
+    obj : TypedAstNode
+        The object whose address should be printed.
+
+    Examples
+    --------
+    >>> CCodePrinter._print(ObjectAddress(Variable('int','a')))
+    '&a'
+    >>> CCodePrinter._print(ObjectAddress(Variable('int','a', memory_handling='alias')))
+    'a'
+    """
+
+    __slots__ = ('_obj', '_rank', '_precision', '_dtype', '_shape', '_order', '_class_type')
     _attribute_nodes = ('_obj',)
 
     def __init__(self, obj):
         if not isinstance(obj, TypedAstNode):
             raise TypeError("object must be an instance of TypedAstNode")
-        self._obj       = obj
-        self._rank      = obj.rank
-        self._shape     = obj.shape
-        self._precision = obj.precision
-        self._dtype     = obj.dtype
-        self._order     = obj.order
+        self._obj        = obj
+        self._rank       = obj.rank
+        self._shape      = obj.shape
+        self._precision  = obj.precision
+        self._dtype      = obj.dtype
+        self._order      = obj.order
+        self._class_type = obj.class_type
         super().__init__()
 
     @property
diff --git a/pyccel/ast/class_defs.py b/pyccel/ast/class_defs.py
index 1f1bf317d2..f4e2da7a8c 100644
--- a/pyccel/ast/class_defs.py
+++ b/pyccel/ast/class_defs.py
@@ -8,11 +8,12 @@
 from .builtins  import PythonImag, PythonReal, PythonConjugate
 from .core      import ClassDef, PyccelFunctionDef
 from .datatypes import (NativeBool, NativeInteger, NativeFloat,
-                        NativeComplex, NativeString, NativeNumeric)
+                        NativeComplex, NativeString, NativeNumeric,
+                        NativeTuple, CustomDataType)
 from .numpyext  import (NumpyShape, NumpySum, NumpyAmin, NumpyAmax,
                         NumpyImag, NumpyReal, NumpyTranspose,
                         NumpyConjugate, NumpySize, NumpyResultType,
-                        NumpyArray)
+                        NumpyArray, NumpyNDArrayType)
 
 __all__ = ('BooleanClass',
         'IntegerClass',
@@ -179,7 +180,7 @@
 
 #=======================================================================================
 
-def get_cls_base(dtype, precision, rank):
+def get_cls_base(dtype, precision, container_type):
     """
     Determine the base class of an object.
 
@@ -193,8 +194,9 @@ def get_cls_base(dtype, precision, rank):
     precision : int
         The precision of the object.
 
-    rank : int
-        The rank of the object.
+    container_type : DataType
+        The Python type of the object. If this is different to the dtype then
+        the object is a container.
 
     Returns
     -------
@@ -206,14 +208,18 @@ def get_cls_base(dtype, precision, rank):
     NotImplementedError
         Raised if the base class cannot be found.
     """
-    if precision in (-1, 0, None) and rank == 0:
+    if isinstance(dtype, CustomDataType) and container_type is dtype:
+        return None
+    if precision in (-1, 0, None) and container_type is dtype:
         return literal_classes[dtype]
-    elif dtype in NativeNumeric:
+    elif dtype in NativeNumeric or isinstance(container_type, NumpyNDArrayType):
         return NumpyArrayClass
+    elif isinstance(container_type, NativeTuple):
+        return TupleClass
     else:
-        type_name = f"{dtype}({precision})"
-        if rank:
-            dims = ','.join(':' for _ in range(rank))
-            type_name += f"[{dims}]"
+        if container_type:
+            type_name = str(container_type)
+        else:
+            type_name = f"{dtype}({precision})"
         raise NotImplementedError(f"No class definition found for type {type_name}")
 
diff --git a/pyccel/ast/cmathext.py b/pyccel/ast/cmathext.py
index 884eeddc55..1d7585f656 100644
--- a/pyccel/ast/cmathext.py
+++ b/pyccel/ast/cmathext.py
@@ -12,7 +12,7 @@
 from pyccel.ast.datatypes import NativeBool, NativeFloat, NativeComplex
 from pyccel.ast.internals import PyccelInternalFunction
 from pyccel.ast.literals  import LiteralInteger
-from pyccel.ast.operators import PyccelOr
+from pyccel.ast.operators import PyccelAnd, PyccelOr
 from pyccel.ast.variable  import Constant
 
 from .mathext import math_constants, MathFunctionBase
@@ -52,14 +52,15 @@ class CmathFunctionBool(MathFunctionBase):
     A super-class from which functions in the `cmath` library which
     return a boolean should inherit.
 
-    Paramters
-    ---------
-    z : TypedAstNode
-        The expression passed as argument to the function.
+    Parameters
+    ----------
+    *args : TypedAstNode
+        The arguments passed to the function.
     """
     __slots__ = ()
     _dtype = NativeBool()
     _precision = -1
+    _class_type = NativeBool()
 
 class CmathFunctionComplex(MathFunctionBase):
     """
@@ -68,8 +69,8 @@ class CmathFunctionComplex(MathFunctionBase):
     A super-class from which functions in the `cmath` library which
     return a complex number should inherit.
 
-    Paramters
-    ---------
+    Parameters
+    ----------
     z : TypedAstNode
         The expression passed as argument to the function.
     """
@@ -79,6 +80,7 @@ class CmathFunctionComplex(MathFunctionBase):
     _shape = None
     _rank  = 0
     _order = None
+    _class_type = NativeComplex()
 
     def __init__(self, z : 'TypedAstNode'):
         super().__init__(z)
@@ -325,11 +327,19 @@ class CmathIsclose (CmathFunctionBool):
 
     Parameters
     ----------
-    z : TypedAstNode
-        The expression passed as argument to the function.
+    a : TypedAstNode
+        The first argument passed to the function.
+    b : TypedAstNode
+        The second argument passed to the function.
+    rel_tol : TypedAstNode
+        The relative tolerance.
+    abs_tol : TypedAstNode
+        The absolute tolerance.
     """
     __slots__ = ()
     name = 'isclose'
+    def __init__(self, a, b, *, rel_tol=1e-09, abs_tol=0.0):
+        super().__init__(a, b, rel_tol, abs_tol)
 
 #==============================================================================
 
@@ -350,7 +360,7 @@ def __new__(cls, z):
         if z.dtype is not NativeComplex():
             return MathIsfinite(z)
         else:
-            return PyccelOr(MathIsfinite(PythonImag(z)), MathIsfinite(PythonReal(z)))
+            return PyccelAnd(MathIsfinite(PythonImag(z)), MathIsfinite(PythonReal(z)))
 
 #==============================================================================
 
@@ -433,6 +443,7 @@ class CmathPolar(PyccelInternalFunction):
     _shape = (LiteralInteger(2),)
     _rank  = 1
     _order = None
+    _class_type = NativeFloat()
 
     def __init__(self, z):
         super().__init__(z)
@@ -457,6 +468,7 @@ class CmathRect(PyccelInternalFunction):
     _shape = None
     _rank  = 0
     _order = None
+    _class_type = NativeComplex()
     def __init__(self, r, phi):
         super().__init__(r, phi)
 
diff --git a/pyccel/ast/core.py b/pyccel/ast/core.py
index 643b2ec16a..4cc60708a7 100644
--- a/pyccel/ast/core.py
+++ b/pyccel/ast/core.py
@@ -15,11 +15,12 @@
 from .basic     import PyccelAstNode, TypedAstNode, iterable, ScopedAstNode
 from .builtins  import (PythonEnumerate, PythonLen, PythonMap, PythonTuple,
                         PythonRange, PythonZip, PythonBool, Lambda)
-from .datatypes import (datatype, DataType, NativeSymbol,
-                        NativeBool, NativeTuple, str_dtype)
+from .datatypes import (datatype, DataType, NativeSymbol, NativeHomogeneousTuple,
+                        NativeBool, NativeTuple, str_dtype, NativeInhomogeneousTuple,
+                        NativeVoid)
 from .internals import PyccelSymbol, PyccelInternalFunction, get_final_precision
 
-from .literals  import Nil, LiteralFalse
+from .literals  import Nil, LiteralFalse, LiteralInteger
 from .literals  import NilArgument, LiteralTrue
 
 from .operators import PyccelAdd, PyccelMinus, PyccelMul, PyccelDiv, PyccelMod
@@ -216,25 +217,32 @@ def __hash__(self):
 
 
 class Duplicate(TypedAstNode):
+    """
+    Class representing the duplicate operator in Python.
 
-    """ this is equivalent to the * operator for python tuples.
+    Class representing the duplicate operator in Python. This is equivalent
+    to the * operator for Python tuples, lists, strings, etc. In other words
+    it represents the * operator when it duplicates the first argument passed
+    to the operator, rather than acting as a numerical operator.
 
     Parameters
     ----------
-    value : TypedAstNode
-           an expression which represents the initilized value of the list
+    val : TypedAstNode
+        The object being duplicated.
 
-    shape : the shape of the array
+    length : TypedAstNode
+        The number of times the val should appear in the final object.
     """
-    __slots__ = ('_val', '_length','_dtype','_precision','_rank','_shape','_order')
+    __slots__ = ('_val', '_length','_dtype','_precision','_rank','_shape','_order','_class_type')
     _attribute_nodes = ('_val', '_length')
 
     def __init__(self, val, length):
-        self._dtype     = val.dtype
-        self._precision = val.precision
-        self._rank      = val.rank
-        self._shape     = tuple(s if i!= 0 else PyccelMul(s, length, simplify=True) for i,s in enumerate(val.shape))
-        self._order     = val.order
+        self._dtype      = val.dtype
+        self._precision  = val.precision
+        self._rank       = val.rank
+        self._shape      = tuple(s if i!= 0 else PyccelMul(s, length, simplify=True) for i,s in enumerate(val.shape))
+        self._order      = val.order
+        self._class_type = val.class_type
 
         self._val       = val
         self._length    = length
@@ -255,24 +263,29 @@ def __repr__(self):
         return '{} * {}'.format(repr(self.val), repr(self.length))
 
 class Concatenate(TypedAstNode):
+    """
+    A class representing the + operator for Python tuples.
 
-    """ this is equivalent to the + operator for python tuples
+    A class representing the + operator for Python tuples.
 
     Parameters
     ----------
-    args : TypedAstNodes
-           The tuples
+    arg1 : TypedAstNodes
+           The first tuple.
+    arg2 : TypedAstNodes
+           The second tuple.
     """
-    __slots__ = ('_args','_dtype','_precision','_rank','_shape','_order')
+    __slots__ = ('_args','_dtype','_precision','_rank','_shape','_order','_class_type')
     _attribute_nodes = ('_args',)
 
     def __init__(self, arg1, arg2):
-        self._dtype     = arg1.dtype
-        self._precision = arg1.precision
-        self._rank      = arg1.rank
-        shape_addition  = arg2.shape[0]
-        self._shape     = tuple(s if i!= 0 else PyccelAdd(s, shape_addition) for i,s in enumerate(arg1.shape))
-        self._order     = arg1.order
+        self._dtype      = arg1.dtype
+        self._precision  = arg1.precision
+        self._rank       = arg1.rank
+        shape_addition   = arg2.shape[0]
+        self._shape      = tuple(s if i!= 0 else PyccelAdd(s, shape_addition) for i,s in enumerate(arg1.shape))
+        self._order      = arg1.order
+        self._class_type = arg1.class_type
 
         self._args = (arg1, arg2)
         super().__init__()
@@ -2001,7 +2014,7 @@ class FunctionCall(TypedAstNode):
         The function where the call takes place.
     """
     __slots__ = ('_arguments','_funcdef','_interface','_func_name','_interface_name',
-                 '_dtype','_precision','_shape','_rank','_order')
+                 '_dtype','_precision','_shape','_rank','_order','_class_type')
     _attribute_nodes = ('_arguments','_funcdef','_interface')
 
     def __init__(self, func, args, current_function=None):
@@ -2070,14 +2083,37 @@ def __init__(self, func, args, current_function=None):
             if len(func.results)>0 and not isinstance(func.results[0], TypedAstNode):
                 errors.report(RECURSIVE_RESULTS_REQUIRED, symbol=func, severity="fatal")
 
-        self._funcdef       = func
-        self._arguments     = args
-        self._dtype         = func.results[0].var.dtype     if len(func.results) == 1 else NativeTuple()
-        self._rank          = func.results[0].var.rank      if len(func.results) == 1 else None
-        self._shape         = func.results[0].var.shape     if len(func.results) == 1 else None
-        self._precision     = func.results[0].var.precision if len(func.results) == 1 else None
-        self._order         = func.results[0].var.order     if len(func.results) == 1 else None
-        self._func_name     = func.name
+        self._funcdef    = func
+        self._arguments  = args
+        self._func_name  = func.name
+        n_results = len(func.results)
+        if n_results == 1:
+            self._dtype      = func.results[0].var.dtype
+            self._rank       = func.results[0].var.rank
+            self._shape      = func.results[0].var.shape
+            self._precision  = func.results[0].var.precision
+            self._order      = func.results[0].var.order
+            self._class_type = func.results[0].var.class_type
+        elif n_results == 0:
+            self._dtype      = NativeVoid()
+            self._rank       = 0
+            self._shape      = None
+            self._precision  = None
+            self._order      = None
+            self._class_type = NativeVoid()
+        else:
+            dtypes = [r.var.dtype for r in func.results]
+            if all(d is dtypes[0] for d in dtypes):
+                dtype = NativeHomogeneousTuple()
+            else:
+                dtype = NativeInhomogeneousTuple(*dtypes)
+            self._dtype      = dtype
+            self._rank       = 1
+            self._shape      = (LiteralInteger(n_results),)
+            self._precision  = None
+            self._order      = None
+            self._class_type = dtype
+
         super().__init__()
 
     @property
@@ -4257,14 +4293,19 @@ def args_var(self):
 
 class InProgram(TypedAstNode):
     """
-    Class representing the boolean:
-    __name__ == '__main__'
+    Class representing the 'in program' test.
+
+    Class representing the test indicating whether the code should
+    only be executed when the file is executed as a program. In
+    other words, a class representing the boolean:
+    `__name__ == '__main__'`
     """
     _dtype = NativeBool()
     _precision = -1
     _rank  = 0
     _shape = None
     _order = None
+    _class_type = NativeBool()
     _attribute_nodes = ()
     __slots__ = ()
 
diff --git a/pyccel/ast/cwrapper.py b/pyccel/ast/cwrapper.py
index e9dc827bff..9b9e607941 100644
--- a/pyccel/ast/cwrapper.py
+++ b/pyccel/ast/cwrapper.py
@@ -229,6 +229,7 @@ class PyBuildValueNode(TypedAstNode):
     _rank = 0
     _precision = 0
     _shape = ()
+    _class_type = PyccelPyObject
     _order = None
 
     def __init__(self, result_args = ()):
@@ -270,6 +271,7 @@ class PyModule_AddObject(TypedAstNode):
     _precision = 4
     _rank = 0
     _shape = None
+    _class_type = NativeInteger()
 
     def __init__(self, mod_name, name, variable):
         if not isinstance(name, LiteralString):
diff --git a/pyccel/ast/datatypes.py b/pyccel/ast/datatypes.py
index 1da00acf13..4fa284b0ff 100644
--- a/pyccel/ast/datatypes.py
+++ b/pyccel/ast/datatypes.py
@@ -8,10 +8,11 @@
 """
 Classes and methods that handle supported datatypes in C/Fortran.
 """
+from functools import lru_cache
 
 import numpy
 
-from pyccel.utilities.metaclasses import Singleton
+from pyccel.utilities.metaclasses import ArgumentSingleton, Singleton
 
 # TODO [YG, 12.03.2020] verify why we need all these types
 # NOTE: symbols not used in pyccel are commented out
@@ -25,12 +26,13 @@
     'NativeComplex',
     'NativeFloat',
     'NativeGeneric',
+    'NativeInhomogeneousTuple',
     'NativeInteger',
-    'NativeTuple',
-    'NativeNil',
-    'NativeRange',
+    'NativeHomogeneousList',
+    'NativeHomogeneousTuple',
     'NativeString',
     'NativeSymbol',
+    'NativeTuple',
     'NativeVoid',
     'UnionType',
     'DataTypeFactory',
@@ -46,7 +48,6 @@
     'Cmplx',
     'Generic',
     'Int',
-    'Nil',
     'Float',
     'String',
     'Void',
@@ -92,11 +93,19 @@
 
 #==============================================================================
 
-class DataType(metaclass=Singleton):
+class DataType(metaclass=ArgumentSingleton):
     """
     Base class representing native datatypes.
 
     The base class from which all data types must inherit.
+
+    Parameters
+    ----------
+    *args : tuple
+        Any arguments required by the class.
+
+    **kwargs : dict
+        Any keyword arguments required by the class.
     """
     __slots__ = ()
     _name = '__UNDEFINED__'
@@ -132,55 +141,208 @@ def __reduce__(self):
         """
         return (self.__class__, ())
 
-class NativeBool(DataType):
-    """Class representing boolean datatype"""
+class NativeBool(DataType, metaclass=Singleton):
+    """
+    Class representing a boolean datatype.
+
+    Class representing a boolean datatype.
+    """
     __slots__ = ()
     _name = 'Bool'
 
-class NativeInteger(DataType):
-    """Class representing integer datatype"""
+    @lru_cache
+    def __add__(self, other):
+        if other in NativeNumeric:
+            return other
+        else:
+            return NotImplemented
+
+class NativeInteger(DataType, metaclass=Singleton):
+    """
+    Class representing an integer datatype.
+
+    Class representing an integer datatype.
+    """
     __slots__ = ()
     _name = 'Int'
 
-class NativeFloat(DataType):
-    """Class representing float datatype"""
+    @lru_cache
+    def __add__(self, other):
+        if other in NativeNumeric:
+            if other is NativeBool():
+                return self
+            else:
+                return other
+        else:
+            return NotImplemented
+
+class NativeFloat(DataType, metaclass=Singleton):
+    """
+    Class representing a float datatype.
+
+    Class representing a float datatype.
+    """
     __slots__ = ()
     _name = 'Float'
 
-class NativeComplex(DataType):
-    """Class representing complex datatype"""
+    @lru_cache
+    def __add__(self, other):
+        if other in NativeNumeric:
+            if other is NativeComplex():
+                return other
+            else:
+                return self
+        else:
+            return NotImplemented
+
+class NativeComplex(DataType, metaclass=Singleton):
+    """
+    Class representing a complex datatype.
+
+    Class representing a complex datatype.
+    """
     __slots__ = ()
     _name = 'Complex'
 
+    @lru_cache
+    def __add__(self, other):
+        if other in NativeNumeric:
+            return self
+        else:
+            return NotImplemented
+
 NativeNumeric = (NativeBool(), NativeInteger(), NativeFloat(), NativeComplex())
+NativeNumericTypes = (NativeBool, NativeInteger, NativeFloat, NativeComplex)
 
-class NativeString(DataType):
-    """Class representing string datatype"""
+class NativeString(DataType, metaclass=Singleton):
+    """
+    Class representing a string datatype.
+
+    Class representing a string datatype.
+    """
     __slots__ = ()
     _name = 'String'
 
-class NativeVoid(DataType):
-    __slots__ = ()
-    _name = 'Void'
+    @lru_cache
+    def __add__(self, other):
+        if isinstance(other, NativeString):
+            return self
+        else:
+            return NotImplemented
+
+class NativeVoid(DataType, metaclass=Singleton):
+    """
+    Class representing a void datatype.
 
-class NativeNil(DataType):
+    Class representing a void datatype. This class is especially useful
+    in the C-Python wrapper when a `void*` type is needed to collect
+    pointers from Fortran.
+    """
     __slots__ = ()
-    _name = 'Nil'
+    _name = 'Void'
 
 class NativeTuple(DataType):
-    """Base class representing native datatypes"""
+    """
+    Base class representing tuple datatypes.
+
+    The class from which tuple datatypes must inherit.
+
+    Parameters
+    ----------
+    *args : tuple
+        Any arguments required by the class.
+
+    **kwargs : dict
+        Any keyword arguments required by the class.
+    """
     __slots__ = ()
     _name = 'Tuple'
 
-class NativeRange(DataType):
+    @lru_cache
+    def __add__(self, other):
+        if isinstance(other, NativeTuple):
+            return self
+        else:
+            return NotImplemented
+
+class NativeHomogeneousTuple(NativeTuple, metaclass = Singleton):
+    """
+    Class representing the homogeneous tuple type.
+
+    Class representing the type of a homogeneous tuple. This
+    is a container type and should be used as the class_type.
+    """
     __slots__ = ()
-    _name = 'Range'
 
-class NativeSymbol(DataType):
+class NativeInhomogeneousTuple(NativeTuple):
+    """
+    Class representing the inhomogeneous tuple type.
+
+    Class representing the type of an inhomogeneous tuple. This is a
+    basic datatype as it cannot be arbitrarily indexed. It is
+    therefore parametrised by the datatypes that it contains.
+
+    Parameters
+    ----------
+    *args : tuple of DataTypes
+        The datatypes stored in the inhomogeneous tuple.
+
+    **kwargs : empty dict
+        Keyword arguments as defined by the ArgumentSingleton class.
+    """
+    __slots__ = ('_dtypes',)
+
+    def __init__(self, *args):
+        self._dtypes = args
+        super().__init__()
+
+    @property
+    def name(self):
+        """
+        The name of the datatype.
+
+        Get the name of the datatype. This name is parametrised by the
+        datatypes in the elements of the tuple.
+
+        Returns
+        -------
+        str
+            The name of the datatype.
+        """
+        datatypes = ', '.join(d.name for d in self._dtypes)
+        return f'tuple[{datatypes}]'
+
+    def __getitem__(self, i):
+        return self._dtypes[i]
+
+class NativeHomogeneousList(DataType, metaclass = Singleton):
+    """
+    Class representing the homogeneous list type.
+
+    Class representing the type of a homogeneous list. This
+    is a container type and should be used as the class_type.
+    """
+    __slots__ = ()
+    _name = 'List'
+
+    @lru_cache
+    def __add__(self, other):
+        if isinstance(other, NativeHomogeneousList):
+            return self
+        else:
+            return NotImplemented
+
+class NativeSymbol(DataType, metaclass=Singleton):
+    """
+    Class representing a symbol datatype.
+
+    Class representing a symbol datatype. This **may** be useful for
+    the sympy decorator and other symbolic manipulations.
+    """
     __slots__ = ()
     _name = 'Symbol'
 
-class CustomDataType(DataType):
+class CustomDataType(DataType, metaclass=Singleton):
     """
     Class from which user-defined types inherit.
 
@@ -189,10 +351,21 @@ class CustomDataType(DataType):
     """
     __slots__ = ()
 
-class NativeGeneric(DataType):
+class NativeGeneric(DataType, metaclass = Singleton):
+    """
+    Class representing a generic datatype.
+
+    Class representing a generic datatype. This datatype is
+    useful for describing the type of an empty container (list/tuple/etc)
+    or an argument which can accept any type (e.g. MPI arguments).
+    """
     __slots__ = ()
     _name = 'Generic'
 
+    @lru_cache
+    def __add__(self, other):
+        return other
+
 # ...
 
 
@@ -202,7 +375,6 @@ class NativeGeneric(DataType):
 Float          = NativeFloat()
 Cmplx          = NativeComplex()
 Void           = NativeVoid()
-Nil            = NativeNil()
 String         = NativeString()
 _Symbol        = NativeSymbol()
 Generic        = NativeGeneric()
@@ -232,7 +404,6 @@ class NativeGeneric(DataType):
                                  'bool' :(Bool,-1),
                                  'b1' :(Bool,-1),
                                  'void' : (Void, 0),
-                                 'nil' : (Nil, 0),
                                  'symbol' : (_Symbol, 0),
                                  '*' : (Generic, 0),
                                  'str' : (String, 0),
diff --git a/pyccel/ast/functionalexpr.py b/pyccel/ast/functionalexpr.py
index b767bc0716..b574305fb5 100644
--- a/pyccel/ast/functionalexpr.py
+++ b/pyccel/ast/functionalexpr.py
@@ -23,21 +23,23 @@
 class FunctionalFor(TypedAstNode):
 
     """
+    Represents a generator expression.
+
     Represents any generator expression e.g:
     a = [i for i in range(10)]
 
     Parameters
     ----------
-    loops   : CodeBlock/For
-              The loops contained in the expression
-    expr    : PyccelAstNode
+    loops : CodeBlock/For
+              The loops contained in the expression.
+    expr : PyccelAstNode
               The expression at the origin of the expression
-              E.g. 'i' for '[i for i in range(10)]'
-    lhs     : Variable
-              The variable to which the result is assigned
+              E.g. 'i' for '[i for i in range(10)]'.
+    lhs : Variable
+              The variable to which the result is assigned.
     indices : list of Variable
-              All iterator targets for the for loops
-    index   : Variable
+              All iterator targets for the for loops.
+    index : Variable
               Index of result in rhs
               E.g.:
               ```
@@ -50,10 +52,10 @@ class FunctionalFor(TypedAstNode):
                   a[Dummy_0]=i
                   Dummy_0 += 1
               ```
-              Index is `Dummy_0`
+              Index is `Dummy_0`.
     """
     __slots__ = ('_loops','_expr', '_lhs', '_indices','_index',
-            '_dtype','_precision','_rank','_shape','_order')
+            '_dtype','_precision','_rank','_shape','_order','_class_type')
     _attribute_nodes = ('_loops','_expr', '_lhs', '_indices','_index')
 
     def __init__(
@@ -72,11 +74,12 @@ def __init__(
         super().__init__()
 
         if pyccel_stage != 'syntactic':
-            self._dtype     = lhs.dtype
-            self._precision = lhs.precision
-            self._rank      = lhs.rank
-            self._shape     = lhs.shape
-            self._order     = lhs.order
+            self._dtype      = lhs.dtype
+            self._precision  = lhs.precision
+            self._rank       = lhs.rank
+            self._shape      = lhs.shape
+            self._order      = lhs.order
+            self._class_type = lhs.class_type
 
     @property
     def loops(self):
diff --git a/pyccel/ast/internals.py b/pyccel/ast/internals.py
index 3ed4292c0d..f674501ba0 100644
--- a/pyccel/ast/internals.py
+++ b/pyccel/ast/internals.py
@@ -90,6 +90,7 @@ class PyccelArraySize(PyccelInternalFunction):
     _rank  = 0
     _shape = None
     _order = None
+    _class_type = NativeInteger()
 
     def __init__(self, arg):
         super().__init__(arg)
@@ -137,6 +138,7 @@ class PyccelArrayShapeElement(PyccelInternalFunction):
     _rank  = 0
     _shape = None
     _order = None
+    _class_type = NativeInteger()
 
     def __init__(self, arg, index):
         if not isinstance(arg, TypedAstNode):
diff --git a/pyccel/ast/literals.py b/pyccel/ast/literals.py
index a3e556595d..ddc6aeb043 100644
--- a/pyccel/ast/literals.py
+++ b/pyccel/ast/literals.py
@@ -98,8 +98,9 @@ class LiteralTrue(LiteralNumeric):
     precision : int
         The precision of the data type.
     """
-    __slots__ = ()
-    _dtype     = NativeBool()
+    __slots__   = ()
+    _dtype      = NativeBool()
+    _class_type = NativeBool()
 
     def __init__(self, precision = -1):
         super().__init__(precision)
@@ -125,8 +126,9 @@ class LiteralFalse(LiteralNumeric):
     precision : int
         The precision of the data type.
     """
-    __slots__ = ()
-    _dtype     = NativeBool()
+    __slots__   = ()
+    _dtype      = NativeBool()
+    _class_type = NativeBool()
 
     def __init__(self, precision = -1):
         super().__init__(precision)
@@ -155,8 +157,9 @@ class LiteralInteger(LiteralNumeric):
     precision : int
         The precision of the integer.
     """
-    __slots__ = ('_value',)
-    _dtype     = NativeInteger()
+    __slots__   = ('_value',)
+    _dtype      = NativeInteger()
+    _class_type = NativeInteger()
 
     def __init__(self, value, precision = -1):
         super().__init__(precision)
@@ -192,8 +195,9 @@ class LiteralFloat(LiteralNumeric):
     precision : int
         The precision of the float.
     """
-    __slots__ = ('_value',)
-    _dtype     = NativeFloat()
+    __slots__   = ('_value',)
+    _dtype      = NativeFloat()
+    _class_type = NativeFloat()
 
     def __init__(self, value, *, precision = -1):
         if not isinstance(value, (int, float, LiteralFloat)):
@@ -232,8 +236,9 @@ class LiteralComplex(LiteralNumeric):
     precision : int
         The precision of the complex number.
     """
-    __slots__ = ('_real_part','_imag_part')
-    _dtype     = NativeComplex()
+    __slots__   = ('_real_part','_imag_part')
+    _dtype      = NativeComplex()
+    _class_type = NativeComplex()
 
     def __new__(cls, real, imag, precision = -1):
         if cls is LiteralImaginaryUnit:
@@ -322,8 +327,9 @@ class LiteralString(Literal):
         The Python literal.
     """
     __slots__ = ('_string',)
-    _dtype     = NativeString()
-    _precision = 0
+    _dtype      = NativeString()
+    _class_type = NativeString()
+    _precision  = 0
 
     def __init__(self, arg):
         super().__init__()
@@ -361,8 +367,9 @@ class Nil(Literal, metaclass=Singleton):
     """
     __slots__ = ()
     _attribute_nodes = ()
-    _dtype = NativeGeneric
+    _dtype = NativeGeneric()
     _precision = 0
+    _class_type = NativeGeneric()
 
     def __str__(self):
         return 'None'
diff --git a/pyccel/ast/macros.py b/pyccel/ast/macros.py
index 9b4bd3dee4..e360d3b0d3 100644
--- a/pyccel/ast/macros.py
+++ b/pyccel/ast/macros.py
@@ -26,14 +26,23 @@
 
 #==============================================================================
 class Macro(TypedAstNode):
-    """."""
+    """
+    Class representing a macro, ie an inline definition.
+
+    Class representing a macro, ie an inline definition. TO BE DEPRECATED.
+
+    Parameters
+    ----------
+    argument : PyccelSymbol | Variable
+        The argument passed to the macro.
+    """
     __slots__ = ('_argument',)
     _name = '__UNDEFINED__'
     _attribute_nodes = ()
 
     def __init__(self, argument):
         if not isinstance(argument, (PyccelSymbol, Variable)):
-            raise TypeError("Argument must be a Pyccelsymbol or a Variable not {}".format(type(argument)))
+            raise TypeError(f"Argument must be a Pyccelsymbol or a Variable not {type(argument)}")
 
         self._argument = argument
         super().__init__()
@@ -48,12 +57,24 @@ def name(self):
 
 #==============================================================================
 class MacroShape(Macro):
-    """."""
+    """
+    A macro indicating the shape of a variable.
+
+    A macro indicating the shape of a variable. TO BE DEPRECATED.
+
+    Parameters
+    ----------
+    argument : PyccelSymbol | Variable
+        The variable whose shape we are accessing.
+    index : int | LiteralInteger
+        The index of the element of the shape we are accessing.
+    """
     __slots__ = ('_index','_rank','_shape')
     _name      = 'shape'
     _dtype     = NativeInteger()
     _precision = -1
     _order     = None
+    _class_type = NativeInteger()
 
     def __init__(self, argument, index=None):
         if index is not None:
@@ -74,14 +95,22 @@ def index(self):
 
     def __str__(self):
         if self.index is None:
-            return 'MacroShape({})'.format(str(self.argument))
+            return f'MacroShape({self.argument})'
         else:
-            return 'MacroShape({}, {})'.format(str(self.argument),
-                                               str(self.index))
+            return f'MacroShape({self.argument}, {self.index})'
 
 #==============================================================================
 class MacroType(Macro):
-    """."""
+    """
+    A macro representing the type of a variable.
+
+    A macro representing the type of a variable. TO BE DEPRECATED.
+
+    Parameters
+    ----------
+    argument : PyccelSymbol | Variable
+        The variable whose datatype we are accessing.
+    """
     __slots__ = ()
     _name      = 'dtype'
     _dtype     = NativeGeneric()
@@ -89,13 +118,24 @@ class MacroType(Macro):
     _rank      = 0
     _shape     = None
     _order     = None
+    _class_type = NativeGeneric()
 
     def __str__(self):
-        return 'MacroType({})'.format(str(self.argument))
+        return f'MacroType({self.argument})'
 
 #==============================================================================
 class MacroCount(Macro):
-    """."""
+    """
+    A macro representing the total number of elements in a variable.
+
+    A macro representing the total number of elements in a variable.
+    TO BE DEPRECATED.
+
+    Parameters
+    ----------
+    argument : PyccelSymbol | Variable
+        The variable whose size we are accessing.
+    """
     __slots__ = ()
     _name      = 'count'
     _dtype     = NativeInteger()
@@ -103,9 +143,10 @@ class MacroCount(Macro):
     _rank      = 0
     _shape     = None
     _order     = None
+    _class_type = NativeInteger()
 
     def __str__(self):
-        return 'MacroCount({})'.format(str(self.argument))
+        return f'MacroCount({self.argument})'
 
 
 
diff --git a/pyccel/ast/mathext.py b/pyccel/ast/mathext.py
index 65d321a7cb..e9b4446583 100644
--- a/pyccel/ast/mathext.py
+++ b/pyccel/ast/mathext.py
@@ -76,7 +76,17 @@
 # Base classes
 #==============================================================================
 class MathFunctionBase(PyccelInternalFunction):
-    """Abstract base class for the Math Functions"""
+    """
+    Abstract base class for the Math Functions.
+
+    A super-class from which all functions in the `math` library
+    should inherit.
+
+    Parameters
+    ----------
+    *args : tuple of TypedAstNode
+        The arguments passed to the function.
+    """
     __slots__ = ()
     _shape = None
     _rank  = 0
@@ -84,24 +94,60 @@ class MathFunctionBase(PyccelInternalFunction):
 
 
 class MathFunctionFloat(MathFunctionBase):
+    """
+    Super-class from which functions returning a float inherit.
+
+    A super-class from which functions in the `math` library which
+    return a float should inherit.
+
+    Parameters
+    ----------
+    *args : tuple of TypedAstNode
+        The arguments passed to the function.
+    """
     __slots__ = ()
     name = 'float'
     _dtype = NativeFloat()
     _precision = -1
+    _class_type = NativeFloat()
 
 
 class MathFunctionInt(MathFunctionBase):
+    """
+    Super-class from which functions returning an integer inherit.
+
+    A super-class from which functions in the `math` library which
+    return an integer should inherit.
+
+    Parameters
+    ----------
+    *args : tuple of TypedAstNode
+        The arguments passed to the function.
+    """
     __slots__ = ()
     name = 'int'
     _dtype = NativeInteger()
     _precision = -1
+    _class_type = NativeInteger()
 
 
 class MathFunctionBool(MathFunctionBase):
+    """
+    Super-class from which functions returning a boolean inherit.
+
+    A super-class from which functions in the `math` library which
+    return a boolean should inherit.
+
+    Parameters
+    ----------
+    *args : tuple of TypedAstNode
+        The arguments passed to the function.
+    """
     __slots__ = ()
     name = 'bool'
     _dtype = NativeBool()
     _precision = -1
+    _class_type = NativeBool()
 
 #==============================================================================
 # Functions that return one value
@@ -121,6 +167,8 @@ class MathAcos(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'acos'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathAcosh(MathFunctionFloat):
@@ -136,6 +184,8 @@ class MathAcosh(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'acosh'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathAsin(MathFunctionFloat):
@@ -151,6 +201,8 @@ class MathAsin(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'asin'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathAsinh(MathFunctionFloat):
@@ -166,6 +218,8 @@ class MathAsinh(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'asinh'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathAtan(MathFunctionFloat):
@@ -181,6 +235,8 @@ class MathAtan(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'atan'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathAtan2   (MathFunctionFloat):
@@ -198,7 +254,8 @@ class MathAtan2   (MathFunctionFloat):
     """
     __slots__ = ()
     name = 'atan2'
-
+    def __init__(self, y, x):
+        super().__init__(y, x)
 
 class MathAtanh(MathFunctionFloat):
     """
@@ -213,6 +270,8 @@ class MathAtanh(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'atanh'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathCopysign(MathFunctionFloat):
@@ -224,10 +283,14 @@ class MathCopysign(MathFunctionFloat):
     Parameters
     ----------
     x : TypedAstNode
-        The expression passed as argument to the function.
+        The first expression passed as argument to the function.
+    y : TypedAstNode
+        The second expression passed as argument to the function.
     """
     __slots__ = ()
     name = 'copysign'
+    def __init__(self, x, y):
+        super().__init__(x, y)
 
 
 class MathCos(MathFunctionFloat):
@@ -243,6 +306,8 @@ class MathCos(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'cos'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathCosh(MathFunctionFloat):
@@ -258,6 +323,8 @@ class MathCosh(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'cosh'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathErf(MathFunctionFloat):
@@ -273,6 +340,8 @@ class MathErf(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'erf'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathErfc(MathFunctionFloat):
@@ -288,6 +357,8 @@ class MathErfc(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'erfc'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathExp(MathFunctionFloat):
@@ -303,12 +374,25 @@ class MathExp(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'exp'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathExpm1   (MathFunctionFloat):
-    """Represent a call to the expm1 function in the Math library"""
+    """
+    Class representing a call to the `math.expm1` function.
+
+    A class which represents a call to the `expm1` function from the `math` library.
+
+    Parameters
+    ----------
+    x : TypedAstNode
+        The expression passed as argument to the function.
+    """
     __slots__ = ()
     name = 'expm1'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathFabs(MathFunctionFloat):
@@ -324,6 +408,8 @@ class MathFabs(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'fabs'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathFmod(MathFunctionFloat):
@@ -339,6 +425,8 @@ class MathFmod(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'fmod'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathFsum(MathFunctionFloat):
@@ -354,6 +442,8 @@ class MathFsum(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'fsum'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathGamma(MathFunctionFloat):
@@ -369,6 +459,8 @@ class MathGamma(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'gamma'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathHypot(MathFunctionFloat):
@@ -379,8 +471,8 @@ class MathHypot(MathFunctionFloat):
 
     Parameters
     ----------
-    x : TypedAstNode
-        The expression passed as argument to the function.
+    *args : TypedAstNode
+        The arguments passed to the function.
     """
     __slots__ = ()
     name = 'hypot'
@@ -395,10 +487,14 @@ class MathLdexp(MathFunctionFloat):
     Parameters
     ----------
     x : TypedAstNode
-        The expression passed as argument to the function.
+        The first expression passed as argument to the function.
+    i : TypedAstNode
+        The second expression passed as argument to the function.
     """
     __slots__ = ()
     name = 'ldexp'
+    def __init__(self, x, i):
+        super().__init__(x, i)
 
 
 class MathLgamma(MathFunctionFloat):
@@ -414,6 +510,8 @@ class MathLgamma(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'lgamma'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathLog(MathFunctionFloat):
@@ -429,6 +527,8 @@ class MathLog(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'log'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathLog10   (MathFunctionFloat):
@@ -444,6 +544,8 @@ class MathLog10   (MathFunctionFloat):
     """
     __slots__ = ()
     name = 'log10'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathLog1p   (MathFunctionFloat):
@@ -459,6 +561,8 @@ class MathLog1p   (MathFunctionFloat):
     """
     __slots__ = ()
     name = 'log1p'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathLog2    (MathFunctionFloat):
@@ -474,6 +578,8 @@ class MathLog2    (MathFunctionFloat):
     """
     __slots__ = ()
     name = 'log2'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathPow(MathFunctionFloat):
@@ -485,10 +591,14 @@ class MathPow(MathFunctionFloat):
     Parameters
     ----------
     x : TypedAstNode
-        The expression passed as argument to the function.
+        The first expression passed as argument to the function.
+    y : TypedAstNode
+        The second expression passed as argument to the function.
     """
     __slots__ = ()
     name = 'pow'
+    def __init__(self, x, y):
+        super().__init__(x, y)
 
 
 class MathSin(MathFunctionFloat):
@@ -504,6 +614,8 @@ class MathSin(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'sin'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathSinh(MathFunctionFloat):
@@ -519,6 +631,8 @@ class MathSinh(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'sinh'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathSqrt(MathFunctionFloat):
@@ -534,6 +648,8 @@ class MathSqrt(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'sqrt'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathTan(MathFunctionFloat):
@@ -549,6 +665,8 @@ class MathTan(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'tan'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathTanh(MathFunctionFloat):
@@ -564,6 +682,8 @@ class MathTanh(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'tanh'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathRemainder(MathFunctionFloat):
@@ -575,10 +695,14 @@ class MathRemainder(MathFunctionFloat):
     Parameters
     ----------
     x : TypedAstNode
-        The expression passed as argument to the function.
+        The first expression passed as argument to the function.
+    y : TypedAstNode
+        The second expression passed as argument to the function.
     """
     __slots__ = ()
     name = 'remainder'
+    def __init__(self, x, y):
+        super().__init__(x, y)
 
 class MathRadians(MathFunctionFloat):
     """
@@ -593,6 +717,8 @@ class MathRadians(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'radians'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathDegrees(MathFunctionFloat):
@@ -608,6 +734,8 @@ class MathDegrees(MathFunctionFloat):
     """
     __slots__ = ()
     name = 'degrees'
+    def __init__(self, x):
+        super().__init__(x)
 
 # Integer result
 class MathFactorial(MathFunctionInt):
@@ -623,6 +751,8 @@ class MathFactorial(MathFunctionInt):
     """
     __slots__ = ()
     name = 'factorial'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathGcd(MathFunctionInt):
@@ -633,8 +763,8 @@ class MathGcd(MathFunctionInt):
 
     Parameters
     ----------
-    x : TypedAstNode
-        The expression passed as argument to the function.
+    *args : tuple of TypedAstNode
+        The arguments passed to the function.
     """
     __slots__ = ()
     name = 'gcd'
@@ -648,8 +778,8 @@ class MathLcm(MathFunctionInt):
 
     Parameters
     ----------
-    x : TypedAstNode
-        The expression passed as argument to the function.
+    *args : tuple of TypedAstNode
+        The arguments passed to the function.
     """
     __slots__ = ()
     name = 'lcm'
@@ -667,6 +797,8 @@ class MathCeil(MathFunctionInt):
     """
     __slots__ = ()
     name = 'ceil'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathFloor(MathFunctionInt):
@@ -682,6 +814,8 @@ class MathFloor(MathFunctionInt):
     """
     __slots__ = ()
     name = 'floor'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathTrunc(MathFunctionInt):
@@ -697,6 +831,8 @@ class MathTrunc(MathFunctionInt):
     """
     __slots__ = ()
     name = 'trunc'
+    def __init__(self, x):
+        super().__init__(x)
 
 # Boolean result
 class MathIsclose(MathFunctionBool):
@@ -707,11 +843,15 @@ class MathIsclose(MathFunctionBool):
 
     Parameters
     ----------
-    x : TypedAstNode
-        The expression passed as argument to the function.
+    a : TypedAstNode
+        The first expression passed as argument to the function.
+    b : TypedAstNode
+        The second expression passed as argument to the function.
     """
     __slots__ = ()
     name = 'isclose'
+    def __init__(self, a, b):
+        super().__init__(a, b)
 
 
 class MathIsfinite(MathFunctionBool):
@@ -727,6 +867,8 @@ class MathIsfinite(MathFunctionBool):
     """
     __slots__ = ()
     name = 'isfinite'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathIsinf(MathFunctionBool):
@@ -742,6 +884,8 @@ class MathIsinf(MathFunctionBool):
     """
     __slots__ = ()
     name = 'isinf'
+    def __init__(self, x):
+        super().__init__(x)
 
 
 class MathIsnan(MathFunctionBool):
@@ -757,6 +901,8 @@ class MathIsnan(MathFunctionBool):
     """
     __slots__ = ()
     name = 'isnan'
+    def __init__(self, x):
+        super().__init__(x)
 
 #==============================================================================
 # Functions that return two values
diff --git a/pyccel/ast/numpy_wrapper.py b/pyccel/ast/numpy_wrapper.py
index 1edf438641..6656aff7c2 100644
--- a/pyccel/ast/numpy_wrapper.py
+++ b/pyccel/ast/numpy_wrapper.py
@@ -123,7 +123,7 @@ def get_numpy_max_acceptable_version_file():
 array_get_data  = FunctionDef(name   = 'nd_data',
                            body      = [],
                            arguments = [FunctionDefArgument(Variable(dtype=NativeVoid(), name = 'o', is_optional=True))],
-                           results   = [FunctionDefResult(Variable(dtype=NativeVoid(), name = 'v', memory_handling='alias', rank = 1))])
+                           results   = [FunctionDefResult(Variable(dtype=NativeVoid(), name = 'v', memory_handling='alias', rank = 1, class_type = NativeVoid()))])
 
 # Basic Array Flags
 # https://numpy.org/doc/stable/reference/c-api/array.html#c.NPY_ARRAY_OWNDATA
diff --git a/pyccel/ast/numpyext.py b/pyccel/ast/numpyext.py
index 6b5c71483d..59874e8761 100644
--- a/pyccel/ast/numpyext.py
+++ b/pyccel/ast/numpyext.py
@@ -6,12 +6,14 @@
 #------------------------------------------------------------------------------------------#
 """ Module containing objects from the numpy module understood by pyccel
 """
+from functools import lru_cache
 
 import numpy
 
 from pyccel.errors.errors import Errors
 from pyccel.errors.messages import WRONG_LINSPACE_ENDPOINT, NON_LITERAL_KEEP_DIMS, NON_LITERAL_AXIS
 
+from pyccel.utilities.metaclasses import Singleton
 from pyccel.utilities.stage import PyccelStage
 
 from .basic          import TypedAstNode
@@ -22,7 +24,7 @@
 from .core           import Module, Import, PyccelFunctionDef, FunctionCall
 
 from .datatypes      import (dtype_and_precision_registry as dtype_registry,
-                             default_precision, NativeInteger,
+                             default_precision, NativeInteger, DataType, NativeNumericTypes,
                              NativeFloat, NativeComplex, NativeBool, NativeNumeric)
 
 from .internals      import PyccelInternalFunction, Slice, get_final_precision
@@ -95,6 +97,7 @@
     'NumpyInt64',
     'NumpyLinspace',
     'NumpyMatmul',
+    'NumpyNDArrayType',
     'NumpyNewArray',
     'NumpyMod',
     'NumpyNonZero',
@@ -115,6 +118,26 @@
     'NumpyShape',
 )
 
+class NumpyNDArrayType(DataType, metaclass=Singleton):
+    """
+    Class representing the NumPy ND array type.
+
+    Class representing the NumPy ND array type.
+    """
+    __slots__ = ()
+    name = 'numpy.ndarray'
+
+    @lru_cache
+    def __add__(self, other):
+        if isinstance(other, (*NativeNumericTypes, NumpyNDArrayType)):
+            return self
+        else:
+            raise NotImplementedError(f"Can't combine {self} and {other}")
+
+    @lru_cache
+    def __radd__(self, other):
+        return self.__add__(other)
+
 #=======================================================================================
 def process_shape(is_scalar, shape):
     """ Modify the input shape to the expected type
@@ -155,12 +178,13 @@ class NumpyFloat(PythonFloat):
     arg : TypedAstNode
         The argument passed to the function.
     """
-    __slots__ = ('_rank','_shape','_order')
+    __slots__ = ('_rank','_shape','_order','_class_type')
     name = 'float'
     def __init__(self, arg):
         self._shape = arg.shape
         self._rank  = arg.rank
         self._order = arg.order
+        self._class_type = arg.class_type
         super().__init__(arg)
 
     @property
@@ -199,12 +223,13 @@ class NumpyBool(PythonBool):
     arg : TypedAstNode
         The argument passed to the function.
     """
-    __slots__ = ('_shape','_rank','_order')
+    __slots__ = ('_shape','_rank','_order','_class_type')
     name = 'bool'
     def __init__(self, arg):
         self._shape = arg.shape
         self._rank  = arg.rank
         self._order = arg.order
+        self._class_type = arg.class_type
         super().__init__(arg)
 
     @property
@@ -230,12 +255,13 @@ class NumpyInt(PythonInt):
     arg : TypedAstNode
         The argument passed to the function.
     """
-    __slots__ = ('_shape','_rank','_order')
+    __slots__ = ('_shape','_rank','_order','_class_type')
     name = 'int'
     def __init__(self, arg=None, base=10):
         self._shape = arg.shape
         self._rank  = arg.rank
         self._order = arg.order
+        self._class_type = arg.class_type
         super().__init__(arg)
 
     @property
@@ -278,13 +304,20 @@ class NumpyInt64(NumpyInt):
 
 #==============================================================================
 class NumpyReal(PythonReal):
-    """Represents a call to  numpy.real for code generation.
+    """
+    Represents a call to numpy.real for code generation.
 
+    Represents a call to the NumPy function real.
     > a = 1+2j
     > np.real(a)
     1.0
+
+    Parameters
+    ----------
+    arg : TypedAstNode
+        The argument passed to the function.
     """
-    __slots__ = ('_precision','_rank','_shape','_order')
+    __slots__ = ('_precision','_rank','_shape','_order','_class_type')
     name = 'real'
     def __new__(cls, arg):
         if isinstance(arg.dtype, NativeBool):
@@ -296,8 +329,9 @@ def __init__(self, arg):
         super().__init__(arg)
         self._precision = arg.precision
         self._order = arg.order
-        self._rank  = self.internal_var.rank
+        self._rank  = arg.rank
         self._shape = process_shape(self._rank == 0, self.internal_var.shape)
+        self._class_type = arg.class_type
 
     @property
     def is_elemental(self):
@@ -309,13 +343,20 @@ def is_elemental(self):
 #==============================================================================
 
 class NumpyImag(PythonImag):
-    """Represents a call to  numpy.imag for code generation.
+    """
+    Represents a call to numpy.imag for code generation.
 
+    Represents a call to the NumPy function imag.
     > a = 1+2j
     > np.imag(a)
     2.0
+
+    Parameters
+    ----------
+    arg : TypedAstNode
+        The argument passed to the function.
     """
-    __slots__ = ('_precision','_rank','_shape','_order')
+    __slots__ = ('_precision','_rank','_shape','_order','_class_type')
     name = 'imag'
     def __new__(cls, arg):
 
@@ -331,8 +372,9 @@ def __init__(self, arg):
         super().__init__(arg)
         self._precision = arg.precision
         self._order = arg.order
-        self._rank  = self.internal_var.rank
+        self._rank  = arg.rank
         self._shape = process_shape(self._rank == 0, self.internal_var.shape)
+        self._class_type = arg.class_type
 
     @property
     def is_elemental(self):
@@ -358,7 +400,7 @@ class NumpyComplex(PythonComplex):
     """
     _real_cast = NumpyReal
     _imag_cast = NumpyImag
-    __slots__ = ('_rank','_shape','_order')
+    __slots__ = ('_rank','_shape','_order','_class_type')
     name = 'complex'
     def __init__(self, arg0, arg1 = None):
         if arg1 is not None:
@@ -366,6 +408,7 @@ def __init__(self, arg0, arg1 = None):
         self._shape = arg0.shape
         self._rank  = arg0.rank
         self._order = arg0.order
+        self._class_type = arg0.class_type
         super().__init__(arg0)
 
     @property
@@ -428,7 +471,7 @@ class NumpyResultType(PyccelInternalFunction):
         Any arrays and dtypes passed to the function (currently only accepts one array
         and no dtypes).
     """
-    __slots__ = ('_dtype','_precision')
+    __slots__ = ('_dtype','_precision','_class_type')
     _rank = 0
     _shape = None
     _order = None
@@ -450,6 +493,7 @@ def __init__(self, *arrays_and_dtypes):
             self._dtype     = NativeInteger()
         else:
             self._dtype     = NativeBool()
+        self._class_type = self._dtype
 
         # An integer cannot be described accurately in a float of the same precision due to the mantissa
         # If the user wants to store the maximum value of an integer they will need the next largest precision
@@ -545,10 +589,11 @@ class NumpyNewArray(PyccelInternalFunction):
     init_dtype : PythonType, PyccelFunctionDef, LiteralString, str
         The actual dtype passed to the NumPy function.
     """
-    __slots__ = ('_init_dtype',)
+    __slots__ = ('_init_dtype','_class_type')
 
     def __init__(self, *args, init_dtype = None):
         self._init_dtype = init_dtype
+        self._class_type = NumpyNDArrayType()
 
         super().__init__(*args)
 
@@ -605,11 +650,11 @@ def __init__(self, arg, dtype=None, order='K', ndmin=None):
         if not isinstance(arg, (PythonTuple, PythonList, Variable)):
             raise TypeError('Unknown type of  %s.' % type(arg))
 
-        is_homogeneous_tuple = isinstance(arg, (PythonTuple, PythonList, HomogeneousTupleVariable)) and arg.is_homogeneous
+        is_homogeneous_tuple = isinstance(arg, (PythonTuple, HomogeneousTupleVariable)) and arg.is_homogeneous
         is_array = isinstance(arg, Variable) and arg.is_ndarray
 
         # TODO: treat inhomogenous lists and tuples when they have mixed ordering
-        if not (is_homogeneous_tuple or is_array):
+        if not (is_homogeneous_tuple or is_array or isinstance(arg, PythonList)):
             raise TypeError('we only accept homogeneous arguments')
 
         if not isinstance(order, (LiteralString, str)):
@@ -749,7 +794,7 @@ class NumpySum(PyccelInternalFunction):
     arg : list , tuple , PythonTuple, PythonList, Variable
         The argument passed to the sum function.
     """
-    __slots__ = ('_dtype','_precision')
+    __slots__ = ('_dtype','_precision','_class_type')
     name = 'sum'
     _rank  = 0
     _shape = None
@@ -763,6 +808,7 @@ def __init__(self, arg):
             self._dtype = NativeInteger()
         else:
             self._dtype = arg.dtype
+        self._class_type = self._dtype
         self._precision = max(arg.precision, default_precision[self._dtype])
 
     @property
@@ -781,7 +827,7 @@ class NumpyProduct(PyccelInternalFunction):
     arg : list , tuple , PythonTuple, PythonList, Variable
         The argument passed to the prod function.
     """
-    __slots__ = ('_arg','_dtype','_precision')
+    __slots__ = ('_arg','_dtype','_precision','_class_type')
     name = 'product'
     _rank  = 0
     _shape = None
@@ -796,6 +842,7 @@ def __init__(self, arg):
                     (isinstance(arg.dtype, NativeInteger) and get_final_precision(self._arg) < default_precision[NativeInteger()]))\
                     else self._arg
         self._dtype = self._arg.dtype
+        self._class_type = self._dtype
         self._precision = get_final_precision(self._arg)
 
     @property
@@ -817,7 +864,7 @@ class NumpyMatmul(PyccelInternalFunction):
     b : TypedAstNode
         The second argument of the matrix multiplication.
     """
-    __slots__ = ('_dtype','_precision','_shape','_rank','_order')
+    __slots__ = ('_dtype','_precision','_shape','_rank','_order','_class_type')
     name = 'matmul'
 
     def __init__(self, a ,b):
@@ -856,6 +903,8 @@ def __init__(self, a ,b):
         else:
             self._order = None if self._rank < 2 else 'C'
 
+        self._class_type = NumpyNDArrayType()
+
     @property
     def a(self):
         return self._args[0]
@@ -985,10 +1034,6 @@ def __init__(self, start, stop, num=None, endpoint=True, dtype=None):
 
         super().__init__(init_dtype = init_dtype)
 
-    @property
-    def dtype(self):
-        return self._dtype
-
     @property
     def endpoint(self):
         """Indicates if the stop must be included or not."""
@@ -1014,10 +1059,6 @@ def index(self):
         """Used in the fortran codegen when there is no for loop created."""
         return self._index
 
-    @property
-    def rank(self):
-        return self._rank
-
     @property
     def step(self):
         """Represent size of spacing between generated elements."""
@@ -1058,7 +1099,7 @@ class NumpyWhere(PyccelInternalFunction):
     """
 
     __slots__ = ('_condition', '_value_true', '_value_false', '_dtype',
-                 '_rank', '_shape', '_order', '_precision')
+                 '_rank', '_shape', '_order', '_precision','_class_type')
     _attribute_nodes = ('_condition','_value_true','_value_false')
     name = 'where'
 
@@ -1086,6 +1127,8 @@ def __init__(self, condition, x, y):
         self._shape = process_shape(False, shape)
         self._rank  = len(shape)
         self._order = None if self._rank < 2 else 'C'
+
+        self._class_type = NumpyNDArrayType()
         super().__init__(condition, x, y)
 
     @property
@@ -1122,7 +1165,7 @@ class NumpyRand(PyccelInternalFunction):
     *args : tuple of TypedAstNode
         The arguments passed to the function.
     """
-    __slots__ = ('_shape','_rank','_order')
+    __slots__ = ('_shape','_rank','_order','_class_type')
     name = 'rand'
     _dtype = NativeFloat()
     _precision = default_precision[NativeFloat()]
@@ -1132,15 +1175,27 @@ def __init__(self, *args):
         self._rank  = len(args)
         self._shape = None if self._rank == 0 else args
         self._order = None if self._rank < 2 else 'C'
+        self._class_type = NumpyNDArrayType()
 
 #==============================================================================
 class NumpyRandint(PyccelInternalFunction):
-
     """
-      Represents a call to  numpy.random.random or numpy.random.rand for code generation.
+    Class representing a call to NumPy's randint function.
 
-    """
-    __slots__ = ('_rand','_low','_high','_shape','_rank','_order')
+    Class representing a call to NumPy's randint function.
+
+    Parameters
+    ----------
+    low : TypedAstNode
+        The first argument passed to the function. The smallest possible value for
+        the generated number.
+    high : TypedAstNode, optional
+        The second argument passed to the function. The largest possible value for
+        the generated number.
+    size : TypedAstNode, optional
+        The size of the array that will be generated.
+    """
+    __slots__ = ('_rand','_low','_high','_shape','_rank','_order','_class_type')
     name = 'randint'
     _dtype     = NativeInteger()
     _precision = -1
@@ -1158,6 +1213,7 @@ def __init__(self, low, high = None, size = None):
         self._rank    = 0 if size is None else len(self.shape)
         self._order   = None if self._rank < 2 else 'C'
         self._rand    = NumpyRand() if size is None else NumpyRand(*size)
+        self._class_type = NumpyNDArrayType()
         self._low     = low
         self._high    = high
         super().__init__()
@@ -1506,8 +1562,20 @@ def __new__(cls, a, dtype=None, order='K', subok=True, shape=None):
 
 #==============================================================================
 class NumpyNorm(PyccelInternalFunction):
-    """ Represents call to numpy.norm"""
-    __slots__ = ('_shape','_rank','_order','_arg','_precision')
+    """
+    Represents call to `numpy.norm`.
+
+    Represents a call to the NumPy function norm.
+
+    Parameters
+    ----------
+    arg : TypedAstNode
+        The first argument passed to the function.
+    axis : TypedAstNode, optional
+        The second argument passed to the function, indicating the axis along
+        which the norm should be calculated.
+    """
+    __slots__ = ('_shape','_rank','_order','_arg','_precision','_class_type')
     name = 'norm'
     _dtype = NativeFloat()
 
@@ -1527,6 +1595,7 @@ def __init__(self, arg, axis=None):
             self._shape = None
             self._order = None
             self._rank  = 0
+        self._class_type = NumpyNDArrayType()
 
     @property
     def arg(self):
@@ -1553,8 +1622,18 @@ def axis(self):
 # https://numpy.org/doc/stable/reference/ufuncs.html#available-ufuncs
 #==============================================================================
 class NumpyUfuncBase(PyccelInternalFunction):
-    """Base class for Numpy's universal functions."""
-    __slots__ = ('_dtype','_precision','_shape','_rank','_order')
+    """
+    Base class for Numpy's universal functions.
+
+    The class from which NumPy's universal functions inherit. All classes which
+    inherit from this class operate on their arguments elementally.
+
+    Parameters
+    ----------
+    *args : tuple of TypedAstNode
+        The arguments passed to the function.
+    """
+    __slots__ = ('_dtype','_precision','_shape','_rank','_order','_class_type')
     @property
     def is_elemental(self):
         return True
@@ -1582,6 +1661,7 @@ def __init__(self, x):
         self._set_dtype_precision(x)
         self._set_shape_rank(x)
         self._set_order(x)
+        self._class_type = NumpyNDArrayType()
         super().__init__(x)
 
     def _set_shape_rank(self, x):
@@ -1631,6 +1711,7 @@ def __init__(self, x1, x2):
         self._set_dtype_precision(x1, x2)
         self._set_shape_rank(x1, x2)
         self._set_order(x1, x2)
+        self._class_type = NumpyNDArrayType()
 
     def _set_shape_rank(self, x1, x2):
         self._shape = broadcast(x1.shape, x2.shape)
@@ -1939,13 +2020,20 @@ def is_elemental(self):
         return False
 
 class NumpyConjugate(PythonConjugate):
-    """Represents a call to  numpy.conj for code generation.
+    """
+    Represents a call to  numpy.conj for code generation.
 
+    Represents a call to the NumPy function conj or conjugate.
     > a = 1+2j
     > np.conj(a)
     1-2j
+
+    Parameters
+    ----------
+    arg : TypedAstNode
+        The argument passed to the function.
     """
-    __slots__ = ('_precision','_rank','_shape','_order')
+    __slots__ = ('_precision','_rank','_shape','_order','_class_type')
     name = 'conj'
 
     def __init__(self, arg):
@@ -1954,6 +2042,7 @@ def __init__(self, arg):
         self._order = arg.order
         self._rank  = self.internal_var.rank
         self._shape = process_shape(self._rank == 0, self.internal_var.shape)
+        self._class_type = NumpyNDArrayType()
 
     @property
     def is_elemental(self):
@@ -2053,20 +2142,23 @@ def __iter__(self):
 
 class NumpyCountNonZero(PyccelInternalFunction):
     """
-    Class representing a call to the numpy size function which
-    returns the shape of an object in a given dimension
+    Class representing a call to the NumPy function `count_nonzero`.
+
+    Class representing a call to the NumPy function `count_nonzero` which
+    counts the number of non-zero elements in an array.
 
     Parameters
-    ==========
-    arg   : TypedAstNode
-            An array for which the non-zero elements should be counted
-    axis  : int
-            The dimension along which the non-zero elements are counted
-    keep_dims : NativeBool
-            Indicates if output arrays should have the same number of dimensions
-            as arg
-    """
-    __slots__ = ('_precision', '_rank', '_shape', '_order', '_arr', '_axis', '_keep_dims')
+    ----------
+    a : TypedAstNode
+        An array for which the non-zero elements should be counted.
+    axis : int, optional
+        The dimension along which the non-zero elements are counted.
+    keepdims : NativeBool
+        Indicates if output arrays should have the same number of dimensions
+        as arg.
+    """
+    __slots__ = ('_precision', '_rank', '_shape', '_order', '_class_type', '_arr',
+                '_axis', '_keep_dims')
     _attribute_nodes = ('_arr','_axis')
     name   = 'count_nonzero'
     _dtype = NativeInteger()
@@ -2097,6 +2189,8 @@ def __init__(self, a, axis = None, *, keepdims = LiteralFalse()):
                 self._shape = None
                 self._order = None
 
+        self._class_type = NumpyNDArrayType()
+
         self._arr = a
         self._axis = axis
         self._keep_dims = keepdims
diff --git a/pyccel/ast/operators.py b/pyccel/ast/operators.py
index 1ef01bf4ba..463d50564f 100644
--- a/pyccel/ast/operators.py
+++ b/pyccel/ast/operators.py
@@ -175,7 +175,7 @@ def _set_dtype(self):
         uses the static method `_calculate_dtype` to set these values. If the
         values are class parameters in a sub-class, this method must be over-ridden.
         """
-        self._dtype, self._precision = self._calculate_dtype(*self._args)  # pylint: disable=no-member
+        self._dtype, self._precision, self._class_type = self._calculate_dtype(*self._args)  # pylint: disable=no-member
 
     def _set_shape_rank(self):
         """
@@ -290,20 +290,38 @@ class PyccelUnaryOperator(PyccelOperator):
     arg : TypedAstNode
         The argument passed to the operator.
     """
-    __slots__ = ('_dtype', '_precision','_shape','_rank','_order')
+    __slots__ = ('_dtype', '_precision','_shape','_rank','_order','_class_type')
 
     def __init__(self, arg):
         super().__init__(arg)
 
     @staticmethod
-    def _calculate_dtype(*args):
-        """ Sets the dtype and precision
-        They are chosen to match the argument
+    def _calculate_dtype(arg):
         """
-        a = args[0]
-        dtype = a.dtype
-        precision = a.precision
-        return dtype, precision
+        Calculate the dtype, precision and class type of the result.
+
+        Calculate the dtype, precision and class type of the result.
+        These are equivalent to the dtype, precision and class type
+        of the only argument.
+
+        Parameters
+        ----------
+        arg : TypedAstNode
+            The argument passed to the operator.
+
+        Returns
+        -------
+        dtype : DataType
+            The underlying datatype of the object.
+        precision : int
+            The precision of the datatype of the object.
+        class_type : DataType
+            The Python type of the object.
+        """
+        dtype = arg.dtype
+        precision = arg.precision
+        class_type = arg.class_type
+        return dtype, precision, class_type
 
     @staticmethod
     def _calculate_shape_rank(*args):
@@ -385,14 +403,31 @@ class PyccelNot(PyccelUnaryOperator):
     _precedence = 6
 
     @staticmethod
-    def _calculate_dtype(*args):
-        """ Sets the dtype and precision
-        They are chosen to match the argument unless the class has
-        a _dtype or _precision member
+    def _calculate_dtype(arg):
+        """
+        Calculate the dtype, precision and class type of the result.
+
+        Calculate the dtype, precision and class type of the result.
+        These are the dtype, precision and class type which represent a
+        boolean.
+
+        Parameters
+        ----------
+        arg : TypedAstNode
+            The argument passed to the operator.
+
+        Returns
+        -------
+        dtype : DataType
+            The underlying datatype of the object.
+        precision : int
+            The precision of the datatype of the object.
+        class_type : DataType
+            The Python type of the object.
         """
         dtype = NativeBool()
         precision = -1
-        return dtype, precision
+        return dtype, precision, dtype
 
     @staticmethod
     def _calculate_shape_rank(*args):
@@ -445,13 +480,13 @@ class PyccelBinaryOperator(PyccelOperator):
     arg2 : TypedAstNode
         The second argument passed to the operator.
     """
-    __slots__ = ('_dtype','_precision','_shape','_rank','_order')
+    __slots__ = ('_dtype','_precision','_shape','_rank','_order','_class_type')
 
     def __init__(self, arg1, arg2):
         super().__init__(arg1, arg2)
 
     @classmethod
-    def _calculate_dtype(cls, *args):
+    def _calculate_dtype(cls, arg1, arg2):
         """
         Sets the dtype and precision.
 
@@ -464,25 +499,38 @@ def _calculate_dtype(cls, *args):
 
         Parameters
         ----------
-        *args : tuple of TypedAstNodes
-            The arguments of the operators.
-        """
-        integers  = [a for a in args if a.dtype in (NativeInteger(),NativeBool())]
-        floats    = [a for a in args if a.dtype is NativeFloat()]
-        complexes = [a for a in args if a.dtype is NativeComplex()]
-        strs      = [a for a in args if a.dtype is NativeString()]
+        arg1 : TypedAstNode
+            The first argument passed to the operator.
+        arg2 : TypedAstNode
+            The second argument passed to the operator.
 
-        if strs:
-            assert len(integers + floats + complexes) == 0
-            return cls._handle_str_type(strs)
-        elif complexes:
-            return cls._handle_complex_type(args)
-        elif floats:
-            return cls._handle_float_type(args)
-        elif integers:
-            return cls._handle_integer_type(args)
+        Returns
+        -------
+        dtype : DataType
+            The underlying datatype of the object.
+        precision : int
+            The precision of the datatype of the object.
+        class_type : DataType
+            The Python type of the object.
+        """
+        try:
+            dtype = arg1.dtype + arg2.dtype
+            class_type = arg1.class_type + arg2.class_type
+        except NotImplementedError:
+            raise TypeError(f'Cannot determine the type of ({arg1}, {arg2})') #pylint: disable=raise-missing-from
+
+        if dtype is NativeString():
+            return *cls._handle_str_type((arg1, arg2)), class_type
+        elif dtype is NativeComplex():
+            return *cls._handle_complex_type((arg1, arg2)), class_type
+        elif dtype is NativeFloat():
+            return *cls._handle_float_type((arg1, arg2)), class_type
+        elif dtype in (NativeInteger(), NativeBool()):
+            if class_type is NativeBool():
+                class_type = NativeInteger()
+            return *cls._handle_integer_type((arg1, arg2)), class_type
         else:
-            raise TypeError(f'cannot determine the type of {args}')
+            raise TypeError(f'Cannot determine the type of ({arg1}, {arg2})')
 
     @staticmethod
     def _handle_str_type(strs):
@@ -634,7 +682,7 @@ def __new__(cls, arg1, arg2, simplify = False):
         if simplify:
             if isinstance(arg2, PyccelUnarySub):
                 return PyccelMinus(arg1, arg2.args[0], simplify = True)
-            dtype, precision = cls._calculate_dtype(arg1, arg2)
+            dtype, precision, _ = cls._calculate_dtype(arg1, arg2)
             if isinstance(arg1, Literal) and isinstance(arg2, Literal):
                 return convert_to_literal(arg1.python_value + arg2.python_value,
                                           dtype, precision)
@@ -700,14 +748,14 @@ def __new__(cls, arg1, arg2, simplify = False):
             if (arg2 == 1):
                 return arg1
             if (arg1 == 0 or arg2 == 0):
-                dtype, precision = cls._calculate_dtype(arg1, arg2)
+                dtype, precision, _ = cls._calculate_dtype(arg1, arg2)
                 return convert_to_literal(0, dtype, precision)
             if (isinstance(arg1, PyccelUnarySub) and arg1.args[0] == 1):
                 return PyccelUnarySub(arg2)
             if (isinstance(arg2, PyccelUnarySub) and arg2.args[0] == 1):
                 return PyccelUnarySub(arg1)
             if isinstance(arg1, Literal) and isinstance(arg2, Literal):
-                dtype, precision = cls._calculate_dtype(arg1, arg2)
+                dtype, precision, _ = cls._calculate_dtype(arg1, arg2)
                 return convert_to_literal(arg1.python_value * arg2.python_value,
                                           dtype, precision)
         return super().__new__(cls)
@@ -748,7 +796,7 @@ def __new__(cls, arg1, arg2, simplify = False):
             if isinstance(arg2, PyccelUnarySub):
                 return PyccelAdd(arg1, arg2.args[0], simplify = True)
             elif isinstance(arg1, Literal) and isinstance(arg2, Literal):
-                dtype, precision = cls._calculate_dtype(arg1, arg2)
+                dtype, precision, _ = cls._calculate_dtype(arg1, arg2)
                 return convert_to_literal(arg1.python_value - arg2.python_value,
                                           dtype, precision)
         if isinstance(arg1, LiteralFloat) and \
@@ -881,9 +929,30 @@ class PyccelComparisonOperator(PyccelBinaryOperator):
     _precedence = 7
     @staticmethod
     def _calculate_dtype(*args):
+        """
+        Calculate the dtype, precision and class type of the result.
+
+        Calculate the dtype, precision and class type of the result.
+        These are the dtype, precision and class type which represent a
+        boolean.
+
+        Parameters
+        ----------
+        *args : TypedAstNode
+            The arguments passed to the operator.
+
+        Returns
+        -------
+        dtype : DataType
+            The underlying datatype of the object.
+        precision : int
+            The precision of the datatype of the object.
+        class_type : DataType
+            The Python type of the object.
+        """
         dtype = NativeBool()
         precision = -1
-        return dtype, precision
+        return dtype, precision, dtype
 
 #==============================================================================
 
@@ -1057,11 +1126,12 @@ class PyccelBooleanOperator(PyccelOperator):
     *args : tuple of TypedAstNode
         The arguments passed to the operator.
     """
-    dtype = NativeBool()
-    precision = -1
-    rank = 0
-    shape = None
-    order = None
+    _dtype = NativeBool()
+    _precision = -1
+    _rank = 0
+    _shape = None
+    _order = None
+    _class_type = NativeBool()
 
     __slots__ = ()
 
@@ -1253,7 +1323,7 @@ class IfTernaryOperator(PyccelOperator):
     >>> IfTernaryOperator(PyccelGt(n > 1),  5,  2)
     IfTernaryOperator(PyccelGt(n > 1),  5,  2)
     """
-    __slots__ = ('_dtype','_precision','_shape','_rank','_order')
+    __slots__ = ('_dtype','_precision','_shape','_rank','_order','_class_type')
     _precedence = 3
 
     def __init__(self, cond, value_true, value_false):
@@ -1276,7 +1346,31 @@ def __init__(self, cond, value_true, value_false):
     @staticmethod
     def _calculate_dtype(cond, value_true, value_false):
         """
-        Sets the dtype and precision for IfTernaryOperator
+        Calculate the dtype, precision and class type of the result.
+
+        Calculate the dtype, precision and class type of the result. The dtype,
+        precision and class type are calculated from the types of the values if
+        true or false.
+
+        Parameters
+        ----------
+        cond : TypedAstNode
+            The first argument passed to the operator representing the condition.
+        value_true : TypedAstNode
+            The second argument passed to the operator representing the result if the
+            condition is true.
+        value_false : TypedAstNode
+            The third argument passed to the operator representing the result if the
+            condition is false.
+
+        Returns
+        -------
+        dtype : DataType
+            The underlying datatype of the object.
+        precision : int
+            The precision of the datatype of the object.
+        class_type : DataType
+            The Python type of the object.
         """
         if value_true.dtype in NativeNumeric and value_false.dtype in NativeNumeric:
             dtype = max([value_true.dtype, value_false.dtype], key = NativeNumeric.index)
@@ -1284,7 +1378,9 @@ def _calculate_dtype(cond, value_true, value_false):
             dtype = value_true.dtype
 
         precision = max_precision([value_true, value_false])
-        return dtype, precision
+
+        class_type = value_true.class_type + value_false.class_type
+        return dtype, precision, class_type
 
     @staticmethod
     def _calculate_shape_rank(cond, value_true, value_false):
diff --git a/pyccel/ast/type_annotations.py b/pyccel/ast/type_annotations.py
index 92ff3d0a66..3718e98b7e 100644
--- a/pyccel/ast/type_annotations.py
+++ b/pyccel/ast/type_annotations.py
@@ -36,8 +36,10 @@ class VariableTypeAnnotation(PyccelAstNode):
     datatype : DataType
         The requested internal data type.
 
-    cls_base : ClassDef
-        The description of the class describing the variable.
+    class_type : DataType
+        The Python type of the variable. In the case of scalars this is equivalent to
+        the datatype. For objects in (homogeneous) containers (e.g. list/ndarray/tuple),
+        this is the type of the container.
 
     precision : int
         The precision of the internal datatype.
@@ -51,13 +53,13 @@ class VariableTypeAnnotation(PyccelAstNode):
     is_const : bool, default=False
         True if the variable cannot be modified, false otherwise.
     """
-    __slots__ = ('_datatype', '_cls_base', '_precision', '_rank',
+    __slots__ = ('_datatype', '_class_type', '_precision', '_rank',
                  '_order', '_is_const')
     _attribute_nodes = ()
-    def __init__(self, datatype : 'DataType', cls_base : 'ClassDef', precision : int = -1,
+    def __init__(self, datatype : 'DataType', class_type : 'DataType', precision : int = -1,
             rank : int = 0, order : str = None, is_const : bool = False):
         self._datatype = datatype
-        self._cls_base = cls_base
+        self._class_type = class_type
         self._precision = precision
         self._rank = rank
         self._order = order
@@ -76,14 +78,15 @@ def datatype(self):
         return self._datatype
 
     @property
-    def cls_base(self):
+    def class_type(self):
         """
-        Get the class description of the object.
+        Get the Python type of the object.
 
-        Get the class def object which describes how the user can interact with the
-        future variable.
+        The Python type of the object. In the case of scalars this is equivalent to
+        the datatype. For objects in (homogeneous) containers (e.g. list/ndarray/tuple),
+        this is the type of the container.
         """
-        return self._cls_base
+        return self._class_type
 
     @property
     def precision(self):
@@ -126,13 +129,13 @@ def is_const(self):
         return self._is_const
 
     def __hash__(self):
-        return hash((self.datatype, self.cls_base, self.precision, self.rank, self.order))
+        return hash((self.datatype, self.class_type, self.precision, self.rank, self.order))
 
     def __eq__(self, other):
         # Needed for set
         if isinstance(other, VariableTypeAnnotation):
             return self.datatype == other.datatype and \
-                   self.cls_base == other.cls_base and \
+                   self.class_type == other.class_type and \
                    self.precision == other.precision and \
                    self.rank == other.rank and \
                    self.order == other.order
diff --git a/pyccel/ast/variable.py b/pyccel/ast/variable.py
index 8d3e00c84a..a60cf6285e 100644
--- a/pyccel/ast/variable.py
+++ b/pyccel/ast/variable.py
@@ -15,7 +15,7 @@
 from .basic     import PyccelAstNode, TypedAstNode
 from .datatypes import (datatype, DataType,
                         NativeInteger, NativeBool, NativeFloat,
-                        NativeComplex)
+                        NativeComplex, NativeHomogeneousTuple, NativeInhomogeneousTuple)
 from .internals import PyccelArrayShapeElement, Slice, get_final_precision, PyccelSymbol
 from .literals  import LiteralInteger, Nil
 from .operators import (PyccelMinus, PyccelDiv, PyccelMul,
@@ -53,6 +53,11 @@ class Variable(TypedAstNode):
         The name of the variable represented. This can be either a string
         or a dotted name, when using a Class attribute.
 
+    class_type : DataType
+        The Python type of the variable. In the case of scalars this is equivalent to
+        the datatype. For objects in (homogeneous) containers (e.g. list/ndarray/tuple),
+        this is the type of the container.
+
     rank : int, default: 0
         The number of dimensions for an array.
 
@@ -112,7 +117,7 @@ class Variable(TypedAstNode):
     __slots__ = ('_name', '_alloc_shape', '_memory_handling', '_is_const',
             '_is_target', '_is_optional', '_allows_negative_indexes',
             '_cls_base', '_is_argument', '_is_temp','_dtype','_precision',
-            '_rank','_shape','_order','_is_private')
+            '_rank','_shape','_order','_is_private','_class_type')
     _attribute_nodes = ()
 
     def __init__(
@@ -120,6 +125,7 @@ def __init__(
         dtype,
         name,
         *,
+        class_type = None,
         rank=0,
         memory_handling='stack',
         is_const=False,
@@ -211,11 +217,20 @@ def __init__(
         if not isinstance(precision,int) and precision is not None:
             raise TypeError('precision must be an integer or None.')
 
+        if rank > 0 and class_type is None:
+            raise TypeError("Multi-dimensional object requires a container type")
+        elif class_type is None:
+            class_type = dtype
+
+        if class_type is not dtype and cls_base is None:
+            raise TypeError(f"Missing class definition for type {dtype}")
+
         self._alloc_shape = shape
         self._dtype = dtype
         self._rank  = rank
         self._shape = self.process_shape(shape)
         self._precision = precision
+        self._class_type = class_type
         if self._rank < 2:
             self._order = None
 
@@ -674,13 +689,19 @@ def is_ndarray(self):
         return False
 
 class HomogeneousTupleVariable(TupleVariable):
+    """
+    Represents a homogeneous tuple variable in the code.
 
-    """Represents a tuple variable in the code.
+    Represents a homogeneous tuple variable in the code.
 
     Parameters
     ----------
-    arg_vars: Iterable
-        Multiple variables contained within the tuple
+    dtype : DataType
+        The data type of the elements of the tuple.
+    *args : tuple
+        See Variable.
+    **kwargs : dict
+        See Variable.
 
     Examples
     --------
@@ -694,6 +715,9 @@ class HomogeneousTupleVariable(TupleVariable):
     __slots__ = ()
     is_homogeneous = True
 
+    def __init__(self, dtype, *args, **kwargs):
+        super().__init__(dtype, *args, **kwargs)
+
     def shape_can_change(self, i):
         """
         Indicates if the shape can change in the i-th dimension
@@ -708,13 +732,25 @@ def __iter__(self):
         return (self[i] for i in range(self.shape[0]))
 
 class InhomogeneousTupleVariable(TupleVariable):
+    """
+    Represents an inhomogeneous tuple variable in the code.
 
-    """Represents a tuple variable in the code.
+    Represents an inhomogeneous tuple variable in the code.
 
     Parameters
     ----------
-    arg_vars: Iterable
-        Multiple variables contained within the tuple
+    arg_vars : tuple of Variable
+        The variables contained within the tuple.
+    name : str
+        The name of the variable.
+    *args : tuple
+        See Variable.
+    class_type : DataType
+        The Python type of the variable. In the case of scalars this is equivalent to
+        the datatype. For objects in (homogeneous) containers (e.g. list/ndarray/tuple),
+        this is the type of the container.
+    **kwargs : dict
+        See Variable.
 
     Examples
     --------
@@ -729,9 +765,9 @@ class InhomogeneousTupleVariable(TupleVariable):
     _attribute_nodes = ('_vars',)
     is_homogeneous = False
 
-    def __init__(self, arg_vars, dtype, name, *args, **kwargs):
+    def __init__(self, arg_vars, name, *args, class_type, **kwargs):
         self._vars = tuple(arg_vars)
-        super().__init__(dtype, name, *args, **kwargs)
+        super().__init__(class_type, name, *args, **kwargs, class_type = class_type)
 
     def get_vars(self):
         """ Get the variables saved internally in the tuple
@@ -879,7 +915,7 @@ class IndexedElement(TypedAstNode):
     >>> IndexedElement(A, i, j) == A[i, j]
     True
     """
-    __slots__ = ('_label', '_indices','_dtype','_precision','_shape','_rank','_order')
+    __slots__ = ('_label', '_indices','_dtype','_precision','_shape','_rank','_order','_class_type')
     _attribute_nodes = ('_label', '_indices')
 
     def __init__(self, base, *indices):
@@ -938,6 +974,11 @@ def __init__(self, base, *indices):
 
         self._order = None if self.rank < 2 else base.order
 
+        if self.rank == 0:
+            self._class_type = self.dtype
+        else:
+            self._class_type = base.class_type
+
     @property
     def base(self):
         """ The object which is indexed
diff --git a/pyccel/codegen/printing/fcode.py b/pyccel/codegen/printing/fcode.py
index 85fe6f2fe5..31dd47a6d2 100644
--- a/pyccel/codegen/printing/fcode.py
+++ b/pyccel/codegen/printing/fcode.py
@@ -37,7 +37,7 @@
 from pyccel.ast.datatypes import NativeInteger, NativeBool, NativeFloat, NativeComplex
 from pyccel.ast.datatypes import iso_c_binding
 from pyccel.ast.datatypes import iso_c_binding_shortcut_mapping
-from pyccel.ast.datatypes import NativeRange, NativeNumeric
+from pyccel.ast.datatypes import NativeNumeric
 from pyccel.ast.datatypes import CustomDataType
 
 from pyccel.ast.internals import Slice, PrecomputedCode, PyccelArrayShapeElement
@@ -667,7 +667,7 @@ def _print_PythonPrint(self, expr):
             elif isinstance(f, PythonType):
                 args_format.append('A')
                 args.append(self._print(f.print_string))
-            elif isinstance(f.rank, int) and f.rank > 0:
+            elif f.rank > 0 and not isinstance(f, FunctionCall):
                 if args_format:
                     code += self._formatted_args_to_print(args_format, args, sep, separator, expr)
                     args_format = []
@@ -1360,9 +1360,6 @@ def _print_Declare(self, expr):
         if isinstance(expr.dtype, NativeSymbol):
             return ''
 
-        if isinstance(expr.dtype, NativeRange):
-            return ''
-
         # meta-variables
         if (isinstance(expr.variable, Variable) and
             expr.variable.name.startswith('__')):
diff --git a/pyccel/parser/semantic.py b/pyccel/parser/semantic.py
index a6f6fa3beb..9a28ddd8f9 100644
--- a/pyccel/parser/semantic.py
+++ b/pyccel/parser/semantic.py
@@ -24,7 +24,7 @@
 from pyccel.ast.basic import PyccelAstNode, TypedAstNode, ScopedAstNode
 
 from pyccel.ast.builtins import PythonPrint
-from pyccel.ast.builtins import PythonInt, PythonBool, PythonFloat, PythonComplex
+from pyccel.ast.builtins import PythonComplex
 from pyccel.ast.builtins import python_builtin_datatype, PythonImag, PythonReal
 from pyccel.ast.builtins import PythonList, PythonConjugate
 from pyccel.ast.builtins import (PythonRange, PythonZip, PythonEnumerate,
@@ -61,12 +61,12 @@
 
 from pyccel.ast.class_defs import NumpyArrayClass, TupleClass, get_cls_base
 
-from pyccel.ast.datatypes import NativeRange, str_dtype
+from pyccel.ast.datatypes import str_dtype
 from pyccel.ast.datatypes import NativeSymbol, DataTypeFactory, CustomDataType
 from pyccel.ast.datatypes import default_precision, dtype_and_precision_registry
-from pyccel.ast.datatypes import (NativeInteger, NativeBool,
-                                  NativeFloat, NativeString,
-                                  NativeGeneric, NativeComplex,
+from pyccel.ast.datatypes import (NativeInteger, NativeBool, NativeHomogeneousList,
+                                  NativeFloat, NativeString, NativeInhomogeneousTuple,
+                                  NativeGeneric, NativeComplex, NativeTuple,
                                   NativeVoid)
 
 from pyccel.ast.functionalexpr import FunctionalSum, FunctionalMax, FunctionalMin, GeneratorComprehension, FunctionalFor
@@ -85,14 +85,10 @@
 from pyccel.ast.mathext  import math_constants, MathSqrt, MathAtan2, MathSin, MathCos
 
 from pyccel.ast.numpyext import NumpyMatmul
-from pyccel.ast.numpyext import NumpyBool
 from pyccel.ast.numpyext import NumpyWhere, NumpyArray
-from pyccel.ast.numpyext import NumpyInt, NumpyInt8, NumpyInt16, NumpyInt32, NumpyInt64
-from pyccel.ast.numpyext import NumpyFloat, NumpyFloat32, NumpyFloat64
-from pyccel.ast.numpyext import NumpyComplex, NumpyComplex64, NumpyComplex128
 from pyccel.ast.numpyext import NumpyTranspose, NumpyConjugate
 from pyccel.ast.numpyext import NumpyNewArray, NumpyNonZero, NumpyResultType
-from pyccel.ast.numpyext import DtypePrecisionToCastFunction
+from pyccel.ast.numpyext import DtypePrecisionToCastFunction, NumpyNDArrayType
 
 from pyccel.ast.omp import (OMP_For_Loop, OMP_Simd_Construct, OMP_Distribute_Construct,
                             OMP_TaskLoop_Construct, OMP_Sections_Construct, Omp_End_Clause,
@@ -594,51 +590,28 @@ def _infer_type(self, expr):
         dict
             Dictionary containing all the type information which was inferred.
         """
-        # TODO - add settings to Errors
-        #      - line and column
-        #      - blocking errors
-
-        errors = Errors()
-
-        d_var = {}
-        # TODO improve => put settings as attribut of Parser
-
-        if expr in (PythonInt, PythonFloat, PythonComplex, PythonBool, NumpyBool, NumpyInt, NumpyInt8, NumpyInt16,
-                      NumpyInt32, NumpyInt64, NumpyComplex, NumpyComplex64,
-                      NumpyComplex128, NumpyFloat, NumpyFloat64, NumpyFloat32):
-
-            d_var['datatype'   ] = '*'
-            d_var['rank'       ] = 0
-            d_var['precision'  ] = 0
-            return d_var
-
-        elif isinstance(expr, Variable):
-            d_var['datatype'       ] = expr.dtype
+        d_var = {
+                'datatype' : expr.dtype,
+                'precision': expr.precision,
+                'shape'    : expr.shape,
+                'rank'     : expr.rank,
+                'order'    : expr.order,
+                'class_type' : expr.class_type,
+            }
+
+        if isinstance(expr, Variable):
             d_var['memory_handling'] = expr.memory_handling
-            d_var['shape'          ] = expr.shape
-            d_var['rank'           ] = expr.rank
+            d_var['class_type'     ] = expr.class_type
             d_var['cls_base'       ] = expr.cls_base or self.scope.find(expr.dtype.name, 'classes')
             d_var['is_target'      ] = expr.is_target
-            d_var['order'          ] = expr.order
-            d_var['precision'      ] = expr.precision
             return d_var
 
         elif isinstance(expr, PythonTuple):
-            d_var['datatype'       ] = expr.dtype
-            d_var['precision'      ] = expr.precision
-            d_var['memory_handling'] = 'heap'
-            d_var['shape'          ] = expr.shape
-            d_var['rank'           ] = expr.rank
-            d_var['order'          ] = expr.order
             d_var['cls_base'       ] = TupleClass
+            d_var['memory_handling'] = 'heap'
             return d_var
 
         elif isinstance(expr, Concatenate):
-            d_var['datatype'      ] = expr.dtype
-            d_var['precision'     ] = expr.precision
-            d_var['shape'         ] = expr.shape
-            d_var['rank'          ] = expr.rank
-            d_var['order'         ] = expr.order
             d_var['cls_base'      ] = TupleClass
             if any(getattr(a, 'on_heap', False) for a in expr.args):
                 d_var['memory_handling'] = 'heap'
@@ -648,13 +621,6 @@ def _infer_type(self, expr):
 
         elif isinstance(expr, Duplicate):
             d = self._infer_type(expr.val)
-
-            # TODO must check that it is consistent with pyccel's rules
-            # TODO improve
-            d_var['datatype'      ] = d['datatype']
-            d_var['rank'          ] = expr.rank
-            d_var['shape'         ] = expr.shape
-            d_var['order'         ] = expr.order
             d_var['cls_base'      ] = TupleClass
             if d.get('on_stack', False) and isinstance(expr.length, LiteralInteger):
                 d_var['memory_handling'] = 'stack'
@@ -663,54 +629,23 @@ def _infer_type(self, expr):
             return d_var
 
         elif isinstance(expr, NumpyNewArray):
-            d_var['datatype'   ] = expr.dtype
-            d_var['memory_handling'] = 'heap' if expr.rank > 0 else 'stack'
-            d_var['shape'      ] = expr.shape
-            d_var['rank'       ] = expr.rank
-            d_var['order'      ] = expr.order
-            d_var['precision'  ] = expr.precision
             d_var['cls_base'   ] = NumpyArrayClass
+            d_var['memory_handling'] = 'heap' if expr.rank > 0 else 'stack'
             return d_var
 
         elif isinstance(expr, NumpyTranspose):
 
             var = expr.internal_var
 
-            d_var['memory_handling'] = 'alias' if isinstance(var, Variable) else 'heap'
-            d_var['datatype'      ] = var.dtype
-            d_var['shape'         ] = tuple(reversed(var.shape))
-            d_var['rank'          ] = var.rank
             d_var['cls_base'      ] = var.cls_base
             d_var['is_target'     ] = var.is_target
-            d_var['order'         ] = 'C' if var.order=='F' else 'F'
-            d_var['precision'     ] = var.precision
+            d_var['memory_handling'] = 'alias' if isinstance(var, Variable) else 'heap'
             return d_var
 
         elif isinstance(expr, TypedAstNode):
 
-            d_var['datatype'   ] = expr.dtype
             d_var['memory_handling'] = 'heap' if expr.rank > 0 else 'stack'
-            d_var['shape'      ] = expr.shape
-            d_var['rank'       ] = expr.rank
-            d_var['order'      ] = expr.order
-            d_var['precision'  ] = expr.precision
-            d_var['cls_base'   ] = get_cls_base(expr.dtype, expr.precision, expr.rank)
-            return d_var
-
-        elif isinstance(expr, PythonRange):
-
-            d_var['datatype'   ] = NativeRange()
-            d_var['memory_handling'] = 'stack' # because rank is 0 and no shape defined
-            d_var['shape'      ] = None
-            d_var['rank'       ] = 0
-            d_var['cls_base'   ] = expr  # TODO: shall we keep it?
-            return d_var
-
-        elif isinstance(expr, Lambda):
-
-            d_var['datatype'   ] = NativeSymbol()
-            d_var['memory_handling'] = 'stack' # because rank is 0 and no shape defined
-            d_var['rank'       ] = 0
+            d_var['cls_base'   ] = get_cls_base(expr.dtype, expr.precision, expr.class_type)
             return d_var
 
         else:
@@ -887,12 +822,12 @@ def _create_Duplicate(self, val, length):
         """
         # Arguments have been visited in PyccelMul
 
-        if not isinstance(val, (TupleVariable, PythonTuple)):
+        if not isinstance(val, (TupleVariable, PythonTuple, PythonList)):
             errors.report("Unexpected Duplicate", symbol=Duplicate(val, length),
                 bounding_box=(self.current_ast_node.lineno, self.current_ast_node.col_offset),
                 severity='fatal')
 
-        if val.is_homogeneous:
+        if getattr(val, 'is_homogeneous', True):
             return Duplicate(val, length)
         else:
             if isinstance(length, LiteralInteger):
@@ -1173,7 +1108,7 @@ def _create_variable(self, name, dtype, rhs, d_lhs, arr_in_multirets=False):
             if is_homogeneous and not (d_lhs['memory_handling'] == 'alias' and isinstance(rhs, PythonTuple)):
                 lhs = HomogeneousTupleVariable(dtype, name, **d_lhs, is_temp=is_temp)
             else:
-                lhs = InhomogeneousTupleVariable(elem_vars, dtype, name, **d_lhs, is_temp=is_temp)
+                lhs = InhomogeneousTupleVariable(elem_vars, name, **d_lhs, is_temp=is_temp)
 
         else:
             new_type = HomogeneousTupleVariable \
@@ -1785,8 +1720,10 @@ def _PyccelAstNode_to_TypeAnnotation(self, expr, input_rank = 0, input_order = N
                 order = None
         if rank > 1 and order is None:
             order = 'C'
-        cls_base = get_cls_base(dtype, prec, rank)
-        return VariableTypeAnnotation(dtype, cls_base, prec, rank, order)
+
+        class_type = dtype if rank == 0 else expr.static_class_type()
+
+        return VariableTypeAnnotation(dtype, class_type, prec, rank, order)
 
     #====================================================
     #                 _visit functions
@@ -1970,13 +1907,15 @@ def _visit_Module(self, expr):
                         for v in headers:
                             types = [self._visit(d).type_list[0] for d in v.dtypes]
                             args = [Variable(t.datatype, PyccelSymbol(f'anon_{i}'), precision = t.precision,
-                                shape = None, rank = t.rank, order = t.order, cls_base = t.cls_base,
+                                shape = None, rank = t.rank, order = t.order, class_type = t.class_type,
                                 is_const = t.is_const, is_optional = False,
+                                cls_base = t.datatype if t.rank == 0 else NumpyNDArrayType(),
                                 memory_handling = 'heap' if t.rank > 0 else 'stack') for i,t in enumerate(types)]
 
                             types = [self._visit(d).type_list[0] for d in v.results]
                             results = [Variable(t.datatype, PyccelSymbol(f'result_{i}'), precision = t.precision,
-                                shape = None, rank = t.rank, order = t.order, cls_base = t.cls_base,
+                                shape = None, rank = t.rank, order = t.order, class_type = t.class_type,
+                                cls_base = t.datatype if t.rank == 0 else NumpyNDArrayType(),
                                 is_const = t.is_const, is_optional = False,
                                 memory_handling = 'heap' if t.rank > 0 else 'stack') for i,t in enumerate(types)]
 
@@ -2036,9 +1975,9 @@ def _visit_PythonTuple(self, expr):
 
     def _visit_PythonList(self, expr):
         ls = [self._visit(i) for i in expr]
-        expr = PythonList(*ls)
-
-        if not expr.is_homogeneous:
+        try:
+            expr = PythonList(*ls)
+        except TypeError:
             errors.report(PYCCEL_RESTRICTION_INHOMOG_LIST, symbol=expr,
                 severity='fatal')
         return expr
@@ -2318,8 +2257,10 @@ def _visit_AnnotatedPyccelSymbol(self, expr):
                 dtype = t.datatype
                 prec  = t.precision
                 rank  = t.rank
-                v = var_class(dtype, name, precision = prec,
-                        shape = None, rank = rank, order = t.order, cls_base = t.cls_base,
+                class_type = t.class_type
+                cls_base = get_cls_base(dtype, prec, class_type) or self.scope.find(class_type.name, 'classes')
+                v = var_class(dtype, name, precision = prec, cls_base = cls_base,
+                        shape = None, rank = rank, order = t.order, class_type = t.class_type,
                         is_const = t.is_const, is_optional = False,
                         memory_handling = array_memory_handling if rank > 0 else 'stack',
                         **kwargs)
@@ -2355,15 +2296,14 @@ def _visit_SyntacticTypeAnnotation(self, expr):
                     types.append(dtype_from_scope)
                 else:
                     types.append(VariableTypeAnnotation(dtype_from_scope.datatype,
-                        dtype_from_scope.cls_base, dtype_from_scope.precision,
+                        dtype_from_scope.class_type, dtype_from_scope.precision,
                         rank, order, dtype_from_scope.is_const))
             elif isinstance(dtype_from_scope, ClassDef):
                 dtype = self.get_class_construct(dtype_name)
                 prec = 0
                 rank = 0
                 order = None
-                cls_base = dtype_from_scope
-                types.append(VariableTypeAnnotation(dtype, cls_base, prec, rank, order))
+                types.append(VariableTypeAnnotation(dtype, dtype, prec, rank, order))
             elif dtype_from_scope is not None:
                 errors.report(PYCCEL_RESTRICTION_TODO + f' Could not deduce type information from {type(dtype_from_scope)} object',
                         severity='fatal', symbol=expr)
@@ -2375,22 +2315,19 @@ def _visit_SyntacticTypeAnnotation(self, expr):
                     errors.report(f'Could not identify type : {dtype_name}',
                             severity='fatal', symbol=expr)
 
-                try:
-                    cls_base = get_cls_base(dtype, prec, rank)
-                except KeyError:
-                    cls_base = self.scope.find(dtype_name, 'classes')
+                # NumPy objects cannot have default precision
+                if prec == -1 and rank > 0:
+                    prec = default_precision[dtype]
 
                 if order is not None and rank < 2:
                     errors.report(f"Ordering is not applicable to objects with rank {rank}",
                             symbol=expr, severity='warning')
                     order = None
 
-                # NumPy objects cannot have default precision
-                if prec == -1 and cls_base is NumpyArrayClass:
-                    prec = default_precision[dtype]
+                cls_type = dtype if rank == 0 else NumpyNDArrayType()
 
                 # Save the potential type
-                types.append(VariableTypeAnnotation(dtype, cls_base, prec, rank, order, is_const))
+                types.append(VariableTypeAnnotation(dtype, cls_type, prec, rank, order, is_const))
 
         # Collect all possible types into a UnionTypeAnnotation
         return UnionTypeAnnotation(*types)
@@ -2462,12 +2399,10 @@ def _visit_DottedName(self, expr):
                     severity='fatal')
 
         d_var = self._infer_type(first)
-        if d_var.get('cls_base', None) is None:
-            errors.report(f'Attribute {rhs_name} not found',
-                bounding_box=(self.current_ast_node.lineno, self.current_ast_node.col_offset),
-                severity='fatal')
-
-        cls_base = d_var['cls_base']
+        dtype = d_var['datatype']
+        cls_base = get_cls_base(dtype, d_var['precision'], d_var['class_type'])
+        if cls_base is None:
+            cls_base = self.scope.find(dtype.name, 'classes')
 
         # look for a class method
         if isinstance(rhs, FunctionCall):
@@ -2526,12 +2461,9 @@ def _visit_PyccelOperator(self, expr):
 
     def _visit_PyccelAdd(self, expr):
         args = [self._visit(a) for a in expr.args]
-        if isinstance(args[0], (TupleVariable, PythonTuple, Concatenate, Duplicate)):
-            is_homogeneous = all((isinstance(a, (TupleVariable, PythonTuple)) and a.is_homogeneous) \
-                                or isinstance(a, (Concatenate, Duplicate)) for a in args)
-            if is_homogeneous:
-                return Concatenate(*args)
-            else:
+        if isinstance(args[0].class_type, (NativeTuple, NativeHomogeneousList)):
+            is_inhomogeneous = any(isinstance(a.class_type, NativeInhomogeneousTuple) for a in args)
+            if is_inhomogeneous:
                 def get_vars(a):
                     if isinstance(a, InhomogeneousTupleVariable):
                         return a.get_vars()
@@ -2548,6 +2480,8 @@ def get_vars(a):
                         raise NotImplementedError(f"Unexpected type {a_type} in tuple addition")
                 tuple_args = [ai for a in args for ai in get_vars(a)]
                 expr_new = PythonTuple(*tuple_args)
+            else:
+                return Concatenate(*args)
         else:
             expr_new = self._create_PyccelOperator(expr, args)
         return expr_new
@@ -2686,6 +2620,8 @@ def _visit_CmathPhase(self, expr):
             return MathAtan2(PythonImag(var), PythonReal(var))
 
     def _visit_Lambda(self, expr):
+        errors.report("Lambda functions are not currently supported",
+                symbol=expr, severity='fatal')
         expr_names = set(str(a) for a in expr.expr.get_attribute_nodes(PyccelSymbol))
         var_names = map(str, expr.variables)
         missing_vars = expr_names.difference(var_names)
@@ -2768,7 +2704,7 @@ def _visit_FunctionCall(self, expr):
             lhs = expr.get_user_nodes(Assign)[0].lhs
             if isinstance(lhs, AnnotatedPyccelSymbol):
                 annotation = self._visit(lhs.annotation)
-                if len(annotation.type_list) != 1 or annotation.type_list[0].cls_base != cls_def:
+                if len(annotation.type_list) != 1 or annotation.type_list[0].class_type.name != method.cls_name:
                     errors.report(f"Unexpected type annotation in creation of {cls_def.name}",
                             symbol=annotation, severity='error')
                 lhs = lhs.name
@@ -3109,7 +3045,7 @@ def _visit_Assign(self, expr):
                             bounding_box=(self.current_ast_node.lineno, self.current_ast_node.col_offset),
                             symbol=li, severity='error')
             else:
-                if l.is_const and (not isinstance(expr.lhs, AnnotatedPyccelSymbol) or len(l.get_all_user_nodes()) > 0):
+                if getattr(l, 'is_const', False) and (not isinstance(expr.lhs, AnnotatedPyccelSymbol) or len(l.get_all_user_nodes()) > 0):
                     # If constant and not the initialising declaration of a constant variable
                     errors.report("Cannot modify 'const' variable",
                         bounding_box=(self.current_ast_node.lineno, self.current_ast_node.col_offset),
@@ -3269,6 +3205,7 @@ def _visit_FunctionalFor(self, expr):
                 if dvar['rank'] == 1:
                     dvar['rank']  = 0
                     dvar['shape'] = None
+                    dvar['class_type'] = dtype
                 if dvar['rank'] > 1:
                     dvar['rank'] -= 1
                     dvar['shape'] = (dvar['shape'])[1:]
diff --git a/pyccel/utilities/metaclasses.py b/pyccel/utilities/metaclasses.py
index c23b631d6e..a02bc81b5b 100644
--- a/pyccel/utilities/metaclasses.py
+++ b/pyccel/utilities/metaclasses.py
@@ -11,20 +11,6 @@
     'ArgumentSingleton',
 )
 
-class Singleton(type):
-    """
-    Metaclass indicating that there is only one instance of the class.
-
-    A metaclass which ensures that only one instance of the class is ever
-    created. Trying to create a second instance will result in accessing
-    the first.
-    """
-    _instances = {}
-    def __call__(cls):
-        if cls not in cls._instances:
-            cls._instances[cls] = super().__call__()
-        return cls._instances[cls]
-
 class ArgumentSingleton(type):
     """ Indicates that there is only one instance of the class
     for any given set of arguments
@@ -35,3 +21,14 @@ def __call__(cls, *args, **kwargs):
         if index not in cls._instances:
             cls._instances[index] = super().__call__(*args, **kwargs)
         return cls._instances[index]
+
+class Singleton(ArgumentSingleton):
+    """
+    Metaclass indicating that there is only one instance of the class.
+
+    A metaclass which ensures that only one instance of the class is ever
+    created. Trying to create a second instance will result in accessing
+    the first.
+    """
+    def __call__(cls):
+        return super().__call__()
diff --git a/pyccel/utilities/stage.py b/pyccel/utilities/stage.py
index 6b769ccc80..c7ff8e3da3 100644
--- a/pyccel/utilities/stage.py
+++ b/pyccel/utilities/stage.py
@@ -49,3 +49,6 @@ def current_stage(self):
         indicating the current stage.
         """
         return self._stage
+
+    def __str__(self):
+        return self._stage


Test patch:
diff --git a/tests/epyccel/recognised_functions/test_cmath_funcs.py b/tests/epyccel/recognised_functions/test_cmath_funcs.py
index 22e1601ba5..c6d1ce49ca 100644
--- a/tests/epyccel/recognised_functions/test_cmath_funcs.py
+++ b/tests/epyccel/recognised_functions/test_cmath_funcs.py
@@ -255,7 +255,8 @@ def tanh_phrase(x : complex, y : complex):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="isfinite not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 @pytest.mark.skipif(os.environ.get('PYCCEL_DEFAULT_COMPILER', None) == 'intel', reason='Nan not correctly passed to intel function')
@@ -285,7 +286,8 @@ def isfinite_call(x : 'T'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="isinf not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 def test_isinf_call(language): # isinf
diff --git a/tests/epyccel/recognised_functions/test_math_funcs.py b/tests/epyccel/recognised_functions/test_math_funcs.py
index b67018118d..f7e138ab91 100644
--- a/tests/epyccel/recognised_functions/test_math_funcs.py
+++ b/tests/epyccel/recognised_functions/test_math_funcs.py
@@ -537,7 +537,8 @@ def copysign_return_type(x : 'float', y : 'int'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="isfinite not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 @pytest.mark.skipif(os.environ.get('PYCCEL_DEFAULT_COMPILER', None) == 'intel', reason='Nan not correctly passed to intel function')
@@ -564,7 +565,8 @@ def isfinite_call(x : 'float'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="isinf not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 def test_isinf_call(language): # isinf
@@ -609,7 +611,8 @@ def isnan_call(x : 'float'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="ldexp not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python),
     )
 )
 def test_ldexp_call(language): # ldexp
@@ -635,7 +638,8 @@ def ldexp_call(x : 'float', exp : 'int'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="ldexp not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 def test_ldexp_return_type(language): # ldexp
@@ -664,7 +668,8 @@ def ldexp_type(x : 'float', exp : 'int'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="remainder not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 def test_remainder_call(language): # remainder
@@ -689,7 +694,8 @@ def remainder_call(x : 'float', y : 'float'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="remainder not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 def test_remainder_return_type(language): # remainder
@@ -755,7 +761,8 @@ def trunc_type(x : 'float'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="expm1 not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 def test_expm1_call(language): # expm1
@@ -772,7 +779,8 @@ def expm1_call(x : 'float'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="expm1 not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 def test_expm1_call_special_case(language): # expm1
@@ -789,7 +797,8 @@ def expm1_call(x : 'float'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="expm1 not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 def test_expm1_phrase(language): # expm1
@@ -808,7 +817,8 @@ def expm1_phrase(x : 'float', y : 'float'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="expm1 not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 def test_expm1_return_type(language): # expm1 # expm1
@@ -829,7 +839,8 @@ def expm1_type(x : 'float'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="log1p not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 def test_log1p_call(language):
@@ -847,7 +858,8 @@ def log1p_call(x : 'float'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="log1p not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 def test_log1p_phrase(language):
@@ -867,7 +879,8 @@ def log1p_phrase(x : 'float', y : 'float'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="log2 not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 def test_log2_call(language):
@@ -887,7 +900,8 @@ def log2_call(x : 'float'):
         pytest.param("fortran", marks = [
             pytest.mark.xfail(reason="log2 not implemented"),
             pytest.mark.fortran]
-        )
+        ),
+        pytest.param("python", marks = pytest.mark.python)
     )
 )
 def test_log2_phrase(language):
diff --git a/tests/epyccel/test_arrays.py b/tests/epyccel/test_arrays.py
index 48fa57f8bc..39ec376a90 100644
--- a/tests/epyccel/test_arrays.py
+++ b/tests/epyccel/test_arrays.py
@@ -302,6 +302,15 @@ def test_array_int_1d_initialization_1(language):
 
     assert f1() == f2()
 
+@pytest.mark.parametrize( 'language', (
+        pytest.param("fortran", marks = pytest.mark.fortran),
+        pytest.param("c", marks = [
+            pytest.mark.skip(reason="Lists not yet supported"),
+            pytest.mark.c]
+        ),
+        pytest.param("python", marks = pytest.mark.python)
+    )
+)
 def test_array_int_1d_initialization_2(language):
 
     f1 = arrays.array_int_1d_initialization_2
diff --git a/tests/epyccel/test_epyccel_IfTernaryOperator.py b/tests/epyccel/test_epyccel_IfTernaryOperator.py
index 29c466a90c..1b55a62781 100644
--- a/tests/epyccel/test_epyccel_IfTernaryOperator.py
+++ b/tests/epyccel/test_epyccel_IfTernaryOperator.py
@@ -117,6 +117,15 @@ def f6wp(x : 'int'):
     # ...
 #------------------------------------------------------------------------------
 
+@pytest.mark.parametrize( 'language', (
+        pytest.param("fortran", marks = pytest.mark.fortran),
+        pytest.param("c", marks = [
+            pytest.mark.skip(reason="Lists not yet supported"),
+            pytest.mark.c]
+        ),
+        pytest.param("python", marks = pytest.mark.python)
+    )
+)
 def test_f7(language):
     def f7(x : 'int'):
         a = [1.,2.,3.] if x < 5 else [1.5,6.5,7.5]
@@ -227,6 +236,15 @@ def f11(x : 'int'):
     # ...
 #------------------------------------------------------------------------------
 
+@pytest.mark.parametrize( 'language', (
+        pytest.param("fortran", marks = pytest.mark.fortran),
+        pytest.param("c", marks = [
+            pytest.mark.skip(reason="Lists not yet supported"),
+            pytest.mark.c]
+        ),
+        pytest.param("python", marks = pytest.mark.python)
+    )
+)
 def test_f12(language):
     def f12(x : 'int'):
         a = [1.,2.,3.,4.] if x < 5 else [1.5,6.5,7.5]
diff --git a/tests/errors/semantic/blocking/LIST_OF_TUPLES.py b/tests/errors/semantic/non_blocking/LIST_OF_TUPLES.py
similarity index 100%
rename from tests/errors/semantic/blocking/LIST_OF_TUPLES.py
rename to tests/errors/semantic/non_blocking/LIST_OF_TUPLES.py