diff --git a/finat/physically_mapped.py b/finat/physically_mapped.py
index 511d3e146..cf202c91e 100644
--- a/finat/physically_mapped.py
+++ b/finat/physically_mapped.py
@@ -1,5 +1,6 @@
 from abc import ABCMeta, abstractmethod
-from collections.abc import Mapping
+from collections.abc import Iterable, Mapping
+from numbers import Number
 
 import gem
 import numpy
@@ -7,6 +8,33 @@
 from finat.citations import cite
 
 
+zero = gem.Zero()
+one = gem.Literal(1.0)
+
+
+def _as_basis_entry(value: object) -> gem.Node:
+    """Normalize a basis-transformation matrix entry.
+
+    Parameters
+    ----------
+    value
+        Scalar numerical or GEM matrix entry.
+
+    Returns
+    -------
+    gem.Node
+        Scalar GEM entry with numerical zero and one represented
+        symbolically.
+
+    """
+    if isinstance(value, Number):
+        if value == 0:
+            return zero
+        if value == 1:
+            return one
+    return gem.as_gem(value)
+
+
 class NeedsCoordinateMappingElement(metaclass=ABCMeta):
     """Abstract class for elements that require physical information
     either to map or construct their basis functions."""
@@ -16,38 +44,81 @@ def dual_transformation(self, Q, coordinate_mapping=None):
 
 
 class MappedTabulation(Mapping):
-    """A lazy tabulation dict that applies the basis transformation only
-    on the requested derivatives.
+    """Apply a sparse basis transformation to reference tabulations.
+
+    Parameters
+    ----------
+    M : gem.ListTensor
+        Basis-transformation matrix.
+    ref_tabulation : Mapping
+        Reference tabulations indexed by derivative order.
+    indices : iterable of int, optional
+        Rows retained by an element restriction.
+
+    Notes
+    -----
+    The transformation is stored in compressed sparse row form.  One
+    constant table holds the column of each nonzero, and a second holds its
+    value.  Rows shorter than the longest are padded with zero values, which
+    contribute nothing to the sum.
+
+    Every row therefore contracts over the same number of entries.  The
+    basis axis stays one loop, so the transformation reaches the quadrature
+    contraction as a linear map over the reference tabulation.
 
-    :arg M: a gem.ListTensor with the basis transformation matrix.
-    :arg ref_tabulation: a dict of tabulations on the reference cell.
-    :kwarg indices: an optional list of restriction indices on the basis functions.
     """
-    def __init__(self, M, ref_tabulation, indices=None):
-        self.M = M
+
+    def __init__(
+            self, M: gem.ListTensor, ref_tabulation: Mapping,
+            indices: Iterable[int] | None = None) -> None:
         self.ref_tabulation = ref_tabulation
         if indices is None:
-            indices = list(range(M.shape[0]))
-        self.indices = indices
-        # we expect M to be sparse with O(1) nonzeros per row
-        # for each row, get the column index of each nonzero entry
-        csr = [[j for j in range(M.shape[1]) if not isinstance(M.array[i, j], gem.Zero)]
-               for i in indices]
-        self.csr = csr
+            indices = range(M.shape[0])
+        self.indices = tuple(indices)
+
+        nonzero_rows = []
+        for source_row in self.indices:
+            row = []
+            for column in range(M.shape[1]):
+                value = _as_basis_entry(M.array[source_row, column])
+                if not isinstance(value, gem.Zero):
+                    row.append((column, value))
+            nonzero_rows.append(row)
+        width = max((len(row) for row in nonzero_rows), default=0)
+        nrows = len(self.indices)
+        columns = numpy.zeros((nrows, width), dtype=gem.uint_type)
+        data = numpy.full((nrows, width), zero, dtype=object)
+        for index, row in enumerate(nonzero_rows):
+            columns[index, :len(row)] = tuple(column for column, _ in row)
+            data[index, :len(row)] = tuple(
+                gem.as_gem(value) for _, value in row)
+        self._width = width
+        self._columns = gem.Literal(columns, dtype=gem.uint_type)
+        self._values = gem.ListTensor(data)
         self._tabulation_cache = {}
 
-    def matvec(self, table):
-        # basis recombination using hand-rolled sparse-dense matrix multiplication
-        ii = gem.indices(len(table.shape)-1)
-        phi = [gem.Indexed(table, (j, *ii)) for j in range(self.M.shape[1])]
-        # the sum approach is faster than calling numpy.dot or gem.IndexSum
-        exprs = [gem.ComponentTensor(gem.Sum(*(self.M.array[i, j] * phi[j] for j in js)), ii)
-                 for i, js in zip(self.indices, self.csr)]
+    def matvec(self, table: gem.Node) -> gem.Node:
+        """Transform one reference tabulation.
+
+        Parameters
+        ----------
+        table
+            Reference tabulation with the basis axis first.
 
-        result = gem.ListTensor(exprs)
-        result, = gem.optimise.unroll_indexsum((result,), lambda index: True)
-        # result = gem.optimise.aggressive_unroll(self.M @ table)
-        return result
+        Returns
+        -------
+        gem.Node
+            Tabulation whose first axis is the transformed basis axis.
+
+        """
+        tail = gem.indices(len(table.shape) - 1)
+        row = gem.Index(extent=len(self.indices))
+        entry = gem.Index(extent=self._width)
+        column = gem.VariableIndex(gem.Indexed(self._columns, (row, entry)))
+        basis = gem.Product(gem.Indexed(self._values, (row, entry)),
+                            gem.Indexed(table, (column, *tail)))
+        mapped = gem.IndexSum(basis, (entry,))
+        return gem.ComponentTensor(mapped, (row, *tail))
 
     def __getitem__(self, alpha):
         try:
@@ -195,10 +266,6 @@ def physical_vertices(self):
                 (gdim, )."""
 
 
-zero = gem.Zero()
-one = gem.Literal(1.0)
-
-
 def identity(*shape):
     V = numpy.eye(*shape, dtype=object)
     for multiindex in numpy.ndindex(V.shape):
diff --git a/gem/coffee.py b/gem/coffee.py
index e9ae20119..49b1391c5 100644
--- a/gem/coffee.py
+++ b/gem/coffee.py
@@ -1,17 +1,28 @@
-"""This module contains an implementation of the COFFEE optimisation
-algorithm operating on a GEM representation.
+"""Eliminate sharing in multilinear GEM expressions.
 
-This file is NOT for code generation as a COFFEE AST.
+The input :class:`MonomialSum` is the normal form of a finite element
+integrand: its atomics are the linear operands and its ``rest`` is
+independent of the multilinear loops.  COFFEE chooses scalar factorizations
+between those operands.  When sum factorization supplies a contraction
+order, the same sharing elimination is applied at every reduction level.
+This realizes sum factorization as generalized code motion while retaining
+the finite element maps selected during argument factorization.
+
+This module transforms GEM; it does not generate a COFFEE AST.
 """
 
-from itertools import chain, repeat
+from collections import OrderedDict, defaultdict
+from collections.abc import Sequence
+from itertools import chain
 import logging
 
 import numpy
 
-from gem.gem import IndexSum, one
-from gem.optimise import make_sum, make_product
-from gem.refactorise import Monomial
+from gem.gem import ComponentTensor, Index, Indexed, IndexSum, Node, one
+from gem.node import MemoizerArg
+from gem.contraction import has_arithmetic, partition_connected
+from gem.optimise import filtered_replace_indices, make_sum, make_product
+from gem.refactorise import Monomial, MonomialSum
 from gem.utils import groupby
 
 
@@ -47,95 +58,129 @@ def index_extent(factor, linear_indices):
     return numpy.prod([i.extent for i in factor.free_indices if i in linear_indices])
 
 
-def sort_monomials(monomials):
-    """Sort monomials to produce a better initial guess for :func:`find_optimal_atomics`.
-
-    :arg monomials: A list of :class:`Monomial`s
+def find_optimal_atomics(
+        monomials: Sequence[Monomial],
+        linear_indices: tuple) -> tuple[Node, ...]:
+    """Find a minimum-cost set of atomics intersecting every monomial.
 
-    :returns: the reordered list of monomials.
-    """
-    if len(monomials) <= 2:
-        return monomials
-    # Construct a monomial subset with non-intersecting atomics
-    head = []
-    rest = []
-    atomics = set()
-    for m in monomials:
-        if atomics.intersection(m.atomics):
-            rest.append(m)
-        else:
-            atomics.update(m.atomics)
-            head.append(m)
-    # Put non-intersecting subset first and recurse on the rest
-    monomials = head + sort_monomials(rest)
-    return monomials
+    Parameters
+    ----------
+    monomials
+        Monomials with equal contraction indices.
+    linear_indices
+        Free indices belonging to form arguments.
 
+    Returns
+    -------
+    tuple of Node
+        Atomics selected for common-subexpression factorization.
 
-def find_optimal_atomics(monomials, linear_indices):
-    """Find optimal atomic common subexpressions, which produce least number of
-    terms in the resultant IndexSum when factorised.
-
-    :arg monomials: A list of :class:`Monomial`s, all of which should have
-                    the same sum indices
-    :arg linear_indices: tuple of linear indices
-
-    :returns: list of atomic GEM expressions
     """
-    monomials = sort_monomials(monomials)
-
-    atomics = tuple(dict.fromkeys(chain.from_iterable(monomial.atomics for monomial in monomials)))
-
-    # Create a list of sets of indices to avoid any hashing during the search
-    monomial_atomics = [set(map(atomics.index, m.atomics)) for m in monomials]
-
-    # Precompute the cost of each atomic
-    atomic_costs = list(map(index_extent, atomics, repeat(linear_indices)))
-
-    def cost(solution):
-        extent = sum(atomic_costs[i] for i in solution)
-        # Prefer shorter solutions, but larger extents
-        return (len(solution), -extent)
-
-    optimal_solution = set(range(len(atomics)))  # pessimal but feasible solution
-    optimal_cost = cost(optimal_solution)
-    solution = set()
-    solution_cost = (0, 0)
-
-    max_it = 1 << 12
-    it = iter(range(max_it))
-
-    def solve(idx):
-        nonlocal solution_cost, optimal_cost
-
-        while idx < len(monomials) and solution.intersection(monomial_atomics[idx]):
-            idx += 1
-
-        if idx < len(monomials):
-            if len(solution) < len(optimal_solution):
-                for atomic in monomial_atomics[idx]:
-                    atomic_cost = atomic_costs[atomic]
-                    old_solution_cost = solution_cost
-                    solution_cost = (solution_cost[0]+1, solution_cost[1]-atomic_cost)
-                    if solution_cost < optimal_cost:
-                        solution.add(atomic)
-                        solve(idx + 1)
-                        solution.remove(atomic)
-                    solution_cost = old_solution_cost
-        else:
-            if solution_cost < optimal_cost:
-                optimal_solution.clear()
-                optimal_solution.update(solution)
-                optimal_cost = solution_cost
-            next(it)
+    atomics = tuple(dict.fromkeys(chain.from_iterable(
+        monomial.atomics for monomial in monomials)))
+    if not atomics:
+        return ()
+
+    positions = {atomic: position
+                 for position, atomic in enumerate(atomics)}
+    constraints = tuple(dict.fromkeys(
+        sum({1 << positions[atomic] for atomic in monomial.atomics})
+        for monomial in monomials))
+    constraints = tuple(
+        constraint for constraint in constraints
+        if not any(
+            other != constraint and other & constraint == other
+            for other in constraints))
+    costs = tuple(int(index_extent(atomic, linear_indices))
+                  for atomic in atomics)
+    covers = tuple(sum(
+        1 << position
+        for position, constraint in enumerate(constraints)
+        if constraint & (1 << atomic))
+        for atomic in range(len(atomics)))
+    full_cover = (1 << len(constraints)) - 1
+
+    covered = 0
+    best_solution = 0
+    while covered != full_cover:
+        atomic = max(
+            range(len(atomics)),
+            key=lambda candidate: (
+                ((covers[candidate] & ~covered).bit_count(),
+                 costs[candidate], -candidate)))
+        best_solution |= 1 << atomic
+        covered |= covers[atomic]
+
+    def solution_cost(solution):
+        selected = [position for position in range(len(atomics))
+                    if solution & (1 << position)]
+        return len(selected), -sum(costs[position] for position in selected)
+
+    best_cost = solution_cost(best_solution)
+    seen = {}
+    states = 0
+    max_states = 1 << 16
+
+    def lower_bound(uncovered):
+        disjoint = 0
+        count = 0
+        for position in sorted(
+                (position for position in range(len(constraints))
+                 if uncovered & (1 << position)),
+                key=lambda position: constraints[position].bit_count()):
+            constraint = constraints[position]
+            if not constraint & disjoint:
+                disjoint |= constraint
+                count += 1
+        return count
+
+    def solve(covered, solution, cardinality, extent):
+        nonlocal best_solution, best_cost, states
+        states += 1
+        if states > max_states:
+            raise StopIteration
+
+        cost = cardinality, -extent
+        if seen.get(covered, (numpy.inf, numpy.inf)) <= cost:
+            return
+        seen[covered] = cost
+        if covered == full_cover:
+            best_solution, best_cost = solution, cost
+            return
+
+        uncovered = full_cover & ~covered
+        if cardinality + lower_bound(uncovered) > best_cost[0]:
+            return
+        constraint = min(
+            (position for position in range(len(constraints))
+             if uncovered & (1 << position)),
+            key=lambda position: constraints[position].bit_count())
+        choices = [atomic for atomic in range(len(atomics))
+                   if constraints[constraint] & (1 << atomic)]
+        choices.sort(
+            key=lambda atomic: (
+                (covers[atomic] & uncovered).bit_count(),
+                costs[atomic], -atomic),
+            reverse=True)
+        for atomic in choices:
+            candidate_cost = cardinality + 1, -(extent + costs[atomic])
+            if candidate_cost < best_cost:
+                solve(
+                    covered | covers[atomic],
+                    solution | (1 << atomic),
+                    cardinality + 1,
+                    extent + costs[atomic])
 
     try:
-        solve(0)
+        solve(0, 0, 0, 0)
     except StopIteration:
-        logger = logging.getLogger('tsfc')
-        logger.warning("Solution to ILP problem may not be optimal: search "
-                       "interrupted after examining %d solutions.", max_it)
+        logging.getLogger("tsfc").warning(
+            "Solution to hitting-set problem may not be optimal: search "
+            "interrupted after examining %d states.", max_states)
 
-    return tuple(atomics[i] for i in optimal_solution)
+    return tuple(
+        atomic for position, atomic in enumerate(atomics)
+        if best_solution & (1 << position))
 
 
 def factorise_atomics(monomials, optimal_atomics, linear_indices):
@@ -158,14 +203,13 @@ def group_key(monomial):
             if oa in monomial.atomics:
                 return oa
         assert False, "Expect at least one optimal atomic per monomial."
-    factor_group = groupby(monomials, key=group_key)
-
-    # We should not drop monomials
-    assert sum(len(ms) for _, ms in factor_group) == len(monomials)
+    factor_groups = OrderedDict()
+    for monomial in monomials:
+        factor_groups.setdefault(group_key(monomial), []).append(monomial)
 
     sum_indices = next(iter(monomials)).sum_indices
     new_monomials = []
-    for oa, monomials in factor_group:
+    for oa, monomials in factor_groups.items():
         # Create new MonomialSum for the factorised out terms
         sub_monomials = []
         for monomial in monomials:
@@ -197,20 +241,204 @@ def group_key(monomial):
     return new_monomials
 
 
-def optimise_monomial_sum(monomial_sum, linear_indices):
-    """Choose optimal common atomic subexpressions and factorise a
-    :class:`MonomialSum` object to create a GEM expression.
+def collect_common_rests(monomials):
+    """Group monomials with a common scalar coefficient.
 
-    :arg monomial_sum: a :class:`MonomialSum` object
-    :arg linear_indices: tuple of linear indices
+    This applies ``r*a + r*b = r*(a + b)`` after argument-factor
+    extraction. It exposes the scalar part of a finite element contraction
+    to COFFEE factorisation.
+
+    Parameters
+    ----------
+    monomials : iterable of Monomial
+        Monomials with equal contraction indices.
+
+    Returns
+    -------
+    list of Monomial
+        Monomials after common coefficients have been collected.
 
-    :returns: factorised GEM expression
     """
+    def group_key(monomial):
+        linear = frozenset(chain.from_iterable(
+            atomic.free_indices for atomic in monomial.atomics))
+        return frozenset(monomial.sum_indices), linear, monomial.rest
+
+    result = []
+    for (_, _, rest), group in groupby(monomials, key=group_key):
+        if len(group) > 1 and rest != one and all(
+                monomial.atomics for monomial in group):
+            sum_indices = group[0].sum_indices
+            node = make_sum(tuple(
+                make_product(monomial.atomics) for monomial in group))
+            result.append(Monomial(sum_indices, (node,), rest))
+        else:
+            result.extend(group)
+    return result
+
+
+def _share_linear_maps(
+        monomial_sum: MonomialSum,
+        linear_indices: tuple[Index, ...]) -> MonomialSum:
+    """Share isomorphic maps of distinct multilinear axes.
+
+    Parameters
+    ----------
+    monomial_sum
+        Sum-of-products representation of a multilinear expression.
+    linear_indices
+        Free indices identifying argument axes.
+
+    Returns
+    -------
+    MonomialSum
+        Representation whose repeated linear maps access one tensor.
+
+    Notes
+    -----
+    Test and trial axes use distinct indices even when they apply the same
+    finite element map.  Renaming each axis to a canonical index exposes
+    that isomorphism without inspecting the element family.  Materializing
+    the canonical map is generalized code motion: the basis transformation
+    is evaluated once and both axes index its result.
+
+    """
+    linear_indices = tuple(linear_indices)
+    linear_set = frozenset(linear_indices)
+    canonical = {
+        index.extent: Index(extent=index.extent)
+        for index in linear_indices
+    }
+    replacer = MemoizerArg(filtered_replace_indices)
+    groups = defaultdict(list)
+    for monomial in monomial_sum:
+        for atomic in monomial.atomics:
+            involved = linear_set.intersection(atomic.free_indices)
+            if len(involved) != 1:
+                continue
+            index, = involved
+            normal = replacer(
+                atomic, ((index, canonical[index.extent]),))
+            groups[normal].append((atomic, index))
+
+    replacements = {}
+    for normal, occurrences in groups.items():
+        indices = {index for _, index in occurrences}
+        if len(indices) < 2 or not has_arithmetic((normal,)):
+            continue
+        index = canonical[next(iter(indices)).extent]
+        tensor = ComponentTensor(normal, (index,))
+        replacements.update(
+            (atomic, Indexed(tensor, (original,)))
+            for atomic, original in occurrences)
+
+    if not replacements:
+        return monomial_sum
+
+    result = MonomialSum()
+    for monomial in monomial_sum:
+        result.add(
+            monomial.sum_indices,
+            tuple(replacements.get(atomic, atomic)
+                  for atomic in monomial.atomics),
+            monomial.rest,
+        )
+    return result
+
+
+def optimise_monomial_sum(
+        monomial_sum: MonomialSum,
+        linear_indices: tuple,
+        contraction_order: tuple = ()) -> Node:
+    """Factor monomial algebra and place ordered contractions.
+
+    Parameters
+    ----------
+    monomial_sum
+        Sum-of-products representation to optimize.
+    linear_indices
+        Free indices identifying argument tabulations.
+    contraction_order
+        Contraction indices, from outermost to innermost stage.
+
+    Returns
+    -------
+    Node
+        Factorized GEM expression.
+
+    Notes
+    -----
+    ``contraction_order`` fixes only reduction-loop placement.  At each
+    level, monomials are partitioned by factors independent of the remaining
+    reductions, the inner reduction is optimized recursively, and COFFEE
+    eliminates sharing in the resulting outer polynomial.  Thus contraction
+    ordering and scalar factorization cooperate without competing planners.
+
+    """
+    monomial_sum = _share_linear_maps(monomial_sum, linear_indices)
+    return _optimise_monomial_sum(
+        monomial_sum, linear_indices, contraction_order)
+
+
+def _optimise_monomial_sum(
+        monomial_sum: MonomialSum,
+        linear_indices: tuple,
+        contraction_order: tuple) -> Node:
+    """Implement recursive monomial and contraction optimization.
+
+    Parameters
+    ----------
+    monomial_sum
+        Sum-of-products representation to optimize.
+    linear_indices
+        Free indices identifying argument tabulations.
+    contraction_order
+        Contraction indices, from outermost to innermost stage.
+
+    Returns
+    -------
+    Node
+        Factorized GEM expression.
+
+    """
+    if contraction_order:
+        grouped = defaultdict(MonomialSum)
+        order = OrderedDict()
+        remaining = frozenset(contraction_order)
+        for monomial in monomial_sum:
+            inner_indices = tuple(index for index in monomial.sum_indices
+                                  if index in remaining)
+            involved = frozenset(inner_indices)
+            inner_atomics = tuple(
+                atomic for atomic in monomial.atomics
+                if involved.intersection(atomic.free_indices))
+            outer_indices = tuple(index for index in monomial.sum_indices
+                                  if index not in remaining)
+            outer_atomics = tuple(
+                atomic for atomic in monomial.atomics
+                if atomic not in inner_atomics)
+            key = outer_indices, outer_atomics
+            order.setdefault(key)
+            grouped[key].add(
+                inner_indices, inner_atomics, monomial.rest)
+
+        outer_sum = MonomialSum()
+        for outer_indices, outer_atomics in order:
+            inner = _optimise_monomial_sum(
+                grouped[(outer_indices, outer_atomics)],
+                linear_indices, contraction_order[1:])
+            outer_sum.add(outer_indices, outer_atomics, inner)
+        monomial_sum = outer_sum
+
     groups = groupby(monomial_sum, key=lambda m: frozenset(m.sum_indices))
-    new_monomials = []
+    optimized = []
     for _, monomials in groups:
-        new_monomials.extend(optimise_monomials(monomials, linear_indices))
-    return monomial_sum_to_expression(new_monomials)
+        old_size = len(monomials) + 1
+        while len(monomials) < old_size:
+            old_size = len(monomials)
+            monomials = optimise_monomials(monomials, linear_indices)
+        optimized.extend(monomials)
+    return monomial_sum_to_expression(optimized)
 
 
 def optimise_monomials(monomials, linear_indices):
@@ -229,24 +457,9 @@ def optimise_monomials(monomials, linear_indices):
     result = [m for m in monomials if not m.atomics]  # skipped monomials
     active_monomials = [m for m in monomials if m.atomics]
 
-    while len(active_monomials) > 0:
-        # Extract a connected component: maximal subset of monomials with intersecting atomics
-        old_size = 0
-        subset = {active_monomials[0]}
-        while len(subset) > old_size:
-            old_size = len(subset)
-            for candidate in active_monomials:
-                if candidate not in subset:
-                    candidate_atomics = frozenset(candidate.atomics)
-                    if any(candidate_atomics.intersection(m.atomics) for m in subset):
-                        subset.add(candidate)
-        connected_monomials = [m for m in active_monomials if m in subset]
-
-        # Optimise the connected component and append to the result
+    for connected_monomials in partition_connected(
+            active_monomials, lambda monomial: monomial.atomics):
         optimal_atomics = find_optimal_atomics(connected_monomials, linear_indices)
         result += factorise_atomics(connected_monomials, optimal_atomics, linear_indices)
 
-        # Discard the connected component
-        active_monomials = [m for m in active_monomials if m not in subset]
-
-    return result
+    return collect_common_rests(result)
diff --git a/gem/contraction.py b/gem/contraction.py
new file mode 100644
index 000000000..8af540223
--- /dev/null
+++ b/gem/contraction.py
@@ -0,0 +1,671 @@
+"""Choose arithmetic- and storage-efficient GEM contraction trees.
+
+This module treats scalar association and index contraction as one
+optimization problem on a tensor-network hypergraph.  Its implementation
+counts rectangular and jagged iteration domains exactly, partitions
+independent networks, and uses subset dynamic programming to place each
+reduction at the earliest legal product subtree.  Plans minimize arithmetic
+work first, then peak live and total materialized storage.
+
+Finite-element rewrites remain in :mod:`gem.optimise`.  The seam here knows
+nothing about element families, tabulations, or quadrature construction; it
+receives contraction indices and scalar factors and returns a GEM expression.
+"""
+
+from collections import defaultdict
+from collections.abc import Callable, Hashable, Iterable
+from functools import lru_cache
+import math
+from typing import TypeVar
+
+import numpy
+
+from gem.gem import (Division, Index, IndexSum, Inverse, Literal,
+                     MathFunction, MaxValue, MinValue, Node, Power, Product,
+                     Solve, Sum)
+from gem.node import traversal
+
+
+T = TypeVar("T")
+
+
+def partition_connected(
+        items: Iterable[T],
+        support: Callable[[T], Iterable[Hashable]]) -> tuple[tuple[T, ...], ...]:
+    """Partition objects connected by transitive support overlap.
+
+    Parameters
+    ----------
+    items
+        Objects to partition.
+    support
+        Function returning the hypergraph vertices incident on an object.
+
+    Returns
+    -------
+    tuple of tuple
+        Connected components in deterministic input order.
+
+    Notes
+    -----
+    Factors joined by contraction indices and monomials joined by common
+    atomics are both incidence hypergraphs.  Traversing their bipartite
+    incidence relation avoids constructing a quadratic pairwise-overlap
+    graph and gives both optimizers the same definition of independence.
+
+    """
+    items = tuple(items)
+    supports = tuple(frozenset(support(item)) for item in items)
+    incidence = defaultdict(list)
+    for position, keys in enumerate(supports):
+        for key in keys:
+            incidence[key].append(position)
+
+    unseen = set(range(len(items)))
+    components = []
+    for seed in range(len(items)):
+        if seed not in unseen:
+            continue
+        component = []
+        pending = [seed]
+        while pending:
+            position = pending.pop()
+            if position not in unseen:
+                continue
+            unseen.remove(position)
+            component.append(position)
+            for key in supports[position]:
+                pending.extend(reversed(incidence[key]))
+        components.append(tuple(items[position]
+                                for position in sorted(component)))
+    return tuple(components)
+
+
+def index_closure(indices: Iterable[Index]) -> frozenset[Index]:
+    """Return indices together with all jagged-loop parents.
+
+    Parameters
+    ----------
+    indices
+        Indices used by an operation.
+
+    Returns
+    -------
+    frozenset of Index
+        The iteration indices required to execute the operation.
+
+    """
+    closure = set(indices)
+    pending = list(closure)
+    while pending:
+        index = pending.pop()
+        for parent in getattr(index, "parents", ()):
+            if parent not in closure:
+                closure.add(parent)
+                pending.append(parent)
+    return frozenset(closure)
+
+
+def _index_components(
+        indices: frozenset[Index]) -> tuple[frozenset[Index], ...]:
+    """Find independent components of an index-parent graph.
+
+    Parameters
+    ----------
+    indices
+        Indices closed under the jagged parent relation.
+
+    Returns
+    -------
+    tuple of frozenset of Index
+        Connected components of the undirected parent graph.
+
+    """
+    neighbours = {index: set() for index in indices}
+    for index in indices:
+        for parent in getattr(index, "parents", ()):
+            neighbours[index].add(parent)
+            neighbours[parent].add(index)
+
+    components = []
+    remaining = set(indices)
+    while remaining:
+        pending = [min(remaining, key=lambda index: index.count)]
+        component = set(pending)
+        while pending:
+            index = pending.pop()
+            new = neighbours[index] - component
+            component.update(new)
+            pending.extend(new)
+        remaining.difference_update(component)
+        components.append(frozenset(component))
+    return tuple(components)
+
+
+def _component_iteration_count(indices: frozenset[Index]) -> int:
+    """Count one connected rectangular or jagged index domain.
+
+    The dynamic-programming state contains only values on the live parent
+    frontier.  Values disappear as soon as no unvisited index depends on
+    them, so equivalent suffixes share one count.
+
+    Parameters
+    ----------
+    indices
+        One connected component of an index-parent graph.
+
+    Returns
+    -------
+    int
+        Number of points in the component domain.
+
+    """
+    parents = {
+        index: frozenset(getattr(index, "parents", ()))
+        for index in indices
+    }
+    ordered = []
+    remaining = set(indices)
+    while remaining:
+        available = sorted(
+            (index for index in remaining
+             if parents[index] <= set(ordered)),
+            key=lambda index: index.count)
+        if not available:
+            raise ValueError("Jagged index parents contain a cycle")
+        ordered.extend(available)
+        remaining.difference_update(available)
+
+    last_use = {
+        index: max(
+            (position for position, child in enumerate(ordered)
+             if index in parents[child]),
+            default=-1,
+        )
+        for index in ordered
+    }
+    frontiers = tuple(
+        tuple(index for index in ordered[:position]
+              if last_use[index] >= position)
+        for position in range(len(ordered) + 1)
+    )
+
+    @lru_cache(maxsize=None)
+    def count(position: int, state: tuple[int, ...]) -> int:
+        if position == len(ordered):
+            return 1
+        values = dict(zip(frontiers[position], state))
+        index = ordered[position]
+        extent = index.iteration_extent(values)
+        next_frontier = frontiers[position + 1]
+        total = 0
+        for value in range(extent):
+            values[index] = value
+            next_state = tuple(values[parent] for parent in next_frontier)
+            total += count(position + 1, next_state)
+        return total
+
+    return count(0, ())
+
+
+@lru_cache(maxsize=1024)
+def _iteration_count(indices: frozenset[Index]) -> int:
+    """Count points in a rectangular or jagged iteration space.
+
+    Independent parent-graph components form a Cartesian product, so their
+    point counts multiply.  Each jagged component is counted by dynamic
+    programming over its live parent frontier.
+
+    Parameters
+    ----------
+    indices
+        Indices on which an operation depends.
+
+    Returns
+    -------
+    int
+        Number of executions of the operation.
+
+    """
+    indices = index_closure(indices)
+    if not indices:
+        return 1
+    if not any(getattr(index, "parents", ()) for index in indices):
+        return int(numpy.prod(
+            [index.extent for index in indices], dtype=int))
+    return math.prod(map(
+        _component_iteration_count, _index_components(indices)))
+
+
+def _storage_count(indices: Iterable[Index]) -> int:
+    """Return the rectangular allocation size for a set of indices.
+
+    Parameters
+    ----------
+    indices
+        Indices retained by an intermediate.
+
+    Returns
+    -------
+    int
+        Number of scalar entries in the intermediate.
+
+    """
+    indices = index_closure(indices)
+    return int(numpy.prod(
+        [index.extent for index in indices], dtype=int))
+
+
+def _operation_count(node: Node) -> int:
+    """Estimate scalar operations performed by one GEM node.
+
+    Parameters
+    ----------
+    node
+        Scalar expression node in a contraction DAG.
+
+    Returns
+    -------
+    int
+        Operations over the node's complete iteration domain.
+
+    """
+    domain = _iteration_count(frozenset(node.free_indices))
+    if isinstance(node, Product):
+        if any(isinstance(child, Literal) and not child.shape
+               and child.value == -1 for child in node.children):
+            return 0
+        return domain
+    if isinstance(node, (Sum, Division, MathFunction, MinValue, MaxValue)):
+        return domain
+    if isinstance(node, Power):
+        _, exponent = node.children
+        if isinstance(exponent, Literal) and not exponent.shape:
+            value = exponent.value
+            if value > 0 and value == math.floor(value):
+                return math.ceil(math.log2(value)) * domain
+        return 5 * domain
+    if isinstance(node, IndexSum):
+        body, = node.children
+        return _iteration_count(frozenset(body.free_indices))
+    if isinstance(node, Inverse):
+        n, _ = node.shape
+        return 2 * n ** 3
+    if isinstance(node, Solve):
+        n, m = node.shape
+        return 2 * n * m + 2 * n ** 3
+    return 0
+
+
+def has_arithmetic(expressions: Iterable[Node]) -> bool:
+    """Does a GEM DAG perform any scalar arithmetic?
+
+    Parameters
+    ----------
+    expressions
+        Roots of a scalar GEM expression DAG.
+
+    Returns
+    -------
+    bool
+        Whether the operation count of :func:`estimate_cost` is nonzero.
+
+    Notes
+    -----
+    A tabulation reference costs nothing to evaluate, so materializing it
+    buys no arithmetic.  Deciding that needs no storage model.
+
+    """
+    return any(map(_operation_count, traversal(tuple(expressions))))
+
+
+def estimate_cost(expressions: Iterable[Node]) -> tuple[int, int, int, int]:
+    """Estimate arithmetic work and contraction storage for a GEM DAG.
+
+    Each structurally shared operation is counted once over the exact
+    rectangular or jagged domain induced by its free indices.  An
+    :class:`IndexSum` contributes one accumulation per point of its body
+    domain.  Storage counts the result domains of contractions, which are the
+    mathematical intermediates exposed to scheduling.
+
+    Parameters
+    ----------
+    expressions
+        Roots of a scalar GEM expression DAG.
+
+    Returns
+    -------
+    tuple of int
+        Operation count, total contraction storage, largest contraction, and
+        expression-node count.
+
+    """
+    nodes = tuple(traversal(tuple(expressions)))
+    sizes = [
+        _storage_count(node.free_indices)
+        for node in nodes if isinstance(node, IndexSum)
+    ]
+    return (
+        sum(map(_operation_count, nodes)),
+        sum(sizes),
+        max(sizes, default=0),
+        len(nodes),
+    )
+
+
+def associate(operator, operands: Iterable[Node]) -> tuple[Node, int]:
+    """Construct a minimum-operation associative expression tree.
+
+    Dynamic programming examines every bipartition of each operand subset.
+    For unusually large expressions a deterministic greedy search bounds
+    compile time.
+
+    Parameters
+    ----------
+    operator
+        Associative binary GEM operator.
+    operands
+        Expressions to combine.
+
+    Returns
+    -------
+    Node
+        Associated GEM expression.
+    int
+        Estimated number of floating-point operations.
+
+    """
+    operands = tuple(operands)
+    if not operands:
+        return operator(), 0
+    if len(operands) == 1:
+        return operands[0], 0
+
+    def combine(left: Node, right: Node) -> tuple[Node, int]:
+        result = operator(left, right)
+        folded = result is left or result is right
+        indices = frozenset(left.free_indices) | frozenset(right.free_indices)
+        return result, 0 if folded else _iteration_count(indices)
+
+    if len(operands) > 8:
+        terms = list(operands)
+        flops = 0
+        while len(terms) > 1:
+            candidates = (
+                (combine(terms[i], terms[j])[1], i, j)
+                for i in range(len(terms))
+                for j in range(i + 1, len(terms)))
+            _, i, j = min(candidates)
+            result, cost = combine(terms[i], terms[j])
+            flops += cost
+            terms = [term for k, term in enumerate(terms)
+                     if k not in (i, j)]
+            terms.append(result)
+        return terms[0], flops
+
+    plans = {
+        1 << position: ((0, 0), operand)
+        for position, operand in enumerate(operands)
+    }
+    for size in range(2, len(operands) + 1):
+        for mask in range(1, 1 << len(operands)):
+            if mask.bit_count() != size:
+                continue
+            anchor = mask & -mask
+            best = None
+            left = (mask - 1) & mask
+            while left:
+                if left & anchor:
+                    right = mask ^ left
+                    if right:
+                        left_score, left_expr = plans[left]
+                        right_score, right_expr = plans[right]
+                        result, cost = combine(left_expr, right_expr)
+                        score = (
+                            left_score[0] + right_score[0] + cost,
+                            max(left_score[1], right_score[1]) + 1,
+                        )
+                        candidate = (score, result)
+                        if best is None or score < best[0]:
+                            best = candidate
+                left = (left - 1) & mask
+            plans[mask] = best
+    score, result = plans[(1 << len(operands)) - 1]
+    return result, score[0]
+
+
+def _ordered_contraction_indices(
+        indices: Iterable[Index],
+        ordering: tuple[Index, ...]) -> tuple[Index, ...]:
+    """Order contractions with jagged parents outside their children.
+
+    Parameters
+    ----------
+    indices
+        Indices to order.
+    ordering
+        Preferred deterministic order.
+
+    Returns
+    -------
+    tuple of Index
+        Legal loop order for an :class:`IndexSum`.
+
+    """
+    indices = frozenset(indices)
+    result = []
+    pending = [index for index in ordering if index in indices]
+    while pending:
+        for position, index in enumerate(pending):
+            parents = set(getattr(index, "parents", ())) & indices
+            if parents <= set(result):
+                result.append(index)
+                pending.pop(position)
+                break
+        else:
+            raise ValueError("Jagged index parents contain a cycle")
+    return tuple(result)
+
+
+@lru_cache(maxsize=4096)
+def _contraction_component(
+        sum_indices: tuple[Index, ...],
+        factors: tuple[Node, ...]) -> Node:
+    """Optimize one connected tensor contraction by subset DP.
+
+    A state is a subset of factors.  A contraction index becomes *closed*
+    when the state contains every factor incident on that index; reducing it
+    at that point is precisely the earliest legal code motion.  The state
+    retains only indices needed by factors outside the subset.  Consequently,
+    its score compares plans lexicographically by FLOPs, peak live
+    intermediate storage, and total materialized storage.
+
+    This dynamic program subsumes scalar reassociation for a connected
+    contraction: each bipartition chooses both a product association and the
+    reductions that become legal there.  :func:`associate` handles expressions
+    for which no contraction-placement decision remains.  Beyond ten factors
+    a deterministic index-at-a-time plan bounds the exponential search.
+
+    Plans are memoized because a loop-ordering search re-plans the same
+    components once per candidate ordering.
+
+    Parameters
+    ----------
+    sum_indices
+        Contracted indices in deterministic order.
+    factors
+        Factors connected through at least one contracted index.
+
+    Returns
+    -------
+    Node
+        Optimized contraction.
+
+    """
+    if len(factors) > 10:
+        terms = list(factors)
+        for index in reversed(sum_indices):
+            contract = [term for term in terms
+                        if index in term.free_indices]
+            if not contract:
+                continue
+            deferred = [term for term in terms
+                        if index not in term.free_indices]
+            product, _ = associate(Product, contract)
+            result = IndexSum(product, (index,))
+            terms = deferred + [result]
+        result, _ = associate(Product, terms)
+        return result
+
+    full_mask = (1 << len(factors)) - 1
+    factor_indices = tuple(
+        index_closure(factor.free_indices) for factor in factors)
+    supports = {
+        index: sum(
+            1 << position
+            for position, indices in enumerate(factor_indices)
+            if index in indices)
+        for index in sum_indices
+    }
+
+    closed_cache = {}
+    live_cache = {}
+
+    def closed(mask: int) -> frozenset[Index]:
+        try:
+            return closed_cache[mask]
+        except KeyError:
+            pass
+        result = frozenset(
+            index for index, support in supports.items()
+            if support and not support & ~mask)
+        closed_cache[mask] = result
+        return result
+
+    def live(mask: int) -> frozenset[Index]:
+        try:
+            return live_cache[mask]
+        except KeyError:
+            pass
+        indices = set().union(*(
+            factor_indices[position]
+            for position in range(len(factors))
+            if mask & (1 << position)))
+        indices.difference_update(closed(mask))
+        result = frozenset(indices)
+        live_cache[mask] = result
+        return result
+
+    def reduce(
+            expression: Node, mask: int,
+            child_closed: frozenset[Index]) -> tuple[Node, int, int]:
+        newly_closed = closed(mask) - child_closed
+        direct = _ordered_contraction_indices(
+            (index for index in newly_closed
+             if index in expression.free_indices),
+            sum_indices)
+        if not direct:
+            return expression, 0, 0
+        extent = _iteration_count(
+            frozenset(expression.free_indices))
+        result = IndexSum(expression, direct)
+        return result, extent, _storage_count(live(mask))
+
+    plans = {}
+    for position, factor in enumerate(factors):
+        mask = 1 << position
+        expression, cost, result_storage = reduce(
+            factor, mask, frozenset())
+        score = (cost, result_storage, result_storage)
+        plans[mask] = (
+            score, expression, result_storage, closed(mask))
+
+    for size in range(2, len(factors) + 1):
+        for mask in range(1, full_mask + 1):
+            if mask.bit_count() != size:
+                continue
+            anchor = mask & -mask
+            best = None
+            left = (mask - 1) & mask
+            while left:
+                if left & anchor:
+                    right = mask ^ left
+                    if right:
+                        left_plan = plans[left]
+                        right_plan = plans[right]
+                        left_score, left_expr, left_storage, left_closed = left_plan
+                        right_score, right_expr, right_storage, right_closed = right_plan
+                        product = Product(left_expr, right_expr)
+                        if product in (left_expr, right_expr):
+                            product_cost = 0
+                        else:
+                            product_cost = _iteration_count(
+                                live(left) | live(right))
+                        expression, reduction_cost, result_storage = reduce(
+                            product, mask, left_closed | right_closed)
+                        flops = (left_score[0] + right_score[0]
+                                 + product_cost + reduction_cost)
+                        live_storage = (left_storage + right_storage
+                                        + result_storage)
+                        peak = max(
+                            left_score[1], right_score[1], live_storage)
+                        total = (left_score[2] + right_score[2]
+                                 + result_storage)
+                        score = (flops, peak, total)
+                        candidate = (
+                            score, expression, result_storage, closed(mask))
+                        if best is None or score < best[0]:
+                            best = candidate
+                left = (left - 1) & mask
+            plans[mask] = best
+    _, expression, _, _ = plans[full_mask]
+    return expression
+
+
+def factor_contraction(
+        sum_indices: Iterable[Index],
+        factors: Iterable[Node]) -> Node:
+    """Plan a normalized tensor contraction.
+
+    Parameters
+    ----------
+    sum_indices
+        Indices to contract in deterministic loop order.
+    factors
+        Scalar factors whose product is contracted.
+
+    Returns
+    -------
+    Node
+        A contraction tree with independent tensor networks reduced
+        separately.
+
+    Notes
+    -----
+    The incidence hypergraph has one vertex for each contracted index and
+    one hyperedge for each factor.  Connected components have no shared
+    contraction and can therefore be planned independently.  Within each
+    component, subset dynamic programming chooses both product association
+    and the earliest legal reductions.
+
+    """
+    sum_indices = tuple(sum_indices)
+    factors = tuple(factors)
+    contraction_set = frozenset(sum_indices)
+    factor_indices = tuple(
+        index_closure(factor.free_indices) & contraction_set
+        for factor in factors)
+    components = partition_connected(
+        range(len(factors)), lambda position: factor_indices[position])
+
+    expressions = []
+    for component in components:
+        component_factors = tuple(factors[position] for position in component)
+        involved = set().union(*(
+            factor_indices[position] for position in component))
+        component_indices = tuple(
+            index for index in sum_indices if index in involved)
+        expressions.append(_contraction_component(
+            component_indices, component_factors))
+    expression, _ = associate(Product, expressions)
+    return expression
diff --git a/gem/flop_count.py b/gem/flop_count.py
index b9595e817..7fed7e9d2 100644
--- a/gem/flop_count.py
+++ b/gem/flop_count.py
@@ -5,6 +5,7 @@
 
 import gem.gem as gem
 import gem.impero as imp
+from contextvars import ContextVar
 from functools import singledispatch
 import numpy
 import math
@@ -25,9 +26,28 @@ def statement_block(tree, temporaries):
 def statement_for(tree, temporaries):
     extent = tree.index.extent
     assert extent is not None
-    child, = tree.children
-    flops = statement(child, temporaries)
-    return flops * extent
+    active_token = _active_indices.set(
+        _active_indices.get() | {tree.index})
+    try:
+        index_values = _index_values.get()
+        if getattr(tree.index, "parents", ()) and all(
+                parent in index_values for parent in tree.index.parents):
+            extent = tree.index.iteration_extent(index_values)
+        child, = tree.children
+        if tree.index in _control_indices.get():
+            flops = 0
+            for value in range(extent):
+                token = _index_values.set(
+                    index_values | {tree.index: value})
+                try:
+                    flops += statement(child, temporaries)
+                finally:
+                    _index_values.reset(token)
+            return flops
+        flops = statement(child, temporaries)
+        return flops * extent
+    finally:
+        _active_indices.reset(active_token)
 
 
 @statement.register(imp.Initialise)
@@ -168,7 +188,37 @@ def flops_solve(expr, temporaries):
 
 @flops.register(gem.ComponentTensor)
 def flops_componenttensor(expr, temporaries):
-    raise ValueError("Not expecting ComponentTensor")
+    body, = expr.children
+    implicit_indices = tuple(
+        index for index in expr.multiindex
+        if index not in _active_indices.get())
+    if not implicit_indices:
+        return expression_flops(body, temporaries)
+    control = _control_indices.get().intersection(implicit_indices)
+    if not control and not any(
+            getattr(index, "parents", ()) for index in implicit_indices):
+        extent = numpy.prod(
+            [index.extent for index in implicit_indices], dtype=int)
+        return extent * expression_flops(body, temporaries)
+
+    def count(position):
+        if position == len(implicit_indices):
+            return expression_flops(body, temporaries)
+        index = implicit_indices[position]
+        values = _index_values.get()
+        extent = index.extent
+        if getattr(index, "parents", ()):
+            extent = index.iteration_extent(values)
+        total = 0
+        for value in range(extent):
+            token = _index_values.set(values | {index: value})
+            try:
+                total += count(position + 1)
+            finally:
+                _index_values.reset(token)
+        return total
+
+    return count(0)
 
 
 def expression_flops(expression, temporaries, top=False):
@@ -192,6 +242,35 @@ def count_flops(impero_c):
     :returns: approximate flop count for the tree.
     """
     try:
-        return statement(impero_c.tree, set(impero_c.temporaries))
+        control_token = _control_indices.set(
+            frozenset(_find_control_indices(impero_c.tree)))
+        index_token = _index_values.set({})
+        active_token = _active_indices.set(frozenset())
+        try:
+            return statement(impero_c.tree, set(impero_c.temporaries))
+        finally:
+            _active_indices.reset(active_token)
+            _index_values.reset(index_token)
+            _control_indices.reset(control_token)
     except (ValueError, NotImplementedError):
         return 0
+
+
+_index_values = ContextVar("flop_count_index_values", default={})
+_active_indices = ContextVar("flop_count_active_indices",
+                             default=frozenset())
+_control_indices = ContextVar("flop_count_control_indices",
+                              default=frozenset())
+
+
+def _find_control_indices(tree):
+    """Find loop indices controlling dependent loop bounds."""
+    result = set()
+    if isinstance(tree, imp.For):
+        if getattr(tree.index, "parents", ()):
+            result.update(tree.index.parents)
+        result.update(_find_control_indices(tree.children[0]))
+    elif isinstance(tree, imp.Block):
+        for child in tree.children:
+            result.update(_find_control_indices(child))
+    return result
diff --git a/gem/gem.py b/gem/gem.py
index 5c9231dc9..d3f972c39 100644
--- a/gem/gem.py
+++ b/gem/gem.py
@@ -16,7 +16,7 @@
 
 from abc import ABCMeta
 from itertools import chain, repeat
-from functools import partial, reduce
+from functools import lru_cache, partial, reduce
 from operator import attrgetter
 from numbers import Integral, Number
 
@@ -32,8 +32,9 @@
            'Variable', 'Sum', 'Product', 'Division', 'FloorDiv', 'Remainder', 'Power',
            'MathFunction', 'MinValue', 'MaxValue', 'Comparison',
            'LogicalNot', 'LogicalAnd', 'LogicalOr', 'Conditional',
-           'Index', 'VariableIndex', 'Indexed', 'ComponentTensor',
-           'IndexSum', 'ListTensor', 'Concatenate', 'Delta', 'OrientationVariableIndex',
+           'Index', 'JaggedIndex', 'VariableIndex', 'Indexed', 'ComponentTensor',
+           'IndexSum', 'ListTensor', 'Concatenate', 'Delta',
+           'OrientationVariableIndex',
            'index_sum', 'partial_indexed', 'reshape', 'view',
            'indices', 'as_gem', 'FlexiblyIndexed',
            'Inverse', 'Solve', 'extract_type', 'uint_type', 'Piecewise']
@@ -311,10 +312,13 @@ def is_equal(self, other):
             return False
         if self.shape != other.shape:
             return False
+        if self.dtype != other.dtype:
+            return False
         return numpy.array_equal(self.array, other.array)
 
     def get_hash(self):
-        return hash((type(self), self.shape, tuple(self.array.flat)))
+        return hash((type(self), self.shape, self.dtype,
+                     tuple(self.array.flat)))
 
     @property
     def value(self):
@@ -614,6 +618,22 @@ def __init__(self, name=None, extent=None):
         self.count = Index._count
         self.extent = extent
 
+    def iteration_extent(self, parent_values: dict) -> int:
+        """Return the loop extent at fixed parent-index values.
+
+        Parameters
+        ----------
+        parent_values
+            Values of indices controlling this index.
+
+        Returns
+        -------
+        int
+            Number of admissible values.
+
+        """
+        return self.extent
+
     def set_extent(self, value):
         # Set extent, check for consistency
         if self.extent is None:
@@ -642,17 +662,66 @@ def __setstate__(self, state):
         self.name, self.extent, self.count = state
 
 
+class JaggedIndex(Index):
+    """Free index whose effective iteration bound depends on the values of
+    other (parent) free indices.
+
+    The iteration bound is ``0 <= i < extent - (p_1 + ... + p_k)`` for
+    parent indices ``p_1, ..., p_k``. The ``extent`` attribute is the static
+    upper bound. Every parent index is zero at this bound. Consumers can
+    treat this as a plain :class:`Index` of extent ``extent``. Expressions
+    indexed by a :class:`JaggedIndex` must evaluate to zero outside the
+    jagged bounds. The jagged bounds only optimize the generated loops.
+
+    Parameters
+    ----------
+    name : str, optional
+        Name of the index.
+    extent : int, optional
+        Static (rectangular) upper bound of the index.
+    parents : tuple of Index
+        The indices whose values reduce the iteration bound.  Loops over
+        this index must nest inside the loops over its parents.
+
+    """
+
+    __slots__ = ('parents',)
+
+    def __init__(self, name: str | None = None, extent: int | None = None,
+                 parents: tuple = ()):
+        super().__init__(name=name, extent=extent)
+        parents = tuple(parents)
+        assert all(isinstance(p, Index) for p in parents)
+        self.parents = parents
+
+    def iteration_extent(self, parent_values: dict) -> int:
+        return self.extent - sum(
+            parent_values[parent] for parent in self.parents)
+
+    def __getstate__(self):
+        return super().__getstate__() + (self.parents,)
+
+    def __setstate__(self, state):
+        super().__setstate__(state[:-1])
+        self.parents = state[-1]
+
+
 class VariableIndex(IndexBase):
     """An index that is constant during a single execution of the
     kernel, but whose value is not known at compile time."""
 
     __slots__ = ('expression',)
 
-    def __init__(self, expression):
+    def __new__(cls, expression):
         assert isinstance(expression, Node)
         assert not expression.shape
         if expression.dtype != uint_type:
             raise ValueError(f"expression.dtype ({expression.dtype}) != uint_type ({uint_type})")
+        if isinstance(expression, Constant):
+            return int(expression.value)
+        return super().__new__(cls)
+
+    def __init__(self, expression):
         self.expression = expression
 
     def __eq__(self, other):
@@ -678,6 +747,15 @@ def __reduce__(self):
         return type(self), (self.expression,)
 
 
+def _index_free_indices(index):
+    """Return the free indices represented by an index expression."""
+    if isinstance(index, Index):
+        return (index,)
+    if isinstance(index, VariableIndex):
+        return index.expression.free_indices
+    return ()
+
+
 class Indexed(Scalar):
     __slots__ = ('children', 'multiindex', 'indirect_children')
     __back__ = ('multiindex',)
@@ -715,7 +793,11 @@ def __new__(cls, aggregate, multiindex):
                 C, = B.children
                 kk = B.multiindex
                 ff = C.free_indices
-                if not any((j in ff) for j in jj):
+                nested = set(chain.from_iterable(
+                    k.expression.free_indices
+                    for k in kk if isinstance(k, VariableIndex)))
+                safe = not set(jj).intersection(set(ff) | nested)
+                if safe:
                     # Only replace indices that are not present in C
                     rep = dict(zip(jj, ii))
                     ll = tuple(rep.get(k, k) for k in kk)
@@ -734,25 +816,14 @@ def __new__(cls, aggregate, multiindex):
         self.multiindex = multiindex
         self.indirect_children = tuple(i.expression for i in self.multiindex if isinstance(i, VariableIndex))
 
-        new_indices = []
-        for i in multiindex:
-            if isinstance(i, Index):
-                new_indices.append(i)
-            elif isinstance(i, VariableIndex):
-                new_indices.extend(i.expression.free_indices)
+        new_indices = tuple(chain.from_iterable(map(_index_free_indices, multiindex)))
         self.free_indices = unique(aggregate.free_indices + tuple(new_indices))
 
         return self
 
     def index_ordering(self):
         """Running indices in the order of indexing in this node."""
-        free_indices = []
-        for i in self.multiindex:
-            if isinstance(i, Index):
-                free_indices.append(i)
-            elif isinstance(i, VariableIndex):
-                free_indices.extend(i.expression.free_indices)
-        return tuple(free_indices)
+        return tuple(chain.from_iterable(map(_index_free_indices, self.multiindex)))
 
 
 class FlexiblyIndexed(Scalar):
@@ -890,6 +961,38 @@ def __new__(cls, expression, multiindex):
         return self
 
 
+def _jagged_layout(multiindex: tuple[Index, ...]) -> tuple:
+    """Return a structural description of a jagged iteration domain."""
+    positions = {}
+    layout = []
+    for position, index in enumerate(multiindex):
+        parents = tuple(positions[parent]
+                        for parent in getattr(index, "parents", ()))
+        layout.append((index.extent, parents))
+        positions[index] = position
+    return tuple(layout)
+
+
+@lru_cache(maxsize=128)
+def _lattice_points(layout: tuple) -> numpy.ndarray:
+    """Enumerate one structural jagged iteration domain."""
+    points = []
+    for alpha in numpy.ndindex(*(extent for extent, _ in layout)):
+        if all(alpha[position] < extent
+               - sum(alpha[parent] for parent in parents)
+               for position, (extent, parents) in enumerate(layout)):
+            points.append(alpha)
+    points = numpy.asarray(points).reshape(len(points), len(layout))
+    points.flags.writeable = False
+    return points
+
+
+def _jagged_lattice(multiindex: tuple[Index, ...]) -> numpy.ndarray:
+    """All lattice points of ``multiindex``'s iteration domain, honouring
+    `JaggedIndex` bounds, as an integer array of shape (npoint, dim)."""
+    return _lattice_points(_jagged_layout(multiindex))
+
+
 class IndexSum(Scalar):
     __slots__ = ('children', 'multiindex')
     __back__ = ('multiindex',)
@@ -901,7 +1004,9 @@ def __new__(cls, summand, multiindex):
             return summand
 
         # Unroll singleton sums
-        unroll = tuple(index for index in multiindex if index.extent <= 1)
+        unroll = tuple(
+            index for index in multiindex
+            if index.extent <= 1 and not getattr(index, "parents", ()))
         if unroll:
             assert numpy.prod([index.extent for index in unroll]) == 1
             summand = Indexed(ComponentTensor(summand, unroll),
@@ -1066,9 +1171,12 @@ def __new__(cls, i, j, dtype=None):
         self = super(Delta, cls).__new__(cls)
         self.i = i
         self.j = j
-        # Set up free indices
-        free_indices = [index for index in (i, j) if isinstance(index, Index)]
-        self.free_indices = tuple(unique(free_indices))
+        # Set up free indices.  A VariableIndex operand is not itself a free
+        # index, but its wrapped expression may be free in other indices
+        # (e.g. a Morton index computed from a lattice multiindex); those
+        # need to propagate here too, exactly as Indexed/FlexiblyIndexed do.
+        self.free_indices = tuple(unique(chain.from_iterable(
+            _index_free_indices(index) for index in (i, j))))
         self._dtype = dtype
         return self
 
diff --git a/gem/impero_utils.py b/gem/impero_utils.py
index 31f9565bb..0e7283638 100644
--- a/gem/impero_utils.py
+++ b/gem/impero_utils.py
@@ -72,8 +72,36 @@ def nonzero(assignment):
 
     get_indices = lambda expr: apply_ordering(expr.free_indices)
 
+    def get_loop_indices(expr: gem.Node) -> tuple[gem.Index, ...]:
+        """Return every explicit loop axis used to evaluate an expression.
+
+        Parameters
+        ----------
+        expr
+            GEM expression being scheduled.
+
+        Returns
+        -------
+        tuple of gem.Index
+            Free indices followed by bound value indices in global loop
+            order.
+
+        Notes
+        -----
+        A ``ComponentTensor`` binds its multi-index in GEM, but evaluating
+        the tensor still executes that index as a value loop.  Exposing the
+        loop to Impero lets several tensor outputs share scalar work within
+        one fused loop instead of materializing that work as arrays.
+
+        """
+        indices = expr.free_indices
+        if isinstance(expr, gem.ComponentTensor):
+            indices = (*indices, *expr.multiindex)
+        return apply_ordering(indices)
+
     # Build operation ordering
-    ops = scheduling.emit_operations(assignments, get_indices, emit_return_accumulate)
+    ops = scheduling.emit_operations(
+        assignments, get_loop_indices, emit_return_accumulate)
 
     # Empty kernel
     if len(ops) == 0:
@@ -83,7 +111,7 @@ def nonzero(assignment):
     ops = inline_temporaries(expressions, ops)
 
     # Build Impero AST
-    tree = make_loop_tree(ops, get_indices)
+    tree = make_loop_tree(ops, get_loop_indices)
 
     # Collect temporaries
     temporaries = collect_temporaries(tree)
@@ -97,9 +125,25 @@ def nonzero(assignment):
 
 def make_prefix_ordering(indices, prefix_ordering):
     """Creates an ordering of ``indices`` which starts with those
-    indices in ``prefix_ordering``."""
+    indices in ``prefix_ordering``.  A `gem.JaggedIndex` is placed after
+    its parents, so that its loop nests inside theirs and the jagged
+    bound can be tightened."""
     # Need to return deterministically ordered indices
-    return tuple(prefix_ordering) + tuple(k for k in indices if k not in prefix_ordering)
+    ordering = tuple(prefix_ordering) + tuple(k for k in indices if k not in prefix_ordering)
+    result = []
+    seen = set()
+
+    def visit(k):
+        if k not in seen:
+            seen.add(k)
+            for parent in getattr(k, 'parents', ()):
+                if parent in ordering:
+                    visit(parent)
+            result.append(k)
+
+    for k in ordering:
+        visit(k)
+    return tuple(result)
 
 
 def make_index_orderer(index_ordering):
@@ -191,11 +235,18 @@ def place_declarations(tree, temporaries, get_indices):
     numbering = {t: n for n, t in enumerate(temporaries)}
     assert len(numbering) == len(temporaries)
 
-    # Collect the total number of temporary references
+    # Collect the total number of temporary references.  Impero is a
+    # tree, so structurally equal subtrees still represent distinct
+    # executions and every occurrence must be visited.  The generic GEM
+    # traversal is DAG-oriented and deliberately skips equal nodes.
     total_refcount = collections.Counter()
-    for node in traversal((tree,)):
+    pending = [tree]
+    while pending:
+        node = pending.pop()
         if isinstance(node, imp.Terminal):
             total_refcount.update(temp_refcount(numbering, node))
+        else:
+            pending.extend(reversed(node.children))
     assert set(total_refcount) == set(temporaries)
 
     # Result
diff --git a/gem/interpreter.py b/gem/interpreter.py
index 13eeb44a2..b2dd609e0 100644
--- a/gem/interpreter.py
+++ b/gem/interpreter.py
@@ -263,8 +263,34 @@ def _evaluate_conditional(e, self):
 def _evaluate_indexed(e, self):
     """Indexing maps shape to free indices"""
     val = self(e.children[0])
-    fids = tuple(i for i in e.multiindex if isinstance(i, gem.Index))
+    variable_indices = {i: self(i.expression) for i in e.multiindex
+                        if isinstance(i, gem.VariableIndex)}
+
+    if any(result.fids for result in variable_indices.values()):
+        # Some variable index depends on free indices: gather entries
+        # one by one over the extent of the free indices.
+        fids = list(val.fids)
+        for i in e.multiindex:
+            new_fids = (i,) if isinstance(i, gem.Index) else \
+                variable_indices[i].fids if isinstance(i, gem.VariableIndex) else ()
+            fids.extend(f for f in new_fids if f not in fids)
+        fids = tuple(fids)
+        out = numpy.empty(tuple(f.extent for f in fids), dtype=val.arr.dtype)
+        for idx in numpy.ndindex(out.shape):
+            env = dict(zip(fids, idx))
+            vidx = [env[f] for f in val.fids]
+            for i in e.multiindex:
+                if isinstance(i, gem.Index):
+                    vidx.append(env[i])
+                elif isinstance(i, gem.VariableIndex):
+                    result = variable_indices[i]
+                    vidx.append(int(result.arr[tuple(env[f] for f in result.fids)]))
+                else:
+                    vidx.append(i)
+            out[idx] = val.arr[tuple(vidx)]
+        return Result(out, fids)
 
+    fids = tuple(i for i in e.multiindex if isinstance(i, gem.Index))
     idx = []
     # First pick up all the existing free indices
     for _ in val.fids:
@@ -275,10 +301,10 @@ def _evaluate_indexed(e, self):
             # Free index, want entire extent
             idx.append(slice(None))
         elif isinstance(i, gem.VariableIndex):
-            # Variable index, evaluate inner expression
-            result, = self(i.expression)
+            # Variable index, constant during kernel execution
+            result = variable_indices[i]
             assert not result.tshape
-            idx.append(result[()])
+            idx.append(int(result.arr[()]))
         else:
             # Fixed index, just pick that value
             idx.append(i)
diff --git a/gem/node.py b/gem/node.py
index 190fe6d40..e5b072df5 100644
--- a/gem/node.py
+++ b/gem/node.py
@@ -37,8 +37,14 @@ def _cons_args(self, children):
 
         Internally used utility function.
         """
-        front_args = (getattr(self, name) for name in self.__front__)
-        back_args = (getattr(self, name) for name in self.__back__)
+        front = self.__front__
+        back = self.__back__
+        if not front and not back:
+            # Operators carry no non-child data, and dominate construction
+            return tuple(children)
+
+        front_args = [getattr(self, name) for name in front]
+        back_args = [getattr(self, name) for name in back]
 
         return (*front_args, *children, *back_args)
 
diff --git a/gem/optimise.py b/gem/optimise.py
index caf254e06..f8ddb8f6f 100644
--- a/gem/optimise.py
+++ b/gem/optimise.py
@@ -1,21 +1,37 @@
-"""A set of routines implementing various transformations on GEM
-expressions."""
+"""Transform GEM expressions while preserving contraction structure.
+
+The contraction optimizer separates three decisions that have different
+mathematical costs.  :func:`associate` reassociates scalar sums and products
+without moving an index contraction.  :func:`sum_factorise` treats a tensor
+contraction as a hypergraph and moves each reduction to the earliest subtree
+containing all factors incident on its index.  A scalar optimizer such as
+COFFEE can then eliminate sharing inside the resulting loops.
+
+This separation makes sum factorization an instance of generalized code
+motion: GEM chooses which index domain computes a value, while COFFEE chooses
+how the scalar value is written.  Contraction plans are ordered first by
+arithmetic work and then by live and total intermediate storage.  The storage
+criteria avoid embedding an architecture-specific cache threshold in this
+symbolic layer.
+"""
 
 from collections import OrderedDict, defaultdict
+from collections.abc import Iterable
 from functools import singledispatch, partial
-from itertools import combinations, permutations, zip_longest
+from itertools import zip_longest
 from numbers import Integral
 
 import numpy
 
+from gem.contraction import associate, factor_contraction, index_closure
 from gem.utils import groupby
 from gem.node import (Memoizer, MemoizerArg, reuse_if_untouched,
                       reuse_if_untouched_arg, traversal)
 from gem.gem import (Node, Failure, Identity, Constant, Literal, Zero,
                      Product, Sum, Comparison, Conditional, Division,
                      Index, VariableIndex, Indexed, FlexiblyIndexed,
-                     IndexSum, ComponentTensor, ListTensor, Delta,
-                     partial_indexed, one)
+                     IndexSum, JaggedIndex, ComponentTensor, ListTensor,
+                     Delta, _jagged_lattice, partial_indexed, one)
 
 
 @singledispatch
@@ -269,6 +285,15 @@ def child(expression):
         elif all(e.j == k and e.i == expr.i for k, e in enumerate(expressions)):
             return expr.reconstruct(expr.i, index)
 
+    if types == {IndexSum}:
+        extents = {tuple(i.extent for i in e.multiindex) for e in expressions}
+        if len(extents) == 1:
+            multiindex = tuple(Index(extent=extent) for extent in extents.pop())
+            summands = [Indexed(ComponentTensor(e.children[0], e.multiindex),
+                                multiindex)
+                        for e in expressions]
+            return IndexSum(_select_expression(summands, index), multiindex)
+
     if len(types) == 1:
         cls, = types
         if cls.__front__ or cls.__back__:
@@ -338,7 +363,6 @@ def substitute(expression, from_, to_):
         to_, = list({delta.i, delta.j} - {from_})
 
         sum_indices.remove(from_)
-
         factors = [substitute(f, from_, to_) for f in factors]
 
         delta_queue = [(f, index)
@@ -348,91 +372,112 @@ def substitute(expression, from_, to_):
     return sum_indices, factors
 
 
-def associate(operator, operands):
-    """Apply associativity rules to construct an operation-minimal expression tree.
-
-    For best performance give factors that have different set of free indices.
-
-    :arg operator: associative binary operator
-    :arg operands: list of operands
+def sum_factorise(
+        sum_indices: Iterable[Index],
+        factors: Iterable[Node],
+        distribute: bool = False) -> Node:
+    """Optimize a tensor contraction using sum factorization.
+
+    The factors form a tensor network whose hyperedges are contraction
+    indices.  Independent connected components are reduced separately.  A
+    subset dynamic program then chooses the contraction tree minimizing
+    arithmetic work, with peak and total intermediate storage as tie-breakers.
+
+    Optional distribution is a competing transformation: it may expose a
+    smaller contraction domain, but it duplicates expression structure.  It
+    is therefore applied only when a summand has strictly smaller contraction
+    support, before the contraction tree is optimized.  Disconnected tensor
+    network components are then planned independently and scalar-associated
+    only after all their reductions have completed.
+
+    Parameters
+    ----------
+    sum_indices
+        Free indices to contract.
+    factors
+        Scalar tensor factors.
+    distribute
+        Split selected sums when doing so exposes contractions over fewer
+        indices.
 
-    :returns: (reduced expression, # of floating-point operations)
-    """
-    if len(operands) > 32:
-        # O(N^3) algorithm
-        raise NotImplementedError("Not expected such a complicated expression!")
-
-    def count(pair):
-        """Operation count to reduce a pair of GEM expressions"""
-        a, b = pair
-        extents = [i.extent for i in set().union(a.free_indices, b.free_indices)]
-        return numpy.prod(extents, dtype=int)
-
-    flops = 0
-    while len(operands) > 1:
-        # Greedy algorithm: choose a pair of operands that are the
-        # cheapest to reduce.
-        a, b = min(combinations(operands, 2), key=count)
-        flops += count((a, b))
-        # Remove chosen factors, append their product
-        operands.remove(a)
-        operands.remove(b)
-        operands.append(operator(a, b))
-    result, = operands
-    return result, flops
-
-
-def sum_factorise(sum_indices, factors):
-    """Optimise a tensor product through sum factorisation.
+    Returns
+    -------
+    Node
+        Optimized GEM expression.
 
-    :arg sum_indices: free indices for contractions
-    :arg factors: product factors
-    :returns: optimised GEM expression
     """
+    sum_indices = tuple(sum_indices)
+    factors = tuple(factors)
     if len(factors) == 0 and len(sum_indices) == 0:
-        # Empty product
         return one
+    if len(sum_indices) == 0:
+        # Without contraction the plan is just an association of the product.
+        expression, _ = associate(Product, factors)
+        return expression
+
+    factor_indices = set().union(*(factor.free_indices for factor in factors))
+    jagged_domain = set().union(*(
+        index_closure((index,))
+        for index in sum_indices if isinstance(index, JaggedIndex)))
+    if jagged_domain - factor_indices:
+        domain_indices = tuple(index for index in sum_indices
+                               if index in jagged_domain)
+        active = index_closure(factor_indices & jagged_domain) \
+            & jagged_domain
+        active_indices = tuple(index for index in domain_indices
+                               if index in active)
+        points = _jagged_lattice(domain_indices)
+        if active_indices:
+            positions = [domain_indices.index(index)
+                         for index in active_indices]
+            multiplicity = numpy.zeros(
+                tuple(index.extent for index in active_indices))
+            numpy.add.at(
+                multiplicity,
+                tuple(points[:, position] for position in positions), 1)
+            domain_factor = Indexed(
+                Literal(multiplicity), active_indices)
+        else:
+            domain_factor = Literal(len(points))
+        factors += (domain_factor,)
+        factor_indices.update(active_indices)
+        marginalised = jagged_domain - active
+        sum_indices = tuple(index for index in sum_indices
+                            if index not in marginalised)
+
+    constant_sum_indices = set(sum_indices) - factor_indices
+    if constant_sum_indices:
+        factors += tuple(Literal(index.extent)
+                         for index in sum_indices
+                         if index in constant_sum_indices)
+        sum_indices = tuple(index for index in sum_indices
+                            if index not in constant_sum_indices)
+
+    if distribute:
+        contraction_indices = frozenset(sum_indices)
+        for position, factor in enumerate(factors):
+            summands = traverse_sum(factor)
+            involved = contraction_indices.intersection(factor.free_indices)
+            if (len(summands) > 1 and involved
+                    and any(any(
+                        contraction_indices.intersection(term.free_indices)
+                        < involved
+                        for term in traverse_product(summand)[1])
+                        for summand in summands)):
+                expressions = []
+                for summand in summands:
+                    extra, summand_factors = traverse_product(summand)
+                    indices = tuple(OrderedDict.fromkeys((*sum_indices, *extra)))
+                    if len(indices) > 6:
+                        break
+                    expressions.append(sum_factorise(
+                        indices,
+                        factors[:position] + tuple(summand_factors)
+                        + factors[position + 1:]))
+                else:
+                    return make_sum(expressions)
 
-    if len(sum_indices) > 6:
-        raise NotImplementedError("Too many indices for sum factorisation!")
-
-    # Form groups by free indices
-    groups = groupby(factors, key=lambda f: f.free_indices)
-    groups = [Product(*terms) for _, terms in groups]
-
-    # Sum factorisation
-    expression = None
-    best_flops = numpy.inf
-
-    # Consider all orderings of contraction indices
-    for ordering in permutations(sum_indices):
-        terms = groups[:]
-        flops = 0
-        # Apply contraction index by index
-        for sum_index in ordering:
-            # Select terms that need to be part of the contraction
-            contract = [t for t in terms if sum_index in t.free_indices]
-            deferred = [t for t in terms if sum_index not in t.free_indices]
-
-            # Optimise associativity
-            product, flops_ = associate(Product, contract)
-            term = IndexSum(product, (sum_index,))
-            flops += flops_ + numpy.prod([i.extent for i in product.free_indices], dtype=int)
-
-            # Replace the contracted terms with the result of the
-            # contraction.
-            terms = deferred + [term]
-
-        # If some contraction indices were independent, then we may
-        # still have several terms at this point.
-        expr, flops_ = associate(Product, terms)
-        flops += flops_
-
-        if flops < best_flops:
-            expression = expr
-            best_flops = flops
-
-    return expression
+    return factor_contraction(sum_indices, factors)
 
 
 def make_sum(summands):
@@ -488,7 +533,8 @@ def applier(expr):
                 else:
                     return expr
         else:
-            applier = lambda expr: expr
+            def applier(expr):
+                return expr
 
         return tuple(renamed), applier
     return partial(_renamer, rename_map, set())
@@ -568,6 +614,156 @@ def traverse_sum(expression, stop_at=None):
     return result
 
 
+def _distribute_sum(expr: Node, predicate=None) -> list[Node]:
+    """Distribute selected sums through products and contractions.
+
+    Parameters
+    ----------
+    expr
+        GEM expression to distribute.
+    predicate
+        Optional predicate selecting operations to distribute.
+
+    Returns
+    -------
+    list of Node
+        Additive terms after distribution.
+
+    Notes
+    -----
+    Memoization uses object identity.  Structurally equal GEM nodes can have
+    deep expression trees, while distribution only needs to reuse actual DAG
+    nodes.
+
+    """
+    if predicate is None:
+        def predicate(node):
+            return True
+
+    results = {}
+    active = {}
+    stack = [(expr, False)]
+    while stack:
+        node, expanded = stack.pop()
+        key = id(node)
+        if key in results:
+            continue
+        if not expanded:
+            stack.append((node, True))
+            stack.extend((c, False) for c in node.children)
+            continue
+        active[key] = predicate(node) or any(
+            active[id(child)] for child in node.children)
+        if active[key] and isinstance(node, (Sum, IndexSum, Product)):
+            if isinstance(node, Sum):
+                results[key] = [
+                    term
+                    for child in node.children
+                    for term in results[id(child)]]
+            elif isinstance(node, IndexSum):
+                body, = node.children
+                results[key] = [
+                    IndexSum(term, tuple(
+                        index for index in node.multiindex
+                        if index in term.free_indices))
+                    for term in results[id(body)]]
+            else:  # Product
+                a, b = node.children
+                ta, tb = results[id(a)], results[id(b)]
+                results[key] = [node] if len(ta) == 1 and len(tb) == 1 \
+                    else [Product(x, y) for x in ta for y in tb]
+        else:
+            results[key] = [node]
+    return results[id(expr)]
+
+
+def preserve_linear_maps(
+        expression: Node,
+        linear_indices: Iterable[Index]) -> tuple[
+            tuple[Node, ...], tuple[Node, ...]]:
+    """Expose multilinear terms and retain each one-axis linear map.
+
+    A sum that depends on one linear index represents a linear map into an
+    argument tabulation. A sum that depends on several linear indices
+    separates multilinear form terms. This function distributes the latter
+    sums and returns the former sums as factors.
+
+    Parameters
+    ----------
+    expression
+        Multilinear GEM expression.
+    linear_indices
+        Free indices identifying the linear axes.
+
+    Returns
+    -------
+    tuple
+        Additive terms and the linear-map factors that they contain.
+
+    Notes
+    -----
+    Polynomial factorization can recover any partial grouping from a fully
+    expanded expression. The map-preserving representation remains useful
+    because it bounds expansion and exposes basis transformation as a
+    separate contraction.
+
+    """
+    linear_indices = frozenset(linear_indices)
+
+    def multilinear_sum(node: Node) -> bool:
+        return isinstance(node, Sum) and len(
+            linear_indices.intersection(node.free_indices)) > 1
+
+    if not any(
+            isinstance(node, Sum)
+            and len(linear_indices.intersection(node.free_indices)) == 1
+            for node in traversal((expression,))):
+        return (expression,), ()
+
+    terms = tuple(_distribute_sum(
+        expression, predicate=multilinear_sum))
+    groups = OrderedDict()
+    for term in terms:
+        _, factors = traverse_product(term)
+        for factor in factors:
+            if (isinstance(factor, Sum)
+                    and len(linear_indices.intersection(
+                        factor.free_indices)) == 1):
+                groups.setdefault(factor)
+
+    if not groups:
+        return (expression,), ()
+    return terms, tuple(groups)
+
+
+def eliminate_deltas(expression):
+    """Cancel contracted deltas without changing other contractions."""
+    replacer = MemoizerArg(filtered_replace_indices)
+    expression = replacer(expression, ())
+    nodes = tuple(traversal((expression,)))
+    contracted = frozenset(
+        index
+        for node in nodes if isinstance(node, IndexSum)
+        for index in node.multiindex)
+
+    def cancellable(node):
+        return isinstance(node, Delta) \
+            and bool({node.i, node.j} & contracted)
+
+    if not any(isinstance(node, Delta) and cancellable(node)
+               for node in nodes):
+        return expression
+
+    terms = []
+    for term in _distribute_sum(expression, predicate=cancellable):
+        indices, factors = traverse_product(term, index_replacer=replacer)
+        indices, factors = delta_elimination(
+            indices, factors, index_replacer=replacer)
+        factors = [replacer(factor, ()) for factor in factors]
+        terms.append(IndexSum(Product(*factors), indices))
+    return make_sum(terms)
+
+
 def contraction(expression, ignore=None):
     """Optimise the contractions of the tensor product at the root of
     the expression, including:
@@ -592,9 +788,11 @@ def contraction(expression, ignore=None):
 
     # Flatten product tree, eliminate deltas, sum factorise
     def rebuild(expression):
+        expression = eliminate_deltas(expression)
         sum_indices, factors = traverse_product(expression, index_replacer=index_replacer)
         sum_indices, factors = delta_elimination(sum_indices, factors, index_replacer=index_replacer)
         factors = [index_replacer(f, ()) for f in factors]
+
         if ignore is not None:
             # TODO: This is a really blunt instrument and one might
             # plausibly want the ignored indices to be contracted on
diff --git a/gem/refactorise.py b/gem/refactorise.py
index 2ca6e4cc0..c9ba02fb9 100644
--- a/gem/refactorise.py
+++ b/gem/refactorise.py
@@ -1,15 +1,25 @@
-"""Data structures and algorithms for generic expansion and
-refactorisation."""
+"""Collect the polynomial structure used by GEM optimizations.
+
+Argument factorization and contraction ordering solve different problems.
+The former identifies the multilinear operands of a finite element form;
+the latter places reduction loops around those operands.  In particular, a
+basis transformation depending on one argument axis remains a linear map
+rather than being expanded into its scalar entries.  Its contractions stay
+visible to GEM, while COFFEE can eliminate sharing between linear maps at
+each reduction level.
+"""
 
-from collections import Counter, OrderedDict, defaultdict, namedtuple
+from collections import Counter, OrderedDict, namedtuple
+from collections.abc import Callable, Iterable
 from functools import singledispatch
 from itertools import product
 from sys import intern
 
 from gem.node import Memoizer, traversal
-from gem.gem import (Node, Conditional, Zero, Product, Sum, Indexed,
+from gem.gem import (Node, Conditional, Zero, Product, Sum, Index, Indexed,
                      ListTensor, one, MathFunction)
-from gem.optimise import (remove_componenttensors, sum_factorise,
+from gem.optimise import (preserve_linear_maps, remove_componenttensors,
+                          sum_factorise,
                           traverse_product, traverse_sum, unroll_indexsum,
                           make_rename_map, make_renamer)
 
@@ -53,7 +63,7 @@ class MonomialSum(object):
     """
     def __init__(self):
         # (unordered sum_indices, unordered atomics) -> rest
-        self.monomials = defaultdict(Zero)
+        self.monomials = {}
 
         # We shall retain ordering for deterministic code generation:
         #
@@ -77,7 +87,8 @@ def add(self, sum_indices, atomics, rest):
         assert isinstance(rest, Node)
 
         key = (sum_indices_set, atomics_set)
-        self.monomials[key] = Sum(self.monomials[key], rest)
+        previous = self.monomials.get(key)
+        self.monomials[key] = rest if previous is None else Sum(previous, rest)
         self.ordering.setdefault(key, (sum_indices, atomics))
 
     def __iter__(self):
@@ -95,7 +106,9 @@ def sum(*args):
             # Optimised implementation: no need to decompose and
             # reconstruct key.
             for key, rest in arg.monomials.items():
-                result.monomials[key] = Sum(result.monomials[key], rest)
+                previous = result.monomials.get(key)
+                result.monomials[key] = \
+                    rest if previous is None else Sum(previous, rest)
             for key, value in arg.ordering.items():
                 result.ordering.setdefault(key, value)
         return result
@@ -235,15 +248,28 @@ def _collect_monomials_mathfunction(expression, self):
 
 
 @_collect_monomials.register(Conditional)
-def _collect_monomials_conditional(expression, self):
-    """Refactorises a conditional expression into a sum-of-products form,
-    pulling only "atomics" out of conditional expressions.
-
-    :arg expression: a GEM expression to refactorise
-    :arg self: function for recursive calls
+def _collect_monomials_conditional(
+        expression: Conditional, self) -> MonomialSum:
+    """Refactorize a compound conditional expression.
+
+    Parameters
+    ----------
+    expression
+        Conditional GEM expression.
+    self
+        Memoized recursive mapper carrying the classifier.
+
+    Returns
+    -------
+    MonomialSum
+        Sum-of-products representation with argument atomics pulled out of
+        the branches.
 
-    :returns: :py:class:`MonomialSum`
     """
+    if self.classifier(expression) != COMPOUND:
+        return _collect_monomials.dispatch(Node.mro()[0])(
+            expression, self)
+
     condition, then, else_ = expression.children
     # Recursively refactorise both branches to `MonomialSum`s
     then_ms = self(then)
@@ -266,20 +292,29 @@ def _collect_monomials_conditional(expression, self):
     return result
 
 
-def collect_monomials(expressions, classifier):
-    """Refactorises expressions into a sum-of-products form, using
-    distributivity rules (i.e. a*(b + c) -> a*b + a*c).  Expansion
-    proceeds until all "compound" expressions are broken up.
+def _collect_monomial_sums(
+        expressions: Iterable[Node],
+        classifier: Callable[[Node], str]) -> list[MonomialSum]:
+    """Collect monomial sums using the supplied node classifier.
 
-    :arg expressions: GEM expressions to refactorise
-    :arg classifier: a function that can classify any GEM expression
-                     as ``ATOMIC``, ``COMPOUND``, or ``OTHER``.  This
-                     classification drives the factorisation.
+    Parameters
+    ----------
+    expressions : iterable of Node
+        GEM expressions to refactorize.
+    classifier : callable
+        Function labeling each node as ``ATOMIC``, ``COMPOUND``, or
+        ``OTHER``.
 
-    :returns: list of :py:class:`MonomialSum`s
+    Returns
+    -------
+    list of MonomialSum
+        Polynomial representations of the expressions.
+
+    Raises
+    ------
+    FactorisationError
+        If a compound expression cannot be expanded.
 
-    :raises FactorisationError: Failed to break up some "compound"
-                                expressions with expansion.
     """
     # Get ComponentTensors out of the way
     expressions = remove_componenttensors(expressions)
@@ -302,3 +337,63 @@ def collect_monomials(expressions, classifier):
     mapper.classifier = classifier
     mapper.rename_map = make_rename_map()
     return list(map(mapper, expressions))
+
+
+def collect_monomials(
+        expressions: Iterable[Node],
+        classifier: Callable[[Node], str],
+        linear_indices: Iterable[Index] = ()) -> list[MonomialSum]:
+    """Collect structure-preserving sum-of-products representations.
+
+    Parameters
+    ----------
+    expressions
+        GEM expressions to refactorize.
+    classifier
+        Function that labels GEM nodes for polynomial collection.
+    linear_indices
+        Free indices identifying the multilinear axes.  Sums depending on
+        exactly one such axis represent linear maps and remain atomic.
+
+    Returns
+    -------
+    list of MonomialSum
+        One polynomial representation for each input expression.
+
+    Notes
+    -----
+    A one-axis sum is a finite element linear operand: examples include a
+    sparse basis transformation and a tensor-product tabulation factor.
+    Preserving it keeps domain-specific basis structure available for
+    contraction optimization.  This does not make the map opaque: its GEM
+    expression remains available for scalar simplification and code motion.
+
+    Sums involving several linear axes separate form monomials and are
+    distributed.  COFFEE subsequently chooses scalar factorizations across
+    the resulting operands.  Keeping these two stages distinct avoids a
+    Cartesian expansion of basis-map entries before loop placement.
+
+    """
+    expressions = tuple(expressions)
+    linear_indices = tuple(linear_indices)
+    if not linear_indices:
+        return _collect_monomial_sums(expressions, classifier)
+
+    result = []
+    for expression in expressions:
+        terms, linear_maps = preserve_linear_maps(
+            expression, linear_indices)
+        if not linear_maps:
+            monomial_sum, = _collect_monomial_sums(
+                (expression,), classifier)
+            result.append(monomial_sum)
+            continue
+
+        map_ids = frozenset(map(id, linear_maps))
+
+        def preserve_map(node: Node) -> str:
+            return ATOMIC if id(node) in map_ids else classifier(node)
+
+        result.append(MonomialSum.sum(*_collect_monomial_sums(
+            terms, preserve_map)))
+    return result
diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py
index 9220dca99..119061b3a 100644
--- a/test/finat/test_zany_mapping.py
+++ b/test/finat/test_zany_mapping.py
@@ -1,11 +1,39 @@
 import FIAT
 import finat
+import gem
 import numpy as np
 import pytest
 import pprint
 
 from gem.interpreter import evaluate
-from finat.physically_mapped import PhysicallyMappedElement
+from gem.node import traversal
+from finat.physically_mapped import MappedTabulation, PhysicallyMappedElement
+
+
+def test_sparse_mapped_tabulation():
+    """Apply a sparse basis map at the cost of its nonzeros."""
+    coefficient = gem.Variable("coefficient", ())
+    matrix = gem.ListTensor(np.asarray([
+        [gem.Literal(1.0), gem.Zero(), coefficient],
+        [gem.Zero(), gem.Literal(1.0), gem.Zero()],
+    ], dtype=object))
+    table_values = np.arange(1.0, 7.0).reshape(3, 2)
+    table = gem.Literal(table_values)
+
+    mapped = MappedTabulation(matrix, {None: table})[None]
+
+    # The three unit entries cost no multiplication, and the one remaining
+    # nonzero costs exactly one. Nothing is selected by a branch.
+    products = [node for node in traversal((mapped,))
+                if isinstance(node, gem.Product)]
+    assert len(products) == 1
+    assert not any(isinstance(node, gem.Conditional)
+                   for node in traversal((mapped,)))
+
+    actual, = evaluate([mapped], {coefficient: np.asarray(2.0)})
+    expected = np.asarray([[1.0, 0.0, 2.0], [0.0, 1.0, 0.0]]) \
+        @ table_values
+    assert np.array_equal(actual.arr, expected)
 
 
 def make_unisolvent_points(element, interior=False):
diff --git a/test/gem/test_simplify.py b/test/gem/test_simplify.py
index fa2242546..ea06d2142 100644
--- a/test/gem/test_simplify.py
+++ b/test/gem/test_simplify.py
@@ -2,6 +2,22 @@
 import gem
 import numpy
 
+from gem import impero
+from gem.coffee import monomial_sum_to_expression
+from gem.flop_count import count_flops
+from gem.impero_utils import (collect_temporaries, compile_gem,
+                              place_declarations)
+from gem.node import traversal
+from gem.interpreter import evaluate
+from gem.optimise import (
+    _distribute_sum,
+    eliminate_deltas,
+    preserve_linear_maps,
+    sum_factorise,
+)
+from gem.refactorise import (ATOMIC, COMPOUND, OTHER,
+                             collect_monomials)
+
 
 @pytest.fixture
 def A():
@@ -55,6 +71,38 @@ def test_componenttensor_from_indexed(A):
     assert A == gem.ComponentTensor(Aij, (i, j))
 
 
+def test_componenttensor_flop_count():
+    i = gem.Index(extent=3)
+    j = gem.Index(extent=3)
+    x = gem.Variable("x", (3,))
+    result = gem.Variable("result", (3,))
+    tensor = gem.ComponentTensor(2 * gem.Indexed(x, (i,)), (i,))
+    expression = gem.Indexed(tensor, (j,))
+    impero_c = compile_gem(
+        [(gem.Indexed(result, (j,)), expression)], (j,))
+
+    assert count_flops(impero_c) == 6
+
+
+def test_componenttensor_sharing_uses_scalar_temporary():
+    """Keep shared work inside a component tensor's value loop."""
+    i = gem.Index(extent=3)
+    j = gem.Index(extent=3)
+    x = gem.Variable("x", (3,))
+    result = gem.Variable("result", (3,))
+    shared = 2 * gem.Indexed(x, (i,))
+    positive = gem.ComponentTensor(shared + 1, (i,))
+    negative = gem.ComponentTensor(shared - 1, (i,))
+    expression = gem.Indexed(positive, (j,)) \
+        + gem.Indexed(negative, (j,))
+
+    impero_c = compile_gem(
+        [(gem.Indexed(result, (j,)), expression)], (j,))
+
+    assert shared in impero_c.temporaries
+    assert impero_c.indices[shared] == ()
+
+
 def test_indexed_transpose(A):
     i, j = gem.indices(2)
     ATij = gem.Indexed(A.T, (i, j))
@@ -84,3 +132,182 @@ def test_flatten_indexsum(A):
     result = gem.IndexSum(gem.IndexSum(Aij, (i,)), (j,))
     expected = gem.IndexSum(Aij, (i, j))
     assert result == expected
+
+
+def test_selective_distribution():
+    a = gem.Variable("a", ())
+    b = gem.Variable("b", ())
+    c = gem.Variable("c", ())
+    i = gem.Index(extent=2)
+    p = gem.Index(extent=1)
+    row = gem.VariableIndex(gem.Indexed(
+        gem.Literal([0], dtype=gem.uint_type), (p,)))
+    delta = gem.Delta(i, row)
+    common = gem.Sum(a, b)
+    expression = gem.Product(common, gem.Sum(c, delta))
+
+    terms = _distribute_sum(
+        expression, predicate=lambda node: isinstance(node, gem.Delta))
+
+    assert len(terms) == 2
+    assert all(common in set(traversal((term,))) for term in terms)
+
+
+def test_preserve_linear_maps_early_exit():
+    """Keep a multilinear sum that contains no separate linear maps."""
+    i = gem.Index(extent=2)
+    j = gem.Index(extent=2)
+    variables = [gem.Variable(f"a{k}", (2, 2)) for k in range(4)]
+    expression = gem.Sum(*(
+        gem.Indexed(variable, (i, j)) for variable in variables))
+
+    terms, linear_maps = preserve_linear_maps(expression, (i, j))
+
+    assert terms == (expression,)
+    assert linear_maps == ()
+
+
+def test_collect_monomials_preserves_linear_maps():
+    """Keep finite element linear maps intact during factorization."""
+    i = gem.Index(extent=2)
+    j = gem.Index(extent=2)
+    left = gem.Sum(
+        gem.Indexed(gem.Literal([1.0, 2.0]), (i,)),
+        gem.Indexed(gem.Literal([3.0, 5.0]), (i,)))
+    right = gem.Sum(
+        gem.Indexed(gem.Literal([7.0, 11.0]), (j,)),
+        gem.Indexed(gem.Literal([13.0, 17.0]), (j,)))
+    expression = left * right
+    linear_indices = frozenset((i, j))
+
+    def classifier(node: gem.Node) -> str:
+        support = linear_indices.intersection(node.free_indices)
+        if not support:
+            return OTHER
+        if isinstance(node, gem.Indexed):
+            return ATOMIC
+        return COMPOUND
+
+    monomial_sum, = collect_monomials(
+        (expression,), classifier, linear_indices)
+
+    monomial, = tuple(monomial_sum)
+    assert frozenset(monomial.atomics) == frozenset((left, right))
+    expected, = evaluate([gem.ComponentTensor(expression, (i, j))])
+    actual, = evaluate([gem.ComponentTensor(
+        monomial_sum_to_expression(monomial_sum), (i, j))])
+    assert numpy.array_equal(actual.arr, expected.arr)
+
+
+def test_constant_variable_index():
+    index = gem.VariableIndex(gem.Literal(1, dtype=gem.uint_type))
+    assert index == 1
+
+
+def test_place_declarations_counts_equal_impero_subtrees():
+    """Equal Impero nodes are distinct occurrences in the loop tree."""
+    expression = gem.Variable("a", ()) * gem.Variable("b", ())
+    tree = impero.Block([
+        impero.Evaluate(expression),
+        impero.Evaluate(expression),
+    ])
+    temporaries = collect_temporaries(tree)
+
+    declare, indices = place_declarations(
+        tree, temporaries, lambda node: node.free_indices)
+
+    assert declare[tree] == []
+    assert indices[expression] == ()
+    assert all(declare[statement] for statement in tree.children)
+
+
+def test_delta_elimination_preserves_indirect_free_index():
+    i = gem.Index(extent=4)
+    k = gem.Index(extent=2)
+    entries = numpy.array([1, 3], dtype=gem.uint_type)
+    indirect = gem.VariableIndex(gem.Indexed(
+        gem.Literal(entries, dtype=gem.uint_type), (k,)))
+    values = gem.Literal([2.0, 3.0, 5.0, 7.0])
+    expression = gem.IndexSum(
+        gem.Delta(i, indirect) * gem.Indexed(values, (i,)), (i,))
+
+    result = eliminate_deltas(expression)
+    assert result.free_indices == (k,)
+    actual, = evaluate([result])
+    assert numpy.array_equal(actual.arr, values.array[entries])
+
+
+def test_sum_factorise_bounded_distribution():
+    indices = tuple(gem.Index(extent=2) for _ in range(6))
+    extra = gem.Index(extent=2)
+
+    def unit(index):
+        return gem.Indexed(gem.Literal(numpy.ones(2)), (index,))
+
+    factor = gem.Sum(gem.IndexSum(unit(indices[0]) * unit(extra), (extra,)),
+                     unit(indices[1]))
+    expression = sum_factorise(
+        indices, [factor, *(unit(index) for index in indices[2:])],
+        distribute=True)
+    value, = evaluate([expression])
+    assert value.arr == 192
+
+
+def test_sum_factorise_distribution():
+    """Preserve rectangular contraction multiplicity after distribution."""
+    indices = tuple(gem.Index(extent=2) for _ in range(2))
+    extra = gem.Index(extent=2)
+
+    def unit(index):
+        return gem.Indexed(gem.Literal(numpy.ones(2)), (index,))
+
+    factor = gem.Sum(gem.IndexSum(unit(indices[0]) * unit(extra), (extra,)),
+                     unit(indices[1]))
+    expression = sum_factorise(indices, [factor], distribute=True)
+    value, = evaluate([expression])
+    assert value.arr == 12
+
+
+def test_sum_factorise_jagged_distribution():
+    """Preserve the joint jagged domain after distribution."""
+    parent = gem.JaggedIndex(extent=3)
+    child = gem.JaggedIndex(extent=3, parents=(parent,))
+
+    def unit(index: gem.Index) -> gem.Node:
+        """Return a unit vector carrying one free index.
+
+        Parameters
+        ----------
+        index
+            Free index of the vector.
+
+        Returns
+        -------
+        gem.Node
+            Indexed unit vector.
+        """
+        return gem.Indexed(gem.Literal(numpy.ones(3)), (index,))
+
+    triangle = numpy.fromfunction(
+        lambda i, j: j < 3 - i, (3, 3), dtype=int)
+    factor = gem.Sum(
+        unit(parent), gem.Indexed(gem.Literal(triangle), (parent, child)))
+    expression = sum_factorise((parent, child), [factor], distribute=True)
+    value, = evaluate([expression])
+    assert value.arr == 12
+
+
+def test_literal_distinguishes_dtypes():
+    """Tell an index literal apart from a value literal.
+
+    An index table holds unsigned integers and a coefficient table holds
+    floats. GEM memoizes on node identity, so the two must not compare
+    equal when they happen to hold the same number.
+    """
+    index = gem.Literal(numpy.uint32(3), dtype=gem.uint_type)
+    value = gem.Literal(3.0)
+
+    assert index.dtype != value.dtype
+    assert index != value
+    assert hash(index) != hash(value)
+    assert {index: "index"}.get(value) is None
diff --git a/test/gem/test_sum_factorise.py b/test/gem/test_sum_factorise.py
new file mode 100644
index 000000000..750a82589
--- /dev/null
+++ b/test/gem/test_sum_factorise.py
@@ -0,0 +1,140 @@
+import numpy
+import pytest
+
+import gem
+from gem.coffee import find_optimal_atomics, optimise_monomial_sum
+from gem.gem import one
+from gem.interpreter import evaluate
+from gem.node import traversal
+from gem.contraction import estimate_cost
+from gem.optimise import sum_factorise
+from gem.refactorise import Monomial, MonomialSum
+
+
+def contraction(nfactors, ndims, extent=2):
+    """Build a product of independent contractions.
+
+    Each factor contracts a table with a coefficient over its own indices.
+    This models one factor of a tensor product coefficient evaluation. No
+    factor carries the indices of another.
+    """
+    numpy.random.seed(0)
+    sum_indices = []
+    factors = []
+    expected = 1.0
+    for _ in range(nfactors):
+        indices = tuple(gem.Index(extent=extent) for _ in range(ndims))
+        table = numpy.random.rand(*(extent,) * ndims)
+        coefficient = numpy.random.rand(*(extent,) * ndims)
+        factors.append(gem.Indexed(gem.Literal(table), indices))
+        factors.append(gem.Indexed(gem.Literal(coefficient), indices))
+        sum_indices.extend(indices)
+        expected *= numpy.sum(table * coefficient)
+    return tuple(sum_indices), factors, expected
+
+
+@pytest.mark.parametrize("nfactors,ndims", [(1, 3), (2, 3), (3, 3), (5, 3), (3, 5)])
+def test_independent_contractions(nfactors, ndims):
+    # Contractions that share no factor are independent, so they are
+    # factorised separately rather than by searching the orderings that
+    # interleave them.  Together they exceed what one exhaustive search
+    # can handle.
+    sum_indices, factors, expected = contraction(nfactors, ndims)
+    assert len(sum_indices) == nfactors * ndims
+
+    expression = sum_factorise(sum_indices, factors)
+    assert expression.free_indices == ()
+
+    result, = evaluate([expression])
+    assert numpy.allclose(result.arr, expected)
+
+
+def test_many_indices_in_one_contraction():
+    indices = tuple(gem.Index(extent=2) for _ in range(7))
+    table = gem.Indexed(gem.Literal(numpy.ones((2,) * 7)), indices)
+    expression = sum_factorise(indices, [table])
+    result, = evaluate([expression])
+    assert result.arr == 2 ** len(indices)
+
+
+def test_optimal_atomics_complete_bipartite():
+    index = gem.Index(extent=3)
+    left = tuple(
+        gem.Indexed(gem.Variable(f"left{i}", (3,)), (index,))
+        for i in range(5))
+    right = tuple(
+        gem.Indexed(gem.Variable(f"right{i}", (3,)), (index,))
+        for i in range(7))
+    monomials = [
+        Monomial((), (a, b), one)
+        for a in left for b in right
+    ]
+
+    selected = find_optimal_atomics(monomials, (index,))
+
+    assert len(selected) == len(left)
+    assert all(any(atomic in monomial.atomics for atomic in selected)
+               for monomial in monomials)
+
+
+def test_share_isomorphic_linear_maps():
+    i = gem.Index(extent=3)
+    j = gem.Index(extent=3)
+    q = gem.Index(extent=2)
+    table = gem.Variable("table", (3, 2))
+    weight = gem.Variable("weight", (3, 2))
+
+    def mapped(index):
+        return 2 * gem.Indexed(table, (index, q)) \
+            + gem.Indexed(weight, (index, q))
+
+    monomial_sum = MonomialSum()
+    monomial_sum.add((), (mapped(i), mapped(j)), one)
+
+    expression = optimise_monomial_sum(monomial_sum, (i, j))
+
+    tensors = [node for node in traversal((expression,))
+               if isinstance(node, gem.ComponentTensor)]
+    tensor, = tensors
+    assert tensor.shape == (3,)
+    assert tensor.free_indices == (q,)
+    accesses = [node for node in traversal((expression,))
+                if isinstance(node, gem.Indexed)
+                and node.children[0] == tensor]
+    assert {access.multiindex for access in accesses} == {(i,), (j,)}
+
+    bindings = {
+        table: numpy.arange(6).reshape(3, 2),
+        weight: numpy.arange(6, 12).reshape(3, 2),
+    }
+    original = mapped(i) * mapped(j)
+    expected, actual = evaluate([original, expression], bindings)
+    assert numpy.array_equal(actual.broadcast(expected.fids), expected.arr)
+
+
+def test_estimate_cost_jagged_contraction():
+    p = gem.JaggedIndex(extent=4)
+    q = gem.JaggedIndex(extent=4, parents=(p,))
+    table = gem.Indexed(gem.Literal(numpy.ones((4, 4))), (p, q))
+    expression = gem.IndexSum(table * table, (p, q))
+
+    operations, storage, largest, _ = estimate_cost((expression,))
+
+    assert operations == 20
+    assert storage == largest == 1
+
+
+def test_estimate_cost_independent_jagged_domains():
+    """Multiply point counts of independent simplex lattices."""
+    p = gem.JaggedIndex(extent=4)
+    q = gem.JaggedIndex(extent=4, parents=(p,))
+    t = gem.JaggedIndex(extent=4, parents=(p, q))
+    r = gem.JaggedIndex(extent=5)
+    s = gem.JaggedIndex(extent=5, parents=(r,))
+    table = gem.Indexed(
+        gem.Literal(numpy.ones((4, 4, 4, 5, 5))), (p, q, t, r, s))
+    expression = gem.IndexSum(table * table, (p, q, t, r, s))
+
+    operations, storage, largest, _ = estimate_cost((expression,))
+
+    assert operations == 2 * 20 * 15