diff --git a/devito/core/cpu.py b/devito/core/cpu.py
index 3435a472e6..8bd9934ad1 100644
--- a/devito/core/cpu.py
+++ b/devito/core/cpu.py
@@ -97,6 +97,7 @@ def _normalize_kwargs(cls, **kwargs):
         o['linearize'] = oo.pop('linearize', cls.LINEARIZE)
         o['mapify-reduce'] = oo.pop('mapify-reduce', cls.MAPIFY_REDUCE)
         o['index-mode'] = oo.pop('index-mode', cls.INDEX_MODE)
+        o['split-modulo-buffers'] = False
         o['place-transfers'] = oo.pop('place-transfers', True)
         o['errctl'] = oo.pop('errctl', cls.ERRCTL)
         o['scalar-min-type'] = oo.pop('scalar-min-type', cls.SCALAR_MIN_TYPE)
diff --git a/devito/core/gpu.py b/devito/core/gpu.py
index 536a06422e..7fe399771f 100644
--- a/devito/core/gpu.py
+++ b/devito/core/gpu.py
@@ -115,6 +115,7 @@ def _normalize_kwargs(cls, **kwargs):
         o['linearize'] = oo.pop('linearize', False)
         o['mapify-reduce'] = oo.pop('mapify-reduce', cls.MAPIFY_REDUCE)
         o['index-mode'] = oo.pop('index-mode', cls.INDEX_MODE)
+        o['split-modulo-buffers'] = False
         o['place-transfers'] = oo.pop('place-transfers', True)
         o['errctl'] = oo.pop('errctl', cls.ERRCTL)
         o['scalar-min-type'] = oo.pop('scalar-min-type', cls.SCALAR_MIN_TYPE)
@@ -151,6 +152,8 @@ def wrapper(expressions, mode='default', options=None, **kwargs1):
             options = options or {}
 
             if mode == 'host':
+                options = {'split-modulo-buffers': options0['split-modulo-buffers'],
+                           **options}
                 target = {
                     'platform': 'cpu64',
                     'language': 'C' if options0['par-disabled'] else 'openmp',
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/nodes.py b/devito/ir/iet/nodes.py
index 7c12b39647..f340b0fea5 100644
--- a/devito/ir/iet/nodes.py
+++ b/devito/ir/iet/nodes.py
@@ -18,7 +18,7 @@
     AFFINE, INBOUND, PARALLEL, PARALLEL_IF_ATOMIC, PARALLEL_IF_PVT, SEQUENTIAL,
     VECTORIZED, Forward, PrefetchUpdate, Property, WithLock
 )
-from devito.symbolics import CallFromPointer, ListInitializer
+from devito.symbolics import CallFromPointer, IndexedPointer, ListInitializer
 from devito.tools import (
     Signer, as_tuple, ctypes_to_cstr, filter_ordered, filter_sorted, flatten
 )
@@ -344,6 +344,10 @@ def _rebuild(self, *args, **kwargs):
             kwargs['arguments'] = [mapper.get(i, i) for i in self.arguments]
         return super()._rebuild(**kwargs)
 
+    def rebuild_args(self, func):
+        """Rebuild the Call, applying `func` to its callee's arguments."""
+        return self._rebuild(arguments=func(self.arguments))
+
     @property
     def children(self):
         return tuple(i for i in self.arguments if isinstance(i, (Call, Lambda)))
@@ -394,7 +398,7 @@ def expr_symbols(self):
         if self.base is not None:
             retval.append(self.base)
 
-        if isinstance(self.retobj, Indexed):
+        if isinstance(self.retobj, (Indexed, IndexedPointer)):
             retval.extend(self.retobj.free_symbols)
         elif self.retobj is not None:
             retval.append(self.retobj)
diff --git a/devito/ir/iet/visitors.py b/devito/ir/iet/visitors.py
index 6186ff4e76..1aec46b257 100644
--- a/devito/ir/iet/visitors.py
+++ b/devito/ir/iet/visitors.py
@@ -23,7 +23,7 @@
 )
 from devito.ir.support.space import Backward
 from devito.symbolics import (
-    FieldFromComposite, FieldFromPointer, ListInitializer, uxreplace
+    FieldFromComposite, FieldFromPointer, IndexedPointer, ListInitializer, uxreplace
 )
 from devito.symbolics.extended_dtypes import NoDeclStruct
 from devito.tools import (
@@ -342,7 +342,7 @@ def _gen_value(self, obj, mode=1, masked=()):
                 strtype = f'{strtype}{self._restrict_keyword}'
         strtype = ' '.join(qualifiers + [strtype])
 
-        if obj.is_LocalObject and obj._C_modifier is not None and mode == 2:
+        if obj._C_modifier is not None and mode == 2:
             strtype += obj._C_modifier
 
         strname = obj._C_name
@@ -615,7 +615,7 @@ def visit_Call(self, o, nested_call=False):
             call = MultilineCall(o.name, arguments, True, o.is_indirect, cast,
                                  o.templates)
             if retobj.is_Indexed or \
-               isinstance(retobj, (FieldFromComposite, FieldFromPointer)):
+               isinstance(retobj, (FieldFromComposite, FieldFromPointer, IndexedPointer)):
                 return c.Assign(self.ccode(retobj), call)
             else:
                 return c.Initializer(c.Value(rettype, retobj._C_name), call)
@@ -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 13927d5761..38b2d4dc28 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
@@ -477,6 +479,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):
@@ -520,15 +527,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..3ed0e1efa7 100644
--- a/devito/passes/iet/engine.py
+++ b/devito/passes/iet/engine.py
@@ -14,15 +14,17 @@
 from devito.mpi.distributed import MPINeighborhood
 from devito.mpi.routines import Gather, HaloUpdate, HaloWait, MPIMsg, Scatter
 from devito.passes import needs_transfer
-from devito.symbolics import FieldFromComposite, FieldFromPointer, IndexedPointer, search
+from devito.symbolics import (
+    FieldFromComposite, FieldFromPointer, IndexedPointer, search, uxreplace
+)
 from devito.tools import (
     DAG, as_hashable, as_tuple, filter_ordered, memoized_func, natural_sort_key,
     sorted_priority, timed_pass
 )
 from devito.types import (
     Array, Auto, Bundle, ComponentAccess, CompositeObject, FunctionMap, IncrDimension,
-    Indirection, ModuloDimension, NPThreads, NThreadsBase, Pointer, SharedData, Symbol,
-    Temp, ThreadArray, Wildcard
+    IndexedData, Indirection, ModuloDimension, NPThreads, NThreadsBase, Pointer,
+    SharedData, Symbol, Temp, ThreadArray, Wildcard
 )
 from devito.types.args import ArgProvider
 from devito.types.dense import DiscreteFunction
@@ -134,6 +136,10 @@ def apply(self, func, **kwargs):
 
         Callable parameters and Call arguments are reconciled before the graph
         walk, after each changed node, and after the pass has completed.
+
+        The returned metadata may include `args`, a mapping from new parameters
+        to expressions in the original parameters. Each call site instantiates
+        these expressions using its own argument bindings.
         """
         _update_args(self)
 
@@ -166,7 +172,7 @@ def apply(self, func, **kwargs):
             efuncs[i] = efunc
             efuncs.update(dict([(i.name, i) for i in new_efuncs]))
 
-            efuncs = _update_args_efunc(efunc, efuncs, dag)
+            efuncs = _update_args_efunc(efunc, efuncs, dag, metadata.get('args'))
 
         # Minimize code size
         if len(efuncs) > len(self.efuncs):
@@ -564,6 +570,8 @@ def abstract_objects(objects0):
     mapper = {}
     sregistry = SymbolRegistry()
     for i in objects:
+        if i in mapper:
+            continue
         abstract_object(i, mapper, sregistry)
 
     return mapper
@@ -723,7 +731,7 @@ def _(i, mapper, sregistry):
     mapper[i] = i._rebuild(name=sregistry.make_name(prefix='nthreads'))
 
 
-def _update_args_efunc(root, efuncs, dag):
+def _update_args_efunc(root, efuncs, dag, args=None):
     """
     Re-derive the parameters of `root` and apply the changes in cascade through
     the `efuncs`.
@@ -775,7 +783,8 @@ def _update_args_efunc(root, efuncs, dag):
     # 4) removed a function that was previously necessary
     functions = FindSymbols('symbolics').visit(root.body)
     drop_params.extend(a for a in root.parameters
-                       if a.is_AbstractFunction and a not in functions)
+                       if (a.is_AbstractFunction or isinstance(a, IndexedData))
+                       and a.function not in functions)
 
     # Must record the index, not the param itself, since a param may be
     # bound to whatever arg, possibly a generic SymPy expr
@@ -808,9 +817,16 @@ def _filter(v, efunc=None):
     efuncs = dict(efuncs)
     efuncs[root.name] = root._rebuild(parameters=_filter(root.parameters, root))
 
+    def _filter_call(values):
+        mapper = {}
+        if args:
+            bindings = dict(zip(root.parameters, values, strict=True))
+            mapper = {p: uxreplace(expr, bindings) for p, expr in args.items()}
+        return [mapper.get(a, a) for a in _filter(values)]
+
     # Update all call sites to use the new signature
     for n in dag.downstream(root.name):
-        mapper = {c: c._rebuild(arguments=_filter(c.arguments))
+        mapper = {c: c.rebuild_args(_filter_call)
                   for c in FindNodes(Call).visit(efuncs[n])
                   if c.name == root.name}
         if not mapper:
diff --git a/devito/passes/iet/misc.py b/devito/passes/iet/misc.py
index 4facf6646d..ddd6fa814c 100644
--- a/devito/passes/iet/misc.py
+++ b/devito/passes/iet/misc.py
@@ -8,8 +8,7 @@
 from devito.finite_differences.differentiable import SafeInv
 from devito.ir import (
     Any, DummyExpr, EmptyList, FindApplications, FindNodes, FindSymbols, Forward,
-    Iteration, Prodder, Transformer, Uxreplace, filter_iterations, pull_dims,
-    retrieve_iteration_tree
+    Iteration, Prodder, Transformer, Uxreplace, filter_iterations, retrieve_iteration_tree
 )
 from devito.ir.iet.efunc import DeviceFunction, EntryFunction
 from devito.passes.iet.engine import iet_pass
@@ -188,8 +187,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 +199,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)
 
@@ -371,10 +371,8 @@ def _(iet):
 
 
 def _rename_subdims(target, dimensions):
-    # Find SubDimensions or SubDimension-derived dimensions used as indices in
-    # the expression
-    indexeds = FindSymbols('indexeds').visit(target)
-    dims = pull_dims(indexeds, flag=False)
+    # Find SubDimensions or SubDimension-derived dimensions used in expressions
+    dims = FindSymbols('dimensions').visit(target)
     dims = [d for d in dims if any(dim.is_AbstractSub for dim in d._defines)]
     dims = [d for d in dims if not d.is_SubIterator]
     names = [d.root.name for d in dims]
diff --git a/devito/symbolics/extended_sympy.py b/devito/symbolics/extended_sympy.py
index eeb290feee..4f8df6d914 100644
--- a/devito/symbolics/extended_sympy.py
+++ b/devito/symbolics/extended_sympy.py
@@ -547,6 +547,8 @@ def base(self):
     def index(self):
         return self._index
 
+    indices = index
+
     def __str__(self):
         indices = ''.join(f'[{i}]' for i in self.index)
         return f"{self.base}{indices}"
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/devito/types/array.py b/devito/types/array.py
index 1988e2b0a8..98f8932813 100644
--- a/devito/types/array.py
+++ b/devito/types/array.py
@@ -645,9 +645,11 @@ class ComponentAccess(Expr, Pickable):
     __rkwargs__ = ('index',)
 
     def __new__(cls, arg, index=0, **kwargs):
-        if not (arg.is_Symbol or arg.is_Indexed):
-            raise ValueError(f"Expected Symbol or Indexed, got `{type(arg)}` "
-                             "instead")
+        from devito.symbolics import IndexedPointer
+
+        if not (arg.is_Symbol or arg.is_Indexed or isinstance(arg, IndexedPointer)):
+            raise ValueError(f"Expected Symbol, Indexed or IndexedPointer, "
+                             f"got `{type(arg)}` instead")
         if not is_integer(index) or index > 3:
             raise ValueError("Expected 0 <= index < 4")
 
diff --git a/devito/types/basic.py b/devito/types/basic.py
index fc8940f5e9..bb3b31a922 100644
--- a/devito/types/basic.py
+++ b/devito/types/basic.py
@@ -58,6 +58,12 @@ class CodeSymbol:
     for a given architecture.
     """
 
+    _C_modifier = None
+    """
+    A suffix added to the type in a function signature. E.g., `&` for a
+    reference parameter, or ` const` for a constant pointer.
+    """
+
     @abc.abstractmethod
     def __init__(self, *args, **kwargs):
         return
@@ -2057,13 +2063,6 @@ def _mem_internal_eager(self):
     def _mem_internal_lazy(self):
         return self._liveness == 'lazy'
 
-    """
-    A modifier added to the declaration of the LocalType when it appears in a
-    function signature. For example, a subclass might define `_C_modifier = '&'`
-    to impose pass-by-reference semantics.
-    """
-    _C_modifier = None
-
     """
     One or more optional keywords added to the declaration of the LocalType
     in between the type and the variable name when it appears in a function
diff --git a/devito/types/dense.py b/devito/types/dense.py
index 30f543a47c..96c99eb33b 100644
--- a/devito/types/dense.py
+++ b/devito/types/dense.py
@@ -23,7 +23,7 @@
 from devito.symbolics import FieldFromPointer, IndexedPointer, normalize_args
 from devito.tools import (
     ReducerMap, as_tuple, c_restrict_void_p, dtype_to_ctype, flatten, humanbytes,
-    is_integer, memoized_meth, mpi4py_mapper
+    is_integer, memoized_meth, mpi4py_mapper, prod
 )
 from devito.types.args import ArgProvider
 from devito.types.basic import AbstractFunction
@@ -1351,11 +1351,17 @@ def _arg_check(self, args, intervals, **kwargs):
         if args.options['index-mode'] == 'int32' and \
            args.options['linearize'] and \
            data.size - 1 >= np.iinfo(np.int32).max:
-            raise InvalidArgument(f"`{self.name}`, with its {data.size} elements, is too "
-                                  "big for int32 pointer arithmetic. Consider using the "
-                                  "'index-mode=int64' option, the save=Buffer(..) "
-                                  "API (TimeFunction only), or domain "
-                                  "decomposition via MPI")
+            size = data.size
+            if args.options['split-modulo-buffers'] and len(data.shape) == self.ndim:
+                size = prod(s for d, s in zip(self.dimensions, data.shape, strict=True)
+                            if not d.is_Stepping)
+            if size - 1 >= np.iinfo(np.int32).max:
+                raise InvalidArgument(
+                    f"`{self.name}`, with its {size} indexed elements, is too "
+                    "big for int32 pointer arithmetic. Consider using the "
+                    "'index-mode=int64' option, the save=Buffer(..) "
+                    "API (TimeFunction only), or domain decomposition via MPI"
+                )
 
         return super()._arg_check(args, intervals, **kwargs)
 
diff --git a/devito/types/misc.py b/devito/types/misc.py
index 61a2860148..e94c1fe3d8 100644
--- a/devito/types/misc.py
+++ b/devito/types/misc.py
@@ -110,6 +110,11 @@ def __new__(cls, base, *args, strides_map=None, accessor=None):
 
         return obj
 
+    @classmethod
+    def class_key(cls):
+        # Keep arithmetic ordering consistent with unlinearized accesses
+        return Indexed.class_key()
+
     def __repr__(self):
         return "{}({})".format(self.name, ", ".join(str(i) for i in self.indices))
 
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..b9452d8025 100644
--- a/tests/test_iet.py
+++ b/tests/test_iet.py
@@ -6,21 +6,22 @@
 import sympy
 
 from devito import (  # noqa
-    Dimension, Eq, Function, Grid, Operator, TimeFunction, switchconfig
+    Dimension, Eq, Function, Grid, Operator, SubDimension, TimeFunction, switchconfig
 )
 from devito.ir import SymbolRegistry
 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,
-    retrieve_iteration_tree
+    Call, Callable, CGen, Conditional, Definition, Dereference, DeviceCall,
+    DeviceFunction, DummyExpr, ElementalFunction, EntryFunction, FindNodes, FindSymbols,
+    Iteration, KernelLaunch, Lambda, List, Switch, Transformer, Uxreplace,
+    filter_iterations, make_callable, make_efunc, retrieve_iteration_tree
 )
 from devito.ir.iet.visitors import sorted_efuncs
 from devito.passes.iet.engine import Graph, reuse_efuncs
 from devito.passes.iet.languages.C import CDataManager
+from devito.passes.iet.misc import abridge_dim_names
 from devito.symbolics import (
-    FLOAT, Byref, Class, FieldFromComposite, InlineIf, ListInitializer, Macro, SizeOf,
-    String
+    FLOAT, Byref, Class, FieldFromComposite, IndexedPointer, InlineIf, ListInitializer,
+    Macro, SizeOf, String
 )
 from devito.symbolics.extended_dtypes import c_complex
 from devito.tools import CustomDtype, as_tuple, dtype_to_ctype
@@ -120,6 +121,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"]')
@@ -430,17 +445,20 @@ def test_call_indexed():
 }"""
 
 
-def test_call_retobj_indexed():
+@pytest.mark.parametrize('pointer', [False, True])
+def test_call_retobj_indexed(pointer):
     grid = Grid(shape=(10, 10))
 
     u = Function(name='u', grid=grid)
     v = Function(name='v', grid=grid)
 
-    call = Call('foo', [u], retobj=v.indexify())
+    retobj = IndexedPointer(v.indexed, grid.dimensions) if pointer else v.indexify()
+    call = Call('foo', [u], retobj=retobj)
 
     assert str(call) == "v[x][y] = foo(u_vec);"
 
     assert not call.defines
+    assert set(grid.dimensions).issubset(call.expr_symbols)
 
 
 def test_call_lambda_transform():
@@ -542,14 +560,69 @@ def test_codegen_quality0():
     assert foo1.parameters[0] is a
 
 
-def test_reuse_efuncs_natural_numbering():
+def test_const_pointer_parameter():
+    class ConstPointer(Symbol):
+        _C_modifier = ' const'
+
+    p = ConstPointer(name='p', dtype=POINTER(dtype_to_ctype(np.float32)))
+    foo = Callable('foo', Call('external', (p,)), 'void', parameters=(p,))
+
+    assert 'void foo(float * const p)' in str(foo)
+
+
+def test_graph_parameter_bindings():
+    x, y, a, b, p = [Symbol(name=n) for n in ('x', 'y', 'a', 'b', 'p')]
+    inner = Callable('inner', Call('external', (x, y)), 'void', parameters=(x, y))
+    root = make_callable('root', [Call('inner', (a + 1, b)),
+                                  Call('inner', (b, a + 2))])
+    graph = Graph(root)
+    graph.apply(lambda iet: (iet, {'efuncs': [inner]}))
+
+    def transform(iet):
+        if iet.name != 'inner':
+            return iet, {}
+        return iet._rebuild(body=Call('external', (p,))), {'args': {p: x + y}}
+
+    graph.apply(transform)
+
+    assert graph.efuncs['inner'].parameters == (p,)
+    assert [c.arguments for c in FindNodes(Call).visit(graph.root)] == [
+        (a + b + 1,), (a + b + 2,)
+    ]
+
+
+@pytest.mark.parametrize('ftype', [EntryFunction, DeviceFunction])
+@pytest.mark.parametrize('pointer', [False, True])
+@pytest.mark.parametrize('clash', [False, True])
+def test_abridge_pointer_dimensions(ftype, pointer, clash):
+    x = Dimension(name='x')
+    xi = SubDimension.left(name='xi', parent=x, thickness=4)
+    a = Array(name='a', dimensions=(x,), shape=(16,))
+    index = xi + x if clash else xi
+    access = IndexedPointer(a, index) if pointer else a[index]
+    loop = Iteration(DummyExpr(access, access + 1), xi,
+                     (xi.symbolic_min, xi.symbolic_max, 1))
+    body = Iteration(loop, x, (0, 15, 1)) if clash else loop
+    foo = ftype('foo', body, 'void', parameters=(a,))
+
+    foo = abridge_dim_names(foo)
+
+    inner = FindNodes(Iteration).visit(foo)[-1]
+    assert inner.dim.name == ('xi' if clash else 'x')
+    assert inner.limits == loop.limits
+
+
+@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_linearize.py b/tests/test_linearize.py
index c78880924d..cfd019f9c8 100644
--- a/tests/test_linearize.py
+++ b/tests/test_linearize.py
@@ -6,6 +6,8 @@
     Eq, Function, Grid, Inc, MatrixSparseTimeFunction, Operator, SparseTimeFunction,
     TimeFunction, configuration, sin, switchconfig
 )
+from devito.data.allocators import DataReference
+from devito.exceptions import InvalidArgument
 from devito.ir import Call, Callable, DummyExpr, Expression, FindNodes, SymbolRegistry
 from devito.passes import Graph, generate_macros, linearize
 from devito.types import Array, Bundle, CustomDimension, DefaultDimension
@@ -161,6 +163,24 @@ def test_interpolation_enforcing_int64_indexing_v2():
              opt=('advanced', {'linearize': True, 'index-mode': 'int32'}))
 
 
+@switchconfig(autopadding=False)
+@pytest.mark.parametrize('index_mode', ['int32', 'int64'])
+def test_large_allocation(index_mode):
+    grid = Grid(shape=(1024, 1024, 1024))
+    # Validate the allocation size without allocating or accessing its storage
+    data = np.broadcast_to(np.zeros(1, dtype=np.float32), (3, *grid.shape))
+    u = TimeFunction(name='u', grid=grid, time_order=2, space_order=0,
+                     allocator=DataReference(data))
+    op = Operator(Eq(u.forward, u + 1), name='large',
+                  opt=('advanced', {'linearize': True, 'index-mode': index_mode}))
+
+    if index_mode == 'int32':
+        with pytest.raises(InvalidArgument, match='int32 pointer arithmetic'):
+            op.arguments(time_m=0, time_M=2)
+    else:
+        op.arguments(time_m=0, time_M=2)
+
+
 def test_interpolation_msf():
     grid = Grid(shape=(4, 4))
 
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..067b05b35f 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, IndexedPointer, 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)"
@@ -627,6 +747,32 @@ def test_findexed():
     assert new_fi.strides_map == strides_map
 
 
+@pytest.mark.parametrize('name', ['a', 'z'])
+@pytest.mark.parametrize('kind', ['sum', 'product'])
+@pytest.mark.parametrize('form', ['indexed', 'linearized', 'pointer'])
+def test_access_arithmetic_order(name, kind, form):
+    grid = Grid(shape=(8,))
+    x, = grid.dimensions
+    f = Function(name=name, grid=grid)
+    s = Array(name='s', dimensions=(x,), shape=grid.shape)
+    r = dSymbol('r', dtype=np.float32)
+    accesses = {
+        'indexed': (f[x], f'{name}[x]'),
+        'linearized': (FIndexed(f.indexed, x), f'{name}(x)'),
+        'pointer': (IndexedPointer(dSymbol(f'{name}_buf0'), x), f'{name}_buf0[x]')
+    }
+    access, code = accesses[form]
+    if kind == 'sum':
+        expr = r*access - 2*r*s[x]
+        expected = f'-2*r*s[x] + r*{code}'
+    else:
+        expr = access*s[x]
+        expected = f's[x]*{code}'
+
+    assert ccode(expr) == expected
+    assert ccode(unevaluate(expr)) == expected
+
+
 def test_component_access():
     grid = Grid(shape=(3, 3, 3))
     x, y, z = grid.dimensions
@@ -645,11 +791,16 @@ def test_component_access():
     assert cf2 == cf1
 
 
-def test_component_access_symbol_printing():
+@pytest.mark.parametrize('pointer,expected', [
+    (False, '1.250F*acc.x'), (True, '1.250F*acc[0].x')
+])
+def test_component_access_symbol_printing(pointer, expected):
     acc = dSymbol(name='acc', dtype=dtypes_vector_mapper[(np.float32, 4)])
+    if pointer:
+        acc = IndexedPointer(acc, 0)
     expr = ComponentAccess(acc, 0)
 
-    assert ccode(1.25*expr, dtype=expr.dtype) == '1.250F*acc.x'
+    assert ccode(1.25*expr, dtype=expr.dtype) == expected
 
 
 def test_vector_access():
@@ -1037,17 +1188,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 +1210,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 +1219,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 +1281,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)