diff --git a/devito/ir/cgen/printer.py b/devito/ir/cgen/printer.py
index 27ee0e3112..f7fcae920b 100644
--- a/devito/ir/cgen/printer.py
+++ b/devito/ir/cgen/printer.py
@@ -280,8 +280,20 @@ def _print_Mod(self, expr):
         return '%'.join(args)
 
     def _print_Mul(self, expr):
-        args = [a for a in expr.args if a != -1]
-        neg = (len(expr.args) - len(args)) % 2
+        args, neg = [], False
+        # SymPy's CodePrinter extracts negative coefficients with the standalone
+        # _keep_coeff helper, which can rebuild a plain SymPy Mul. Since this
+        # is not a method, UnevalMul cannot override it. Extract signs before
+        # delegating to preserve the product's class without patching SymPy or
+        # copying its multiplication printer.
+        for a in expr.args:
+            if a.is_Number and a.is_negative:
+                a = -a
+                neg = not neg
+                if a == 1:
+                    continue
+            args.append(a)
+        args = args or [S.One]
 
         if len(args) > 1:
             term = super()._print_Mul(expr.func(*args, evaluate=False))
@@ -289,7 +301,7 @@ def _print_Mul(self, expr):
             term = self.parenthesize(args[0], precedence(expr))
 
         if neg:
-            return f'-{term}'
+            return term[1:] if term.startswith('-') else f'-{term}'
         else:
             return term
 
diff --git a/devito/ir/iet/visitors.py b/devito/ir/iet/visitors.py
index 6186ff4e76..8a825c7619 100644
--- a/devito/ir/iet/visitors.py
+++ b/devito/ir/iet/visitors.py
@@ -1506,7 +1506,8 @@ def visit_PointerCast(self, o):
     def visit_Dereference(self, o):
         pointee = self.mapper.get(o.pointee, o.pointee)
         pointer = self.mapper.get(o.pointer, o.pointer)
-        return reuse_if_unchanged(o, pointee=pointee, pointer=pointer)
+        offset = uxreplace(o.offset, self.mapper)
+        return reuse_if_unchanged(o, pointee=pointee, pointer=pointer, offset=offset)
 
     def visit_Pragma(self, o):
         arguments = [uxreplace(i, self.mapper) for i in o.arguments]
diff --git a/devito/ir/support/symregistry.py b/devito/ir/support/symregistry.py
index d7fca20f7c..e32e37acd6 100644
--- a/devito/ir/support/symregistry.py
+++ b/devito/ir/support/symregistry.py
@@ -1,4 +1,5 @@
-from devito.tools import generator
+from devito.symbolics import retrieve_functions
+from devito.tools import as_tuple, generator
 from devito.types import (
     DeviceID, NPThreads, NThreads, NThreadsNested, NThreadsNonaffine, ThreadID
 )
@@ -15,6 +16,7 @@ class SymbolRegistry:
     def __init__(self):
         # {name -> generator()} -- to create unique names for symbols, functions, ...
         self.counters = {}
+        self.reserved = set()
 
         # Special symbols
         self.nthreads = NThreads()
@@ -31,6 +33,14 @@ def __init__(self):
         # passes, to maximize symbol (especially Dimension) reuse
         self.caches = {}
 
+    def reserve(self, exprs):
+        """Reserve identifiers already used by the input expressions."""
+        exprs = as_tuple(exprs)
+        self.reserved.update(i.name for e in exprs for i in e.free_symbols)
+
+        # Before indexification, a Function's free symbols omit its own name
+        self.reserved.update(i.name for i in retrieve_functions(exprs, mode='unique'))
+
     def make_name(self, prefix=None, increment_first=True):
         # By default we're creating a new symbol
         if prefix is None:
@@ -42,10 +52,15 @@ def make_name(self, prefix=None, increment_first=True):
             counter = self.counters.setdefault(prefix, generator())
 
             # Only increment symbol names after the first encountered
-            if not increment_first:
+            if not increment_first and prefix not in self.reserved:
+                self.reserved.add(prefix)
                 return prefix
 
-        return f'{prefix}{counter()}'
+        name = f'{prefix}{counter()}'
+        while name in self.reserved:
+            name = f'{prefix}{counter()}'
+        self.reserved.add(name)
+        return name
 
     def make_npthreads(self, size):
         name = self.make_name(prefix='npthreads')
diff --git a/devito/operator/operator.py b/devito/operator/operator.py
index 8f3c02f983..e36c48a351 100644
--- a/devito/operator/operator.py
+++ b/devito/operator/operator.py
@@ -293,6 +293,8 @@ def _lower(cls, expressions, **kwargs):
 
         expressions = as_tuple(expressions)
 
+        kwargs['sregistry'].reserve(expressions)
+
         # Enable recursive lowering
         # This may be used by a compilation pass that constructs a new
         # expression for which a partial or complete lowering is desired
@@ -472,6 +474,11 @@ def _specialize_iet(cls, graph, **kwargs):
         """
         return graph
 
+    @classmethod
+    def _finalize_iet(cls, graph, **kwargs):
+        """Backend hook after lowering, before macro and signature extraction."""
+        return graph
+
     @classmethod
     @timed_pass(name='lowering.uIET')
     def _lower_uiet(cls, stree, profiler=None, **kwargs):
@@ -515,15 +522,17 @@ def _lower_iet(cls, uiet, **kwargs):
         # specialization further Sections may be introduced
         cls._Target.instrument(graph, **kwargs)
 
-        # Extract the necessary macros from the symbolic objects
-        generate_macros(graph, **kwargs)
-
         # Target-specific lowering
         lower_dtypes(graph, **kwargs)
 
         # Target-independent optimizations
         minimize_symbols(graph)
 
+        graph = cls._finalize_iet(graph, **kwargs)
+
+        # Extract macros only after the symbolic objects have their final names
+        generate_macros(graph, **kwargs)
+
         # Finalize helper signatures after all IET transformations have settled.
         finalize_args(graph)
 
diff --git a/devito/passes/clusters/cse.py b/devito/passes/clusters/cse.py
index 9c899eb23e..f00bf05e04 100644
--- a/devito/passes/clusters/cse.py
+++ b/devito/passes/clusters/cse.py
@@ -293,15 +293,14 @@ def _toposort(exprs):
             if e0.lhs in e1.rhs.free_symbols:
                 dag.add_edge(e0, e1, force_add=True)
 
+    # Cache positions to avoid repeated linear scans while choosing ready nodes
+    ordering = {e: n for n, e in enumerate(exprs)}
+
     def choose_element(queue, scheduled):
+        # Preserve expression order among ready temporaries, independent of names
         tmps = [i for i in queue if isinstance(i.lhs, CTemp)]
-        if tmps:
-            # Try to honor temporary names as much as possible
-            first = sorted(tmps, key=lambda i: i.lhs.name).pop(0)
-            queue.remove(first)
-        else:
-            first = sorted(queue, key=lambda i: exprs.index(i)).pop(0)
-            queue.remove(first)
+        first = min(tmps or queue, key=ordering.__getitem__)
+        queue.remove(first)
         return first
 
     processed = dag.topological_sort(choose_element)
diff --git a/devito/passes/iet/engine.py b/devito/passes/iet/engine.py
index fe9a78b65b..8615582da5 100644
--- a/devito/passes/iet/engine.py
+++ b/devito/passes/iet/engine.py
@@ -564,6 +564,8 @@ def abstract_objects(objects0):
     mapper = {}
     sregistry = SymbolRegistry()
     for i in objects:
+        if i in mapper:
+            continue
         abstract_object(i, mapper, sregistry)
 
     return mapper
diff --git a/devito/passes/iet/misc.py b/devito/passes/iet/misc.py
index 1b404daca1..a94a754124 100644
--- a/devito/passes/iet/misc.py
+++ b/devito/passes/iet/misc.py
@@ -188,8 +188,9 @@ def _generate_macros_findexeds(iet, sregistry=None, tracker=None, **kwargs):
 
     subs = {}
     for i in findexeds:
+        key = (i.base, i.strides_map)
         try:
-            v = tracker[i.base].v
+            v = tracker[key].v
             subs[i] = v.func(v.base, *i.indices)
             continue
         except KeyError:
@@ -199,7 +200,7 @@ def _generate_macros_findexeds(iet, sregistry=None, tracker=None, **kwargs):
         header, v = i.bind(pname)
 
         subs[i] = v
-        tracker[i.base] = Bunch(header=header, v=v)
+        tracker[key] = Bunch(header=header, v=v)
 
     iet = Uxreplace(subs).visit(iet)
 
diff --git a/devito/symbolics/manipulation.py b/devito/symbolics/manipulation.py
index b545290a15..11e085722a 100644
--- a/devito/symbolics/manipulation.py
+++ b/devito/symbolics/manipulation.py
@@ -3,7 +3,7 @@
 from functools import singledispatch
 
 import numpy as np
-from sympy import Add, Max, Min, Mul, Pow, S, SympifyError, Tuple, sympify
+from sympy import Add, Max, Min, Mul, Pow, S, SympifyError, Tuple, ceiling, floor, sympify
 from sympy.core import Basic as SympyBasic
 from sympy.core.add import _addsort
 from sympy.core.mul import _mulsort
@@ -12,7 +12,7 @@
     EvalDerivative, IndexDerivative, LocalSum
 )
 from devito.symbolics.extended_dtypes import LONG
-from devito.symbolics.extended_sympy import DefFunction, rfunc
+from devito.symbolics.extended_sympy import CallFromPointer, DefFunction, rfunc
 from devito.symbolics.queries import q_leaf
 from devito.symbolics.search import (
     retrieve_functions, retrieve_indexed, retrieve_terminals
@@ -21,12 +21,15 @@
 from devito.symbolics.unevaluation import Mul as UnevalMul
 from devito.symbolics.unevaluation import Pow as UnevalPow
 from devito.symbolics.unevaluation import UnevaluableMixin
+from devito.symbolics.unevaluation import ceiling as UnevalCeil
+from devito.symbolics.unevaluation import floor as UnevalFloor
 from devito.tools import (
-    EnrichedTuple, as_list, as_tuple, flatten, split, transitive_closure
+    EnrichedTuple, as_list, as_tuple, flatten, frozendict, split, transitive_closure
 )
 from devito.types.array import ComponentAccess
 from devito.types.basic import Basic, Indexed
 from devito.types.equation import Eq
+from devito.types.misc import FIndexed
 from devito.types.relational import Ge, Gt, Le, Lt
 
 __all__ = [
@@ -53,8 +56,9 @@ def uxreplace(expr, rule):
     be avoided at all costs (e.g., to prevent SymPy from unpicking Devito
     transformations, such as factorization).
 
-    The canonical ordering of the arguments is however guaranteed; where this is
-    not possible, a re-evaluation will be enforced.
+    Evaluated expressions retain canonical argument ordering. Lowered,
+    unevaluated arithmetic retains its operand order and grouping unless
+    substitution makes it completely numeric.
 
     By avoiding re-evaluations, this function is typically much quicker than
     SymPy's xreplace.
@@ -77,6 +81,10 @@ def uxreplace(expr, rule):
 
 @singledispatch
 def _uxreplace(expr, rule):
+    return _uxreplace_default(expr, rule)
+
+
+def _uxreplace_default(expr, rule):
     if expr in rule:
         v = rule[expr]
         if not isinstance(v, dict):
@@ -117,6 +125,17 @@ def _uxreplace(expr, rule):
     return expr, False
 
 
+@_uxreplace.register(CallFromPointer)
+def _(expr, rule):
+    if expr in rule or expr.__class__ in rule:
+        return _uxreplace_default(expr, rule)
+
+    args, changed = _uxreplace_dispatch((expr.pointer, expr.params), rule)
+    if changed:
+        return expr.func(expr.call, *args), True
+    return expr, False
+
+
 @_uxreplace.register(LocalSum)
 def _(expr, rule):
     # Only the sum owns these guards; ordinary ConditionalDimensions stay atomic
@@ -128,7 +147,7 @@ def _(expr, rule):
             mapper[d] = d._rebuild(**kwargs)
 
     # Apply the same dimension substitutions to the summand and its bindings
-    return _uxreplace.dispatch(object)(expr, {**mapper, **rule})
+    return _uxreplace_default(expr, {**mapper, **rule})
 
 
 @singledispatch
@@ -165,6 +184,7 @@ def _(iterable, rule):
 
 
 @_uxreplace_dispatch.register(dict)
+@_uxreplace_dispatch.register(frozendict)
 def _(mapper, rule):
     ret = {}
     changed = False
@@ -172,7 +192,7 @@ def _(mapper, rule):
         vx, flag = _uxreplace_dispatch(v, rule)
         ret[k] = vx
         changed |= flag
-    return (ret, True) if changed else (mapper, False)
+    return (mapper.__class__(ret), True) if changed else (mapper, False)
 
 
 @singledispatch
@@ -183,6 +203,12 @@ def _uxreplace_handle(expr, args, kwargs):
         return expr.func(*args)
 
 
+@_uxreplace_handle.register(CallFromPointer)
+def _(expr, args, kwargs):
+    # C member names belong to their compound type, rather than the local scope
+    return expr.func(expr.call, *args[1:])
+
+
 @_uxreplace_handle.register(Min)
 @_uxreplace_handle.register(Max)
 @_uxreplace_handle.register(Pow)
@@ -191,6 +217,16 @@ def _(expr, args, kwargs):
     return expr.func(*args, evaluate=evaluate)
 
 
+# The SymPy bases precede UnevaluableMixin in these classes' dispatch MRO
+@_uxreplace_handle.register(UnevalAdd)
+@_uxreplace_handle.register(UnevalMul)
+def _(expr, args, kwargs):
+    # Lowered arithmetic must survive substitutions without reassociation
+    if all(i.is_Number for i in args):
+        return expr.func.__base__(*args)
+    return expr.func(*args)
+
+
 @_uxreplace_handle.register(Add)
 def _(expr, args, kwargs):
     args = [i for i in args if i != 0]
@@ -252,6 +288,7 @@ def dispatchable(self, obj):
 _uxreplace_registry.register(Eq)
 _uxreplace_registry.register(DefFunction)
 _uxreplace_registry.register(ComponentAccess)
+_uxreplace_registry.register(FIndexed)
 _uxreplace_registry.register(LocalSum)
 
 
@@ -585,10 +622,14 @@ def evalrel(func=min, input=None, assumptions=None):
     return rfunc(func, *input)
 
 
-uneval_mapper = {Add: UnevalAdd, Mul: UnevalMul, Pow: UnevalPow}
+# Rounding functions also need freezing: reconstruction otherwise reevaluates
+# their lowered arguments, including during Operator unpickling.
+uneval_mapper = {Add: UnevalAdd, Mul: UnevalMul, Pow: UnevalPow,
+                 ceiling: UnevalCeil, floor: UnevalFloor}
 
 
 def unevaluate(expr):
+    """Protect lowered arithmetic from reevaluation during reconstruction."""
     if q_leaf(expr):
         return expr
 
diff --git a/devito/symbolics/unevaluation.py b/devito/symbolics/unevaluation.py
index 6a860c70b5..7e348ec518 100644
--- a/devito/symbolics/unevaluation.py
+++ b/devito/symbolics/unevaluation.py
@@ -1,7 +1,7 @@
 import sympy
 from sympy import sympify
 
-__all__ = ['Add', 'Mod', 'Mul', 'Pow']
+__all__ = ['Add', 'Mod', 'Mul', 'Pow', 'ceiling', 'floor']
 
 
 class UnevaluableMixin:
@@ -38,3 +38,11 @@ def __new__(cls, base, exp, evaluate=None, **kwargs):
             return sympy.S.One
         else:
             return cls.__base__.__new__(cls, base, exp, evaluate=False, **kwargs)
+
+
+class ceiling(sympy.ceiling, UnevaluableMixin):
+    __new__ = UnevaluableMixin.__new__
+
+
+class floor(sympy.floor, UnevaluableMixin):
+    __new__ = UnevaluableMixin.__new__
diff --git a/examples/performance/00_overview.ipynb b/examples/performance/00_overview.ipynb
index 917b75df65..3345a1ac4e 100644
--- a/examples/performance/00_overview.ipynb
+++ b/examples/performance/00_overview.ipynb
@@ -208,7 +208,7 @@
       "      {\n",
       "        for (int z = z_m; z <= z_M; z += 1)\n",
       "        {\n",
-      "          u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
+      "          u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
       "        }\n",
       "      }\n",
       "    }\n",
@@ -283,7 +283,7 @@
       "      {\n",
       "        for (int z = z_m; z <= z_M; z += 1)\n",
       "        {\n",
-      "          u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
+      "          u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
       "        }\n",
       "      }\n",
       "    }\n",
@@ -341,7 +341,7 @@
       "          {\n",
       "            for (int z = z_m; z <= z_M; z += 1)\n",
       "            {\n",
-      "              u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
+      "              u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
       "            }\n",
       "          }\n",
       "        }\n",
@@ -431,7 +431,7 @@
       "                  {\n",
       "                    for (int z = z0_blk1; z <= MIN(z_M, z0_blk1 + z0_blk1_size - 1); z += 1)\n",
       "                    {\n",
-      "                      u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
+      "                      u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
       "                    }\n",
       "                  }\n",
       "                }\n",
@@ -567,7 +567,7 @@
       "         {\n",
       "           for (int z = z_m; z <= z_M; z += 1)\n",
       "           {\n",
-      "-            u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
+      "-            u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
       "+            float r0 = 1.0F/h_y;\n",
       "+            u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 5][z + 4]) + (-8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 8][z + 4]) + (8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 7][z + 4]))*sinf(f[x + 1][y + 1][z + 1]);\n",
       "           }\n",
@@ -882,7 +882,7 @@
       "      {\n",
       "        for (int z = z_m; z <= z_M; z += 1)\n",
       "        {\n",
-      "          u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
+      "          u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
       "        }\n",
       "      }\n",
       "    }\n",
@@ -1020,7 +1020,7 @@
       "      {\n",
       "        for (int z = z_m; z <= z_M; z += 1)\n",
       "        {\n",
-      "          u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + ((-6.66666667e-1F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + ((-8.33333333e-2F)/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
+      "          u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 1][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 3][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y) + (6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 3][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 6][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 7][z + 4]/h_y) + (-6.66666667e-1F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 1][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 2][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 4][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 5][z + 4]/h_y) + (-8.33333333e-2F/h_y)*(8.33333333e-2F*u[t0][x + 4][y + 4][z + 4]/h_y - 6.66666667e-1F*u[t0][x + 4][y + 5][z + 4]/h_y + 6.66666667e-1F*u[t0][x + 4][y + 7][z + 4]/h_y - 8.33333333e-2F*u[t0][x + 4][y + 8][z + 4]/h_y))*sinf(f[x + 1][y + 1][z + 1]);\n",
       "        }\n",
       "      }\n",
       "    }\n",
diff --git a/tests/test_cse.py b/tests/test_cse.py
index 6ab68f05a0..4fc439c49b 100644
--- a/tests/test_cse.py
+++ b/tests/test_cse.py
@@ -8,7 +8,7 @@
     ConditionalDimension, Eq, Function, Grid, Operator, TimeFunction, cos, sin
 )
 from devito.finite_differences.differentiable import diffify
-from devito.ir import Conditional, DummyEq, FindNodes, FindSymbols
+from devito.ir import Conditional, DummyEq, Expression, FindNodes, FindSymbols
 from devito.ir.support import generator
 from devito.passes.clusters.cse import CTemp, _cse
 from devito.symbolics import DefFunction, FieldFromPointer, IndexedPointer, indexify
@@ -131,6 +131,23 @@ def test_temp_order():
     assert type(args[2]) is CTemp
 
 
+@pytest.mark.parametrize('names', [('r2', 'r3'), ('r9', 'r10'), ('later', 'earlier')])
+@pytest.mark.parametrize('mode', ['smartsort', 'advanced'])
+def test_schedule_independent_of_temp_names(names, mode):
+    grid = Grid(shape=(4, 4))
+    f = Function(name='f', grid=grid)
+    g = Function(name='g', grid=grid)
+    u = Function(name='u', grid=grid)
+    v = Function(name='v', grid=grid)
+    a, b = [CTemp(name=name, dtype=np.float32) for name in names]
+
+    op = Operator([Eq(a, cos(f)), Eq(b, sin(g)), Eq(u, a + b), Eq(v, a*b)],
+                  opt=('advanced', {'cse-algo': mode}))
+
+    exprs = [e.expr for e in FindNodes(Expression).visit(op)]
+    assert [e.lhs for e in exprs if isinstance(e.lhs, CTemp)][:2] == [a, b]
+
+
 def test_deffunction_not_captured():
     """
     Opaque C-level calls (DefFunction) must not be CSE-captured: their return
diff --git a/tests/test_iet.py b/tests/test_iet.py
index 99a94ed24c..2834ac2712 100644
--- a/tests/test_iet.py
+++ b/tests/test_iet.py
@@ -12,7 +12,7 @@
 from devito.ir.iet import (
     Call, Callable, CGen, Conditional, Definition, Dereference, DeviceCall, DummyExpr,
     ElementalFunction, FindNodes, FindSymbols, Iteration, KernelLaunch, Lambda, List,
-    Switch, Transformer, filter_iterations, make_callable, make_efunc,
+    Switch, Transformer, Uxreplace, filter_iterations, make_callable, make_efunc,
     retrieve_iteration_tree
 )
 from devito.ir.iet.visitors import sorted_efuncs
@@ -120,6 +120,20 @@ def test_nested_calls_cgen():
     assert str(code) == 'foo(bar());'
 
 
+def test_uxreplace_dereference_offset():
+    d = CustomDimension(name='d', symbolic_size=4)
+    array = Array(name='a0', dimensions=(d,), shape=(4,), scope='shared')
+    pointer = Pointer(name='sm7')
+    old, new = Symbol(name='off7'), Symbol(name='off0')
+    node = Dereference(array, pointer, offset=old + 32)
+
+    replaced = Uxreplace({old: new}).visit(node)
+
+    assert replaced.offset == new + 32
+    assert replaced.pointee is array
+    assert replaced.pointer is pointer
+
+
 @pytest.mark.parametrize('mode,expected', [
     ('basics', '["x"]'),
     ('symbolics', '["f"]')
@@ -542,14 +556,17 @@ def test_codegen_quality0():
     assert foo1.parameters[0] is a
 
 
-def test_reuse_efuncs_natural_numbering():
+@pytest.mark.parametrize('copy_symbols', [False, True])
+def test_reuse_efuncs_natural_numbering(copy_symbols):
 
-    def make_foo(name, start):
+    def make_foo(name, start, copies=False):
         pointers = [Pointer(name=f'p{i}') for i in range(start, start + 3)]
-        return make_callable(name, Call('bar', pointers))
+        aliases = [p._rebuild() for p in pointers] if copies else pointers
+        return make_callable(name, List(body=[Call('bar', pointers),
+                                              Call('bar', aliases)]))
 
     foo0 = make_foo('foo0', 0)
-    foo1 = make_foo('foo1', 9)
+    foo1 = make_foo('foo1', 9, copies=copy_symbols)
 
     root = make_callable('root', List(body=[
         Call(foo0.name, foo0.parameters),
diff --git a/tests/test_operator.py b/tests/test_operator.py
index e6fc86463a..a2cf17861a 100644
--- a/tests/test_operator.py
+++ b/tests/test_operator.py
@@ -198,6 +198,19 @@ def test_jit_compiler_flags(self, envvar, flags, monkeypatch):
         assert op0._soname != op1._soname
         assert op1._soname == op2._soname
 
+    @pytest.mark.parametrize('name', ['r0', 'r1'])
+    def test_generated_names_avoid_inputs(self, name):
+        grid = Grid(shape=(8, 8), extent=(7., 7.))
+        f = Function(name=name, grid=grid, space_order=4)
+        u = Function(name='u', grid=grid)
+        x, y = np.indices(f.data_with_halo.shape, dtype=np.float32)
+        f.data_with_halo[:] = x + 2*y
+
+        op = Operator(Eq(u, f.dx + f.dy))
+        op.apply()
+
+        np.testing.assert_allclose(u.data, 3., rtol=1.e-5, atol=1.e-6)
+
 
 class TestCodeGen:
 
diff --git a/tests/test_pickle.py b/tests/test_pickle.py
index 8bd989d31f..838dd2ea80 100644
--- a/tests/test_pickle.py
+++ b/tests/test_pickle.py
@@ -9,7 +9,7 @@
 from devito import (
     MPI, ConditionalDimension, Constant, Dimension, Eq, Function, Grid, IncrDimension,
     Min, Operator, PrecomputedSparseTimeFunction, SparseFunction, SteppingDimension,
-    SubDimension, SubDomain, TimeDimension, TimeFunction, floor, solve
+    SubDimension, SubDomain, TimeDimension, TimeFunction, ceiling, floor, solve
 )
 from devito.data import LEFT, OWNED
 from devito.finite_differences.tools import centered, direct, left, right, transpose
@@ -905,14 +905,15 @@ def test_collected_coeffs(self, pickle):
 
         assert str(op) == str(new_op)
 
-    def test_elemental(self, pickle):
+    @pytest.mark.parametrize('func', [ceiling, floor])
+    def test_elemental(self, pickle, func):
         """
         Tests that elementary functions don't get reconstructed differently.
         """
         grid = Grid(shape=(4, 4))
         f = Function(name='f', grid=grid)
 
-        op = Operator(Eq(f, floor(f + 0.5)))
+        op = Operator(Eq(f, func(f + 0.5)))
 
         pkl_op = pickle.dumps(op)
         new_op = pickle.loads(pkl_op)
diff --git a/tests/test_symbolics.py b/tests/test_symbolics.py
index 51a759a705..7681fef792 100644
--- a/tests/test_symbolics.py
+++ b/tests/test_symbolics.py
@@ -15,17 +15,20 @@
 from devito.ir.support.guards import GuardExpr, pairwise_or, simplify_and
 from devito.mpi.halo_scheme import HaloTouch
 from devito.symbolics import (  # noqa
-    INT, BaseCast, CallFromPointer, Cast, DefFunction, FieldFromComposite,
-    FieldFromPointer, IntDiv, ListInitializer, Namespace, ReservedWord, RoundUp, Rvalue,
-    SizeOf, VectorAccess, evalrel, pow_to_mul, q_negative, q_positive,
-    retrieve_derivatives, retrieve_functions, retrieve_indexed, subs_if_composite,
-    uxreplace, xreplace_indices
+    INT, BaseCast, CallFromComposite, CallFromPointer, Cast, DefFunction,
+    FieldFromComposite, FieldFromPointer, IntDiv, ListInitializer, Namespace,
+    ReservedWord, RoundUp, Rvalue, SizeOf, VectorAccess, evalrel, pow_to_mul, q_negative,
+    q_positive, retrieve_derivatives, retrieve_functions, retrieve_indexed,
+    subs_if_composite, unevaluate, uxreplace, xreplace_indices
 )
-from devito.tools import CustomDtype, as_tuple, dtypes_vector_mapper
+from devito.symbolics.unevaluation import Add as UnevalAdd
+from devito.symbolics.unevaluation import Mul as UnevalMul
+from devito.tools import CustomDtype, as_tuple, dtypes_vector_mapper, frozendict
 from devito.types import (
     Array, Bundle, ComponentAccess, FIndexed, LocalObject, Object, StencilDimension
 )
 from devito.types import Symbol as dSymbol
+from devito.types import Temp
 from devito.types.basic import AbstractSymbol
 
 
@@ -336,6 +339,123 @@ def test_extended_sympy_arithmetic():
     assert ccode(-1 + bar) == '-1 + o->bar'
 
 
+@pytest.mark.parametrize('cls, regular, separator', [
+    (UnevalAdd, sympy.Add, ' + '),
+    (UnevalMul, sympy.Mul, '*')
+])
+def test_lowered_arithmetic_printing(cls, regular, separator):
+    a, b, c = [dSymbol(i, dtype=np.float32) for i in ('a', 'b', 'c')]
+
+    # Lowering retains ordinary canonical printing
+    assert ccode(cls(c, a, b)) == separator.join(('a', 'b', 'c'))
+    renamed = uxreplace(cls(c, a, b), {c: a, a: c})
+    assert ccode(renamed) == separator.join(('a', 'b', 'c'))
+
+    # Ordinary index arithmetic retains its usual printing inside lowered math
+    grid = Grid(shape=(8,))
+    x, = grid.dimensions
+    f = Function(name='f', grid=grid)
+    assert ccode(cls(c, f.indexed[x + 1], a)) == \
+        ccode(regular(c, f.indexed[x + 1], a, evaluate=False))
+    assert ccode(x + 1) == 'x + 1'
+
+    # Nested sums retain main's printing policy
+    assert ccode(UnevalMul(b, UnevalAdd(c, a))) == 'b*(a + c)'
+    assert ccode(cls(b, cls(c, a))) == \
+        ccode(regular(b, regular(c, a, evaluate=False), evaluate=False))
+    assert ccode(UnevalAdd(b, sympy.Add(c, a, evaluate=False))) == 'b + a + c'
+
+
+@pytest.mark.parametrize('coefficient', [-2, sympy.Float(-2), sympy.Rational(-1, 2)])
+@pytest.mark.parametrize('position', [0, 1, 3])
+def test_lowered_negative_product(coefficient, position):
+    a, b, c = [dSymbol(i, dtype=np.float32) for i in ('a', 'b', 'c')]
+    args = [c, a, b]
+    args.insert(position, coefficient)
+
+    assert ccode(UnevalMul(*args)) == ccode(coefficient*a*b*c)
+    expr = coefficient*c*a*b
+    assert ccode(unevaluate(expr)) == ccode(expr)
+
+
+@pytest.mark.parametrize('coefficient', [-1, -2, sympy.Float(-2), sympy.Rational(-1, 2)])
+def test_lowered_product_sign(coefficient):
+    a, b, c = [dSymbol(i, dtype=np.float32) for i in ('a', 'b', 'c')]
+    expr = UnevalMul(-1, coefficient, c, a, b)
+
+    assert ccode(expr) == ccode(-coefficient*a*b*c)
+
+
+def test_lowered_sign_identity():
+    assert ccode(UnevalMul(-1, -1)) == '1'
+
+
+@pytest.mark.parametrize('cls, lowered', [
+    (sympy.Add, UnevalAdd),
+    (sympy.Mul, UnevalMul)
+])
+def test_freeze_grouping(cls, lowered):
+    a, b, c = [Temp(name=n, dtype=np.float32) for n in ('r9', 'r10', 'r11')]
+    expr = cls(c, cls(b, a, evaluate=False), evaluate=False)
+
+    frozen = unevaluate(expr)
+
+    assert frozen.args == (c, lowered(b, a))
+    assert ccode(frozen) == ccode(expr)
+    assert unevaluate(frozen) is frozen
+
+
+def test_freeze_equal_instances():
+    a = Symbol.__xnew__(Symbol, 'equal_frozen_operand')
+    b = Symbol.__xnew__(Symbol, 'equal_frozen_operand')
+    expr = sympy.Add(a, b, evaluate=False)
+
+    frozen = unevaluate(expr)
+
+    assert a == b and a is not b
+    assert frozen.args[0] is a
+    assert frozen.args[1] is b
+
+
+@pytest.mark.parametrize('dtype', [np.int32, np.int64])
+@pytest.mark.parametrize('offset, expected', [(-1, ' - 1'), (1, ' + 1')])
+def test_lowered_integer_offset(dtype, offset, expected):
+    size = dSymbol('bx_size', dtype=dtype)
+    other = dSymbol('a_size', dtype=dtype)
+
+    expr = size + offset
+    lowered = unevaluate(expr)
+    assert lowered.args == expr.args
+    assert ccode(lowered) == 'bx_size' + expected
+
+    # Renaming preserves the tree, while printing remains canonical
+    lowered = unevaluate(size + other + offset)
+    renamed = uxreplace(lowered, {size: other, other: size})
+    assert ccode(renamed) == 'a_size + bx_size' + expected
+
+    # Result types do not affect how the arithmetic order is frozen
+    value = dSymbol('value', dtype=np.float32)
+    for addend in (value, sympy.sin(size)):
+        expr = size + addend + offset
+        assert ccode(unevaluate(expr)) == ccode(expr)
+
+
+@pytest.mark.parametrize('cls', [sympy.Add, sympy.Mul])
+@pytest.mark.parametrize('names', [('r9', 'r10', 'r11'), ('r21', 'r20', 'r19'),
+                                   ('later', 'middle', 'earlier')])
+def test_noop_temporary_order(cls, names):
+    grid = Grid(shape=(4,))
+    f, g, h, out = [Function(name=n, grid=grid) for n in ('f', 'g', 'h', 'out')]
+    a, b, c = [Temp(name=n, dtype=np.float32) for n in names]
+    rhs = cls(a, b, c)
+    op = Operator([Eq(a, cos(f)), Eq(b, sin(g)), Eq(c, cos(h)),
+                   Eq(out, rhs)], opt='noop')
+    expr, = [i for i in FindNodes(Expression).visit(op) if i.write == out]
+
+    assert tuple(i.name for i in expr.expr.rhs.args) == tuple(i.name for i in rhs.args)
+    assert ccode(expr.expr.rhs) == ccode(rhs)
+
+
 def test_integer_abs():
     i1 = Dimension(name="i1")
     assert ccode(Abs(i1 - 1)) == "abs(i1 - 1)"
@@ -1037,17 +1157,20 @@ def test_xreplace_indices_reuses_untouched_sequence(self):
         assert xreplace_indices(exprs, {z: z + 1}) is exprs
         assert xreplace_indices(exprs, {x: x + 1}) is not exprs
 
-    def test_custom_reconstructable(self):
+    @pytest.mark.parametrize('mapping_type', [dict, frozendict])
+    def test_custom_reconstructable(self, mapping_type):
 
         class MyDefFunction(DefFunction):
             __rargs__ = ('name', 'arguments')
-            __rkwargs__ = ('p0', 'p1', 'p2')
+            __rkwargs__ = ('p0', 'p1', 'p2', 'p3')
 
-            def __new__(cls, name=None, arguments=None, p0=None, p1=None, p2=None):
+            def __new__(cls, name=None, arguments=None, p0=None, p1=None, p2=None,
+                        p3=None):
                 obj = super().__new__(cls, name=name, arguments=arguments)
                 obj.p0 = p0
                 obj.p1 = as_tuple(p1)
                 obj.p2 = p2
+                obj.p3 = p3
                 return obj
 
         grid = Grid(shape=(4, 4))
@@ -1056,7 +1179,7 @@ def __new__(cls, name=None, arguments=None, p0=None, p1=None, p2=None):
         g = Function(name='g', grid=grid)
 
         func = MyDefFunction(name='foo', arguments=f.indexify(),
-                             p0=f, p1=f, p2='bar')
+                             p0=f, p1=f, p2='bar', p3=mapping_type({'field': f}))
 
         mapper = {f: g, f.indexify(): g.indexify()}
         func1 = uxreplace(func, mapper)
@@ -1065,6 +1188,10 @@ def __new__(cls, name=None, arguments=None, p0=None, p1=None, p2=None):
         assert func1.p0 is g
         assert func1.p1 == (g,)
         assert func1.p2 == 'bar'
+        assert type(func1.p3) is mapping_type
+        assert func1.p3['field'] is g
+        assert func.p3['field'] is f
+        assert uxreplace(func1, mapper) is func1
 
     def test_custom_def_function_reconstruction_no_aliasing(self):
 
@@ -1123,6 +1250,44 @@ def test_halo_touch(self):
 
         assert ht1.args == (g[x, y],)
 
+    @pytest.mark.parametrize('op', [UnevalAdd, UnevalMul])
+    def test_lowered_arithmetic_order(self, op):
+        a, b, c = [Temp(name=f'r{i}') for i in (9, 10, 11)]
+        x, y, z = [Temp(name=f'r{i}') for i in (0, 1, 2)]
+        expr = op(c, op(b, a))
+
+        rebuilt = uxreplace(expr, {a: x, b: y, c: z})
+
+        assert rebuilt.args == (z, op(y, x))
+
+    @pytest.mark.parametrize('cls,converted', [(CallFromPointer, CallFromComposite),
+                                               (CallFromComposite, CallFromPointer),
+                                               (FieldFromPointer, FieldFromComposite),
+                                               (FieldFromComposite, FieldFromPointer)])
+    def test_bound_member_names(self, cls, converted):
+        a, b = [dSymbol(name=f'r{i}') for i in (9, 0)]
+        p, q = [dSymbol(name=f'p{i}') for i in (9, 0)]
+        expr = cls(a, p, [a])
+
+        assert uxreplace(expr, {a: b, p: q}) == cls(a, q, [b])
+        assert uxreplace(expr, {expr: cls(b, q, [b])}) == cls(b, q, [b])
+        assert uxreplace(expr, {cls: converted, a: b, p: q}) == converted(a, q, [b])
+
+        bound = cls(a, p)
+        assert uxreplace(bound, {a: b}) is bound
+
+    def test_findexed_strides(self):
+        grid = Grid(shape=(4, 4))
+        x, y = grid.dimensions
+        f = Function(name='f', grid=grid)
+        a, b = [dSymbol(name=f'y_stride{i}') for i in (9, 0)]
+        expr = FIndexed(f.indexed, x, y, strides_map={y: a})
+
+        rebuilt = uxreplace(expr, {a: b})
+
+        assert rebuilt.strides == (b,)
+        assert expr.strides == (a,)
+
 
 def test_minmax():
     grid = Grid(shape=(5, 5))
diff --git a/tests/test_visitors.py b/tests/test_visitors.py
index bdbd5dd6da..e59d8798fb 100644
--- a/tests/test_visitors.py
+++ b/tests/test_visitors.py
@@ -11,7 +11,7 @@
     FindNodes, FindSections, FindSymbols, FindWithin, IsPerfectIteration, Iteration, List,
     MapNodes, Node, Transformer, Uxreplace, printAST
 )
-from devito.symbolics import ListInitializer
+from devito.symbolics import ListInitializer, retrieve_indexed
 from devito.types import Array, LocalObject, SpaceDimension, Symbol
 
 
@@ -422,7 +422,8 @@ def test_find_symbols_with_duplicates():
     # objects because they are constructed in two different places during
     # the CIRE pass
     r0a = exprs[0].output
-    r0b = exprs[1].expr.rhs.args[0].args[0]
+    r0b, = [i for i in retrieve_indexed(exprs[1].expr.rhs)
+            if i.name == r0a.name]
     assert r0a.function is r0b.function
     assert r0a.name == r0b.name == "r0"
     assert hash(r0a) == hash(r0b)