diff --git a/benchmarks/johnson_mercier.py b/benchmarks/johnson_mercier.py
new file mode 100644
index 0000000000..cbe88783ca
--- /dev/null
+++ b/benchmarks/johnson_mercier.py
@@ -0,0 +1,56 @@
+#!/usr/bin/env python
+"""Measure simplex Johnson--Mercier code generation."""
+
+import argparse
+import hashlib
+
+import numpy
+
+from finat.ufl import FiniteElement, VectorElement
+from tsfc import compile_form
+from ufl import FunctionSpace, Mesh, TestFunction, TrialFunction, div, dx, inner
+from ufl.cell import Cell
+
+
+def compile_target(dim: int):
+    """Compile the JM mass-plus-divergence target on a simplex."""
+    cell = Cell(("triangle", "tetrahedron")[dim - 2])
+    mesh = Mesh(VectorElement("CG", cell, 1))
+    element = FiniteElement("Johnson-Mercier", cell, 1)
+    space = FunctionSpace(mesh, element)
+    u = TrialFunction(space)
+    v = TestFunction(space)
+    form = (inner(u, v) + inner(div(u), div(v))) * dx
+    return compile_form(form, parameters={"mode": "spectral"})[0]
+
+
+def temporary_metrics(kernel):
+    """Return counts for the statically allocated Loopy temporaries."""
+    temporaries = kernel.ast.default_entrypoint.temporary_variables.values()
+    sizes = []
+    for temporary in temporaries:
+        if temporary.shape:
+            sizes.append(numpy.prod(temporary.shape, dtype=int))
+    return len(sizes), sum(sizes), max(sizes, default=0)
+
+
+def main():
+    """Print compiler metrics as copyable Markdown."""
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--dims", nargs="+", type=int, default=(2, 3))
+    args = parser.parse_args()
+    print("<!-- generated by benchmarks/johnson_mercier.py -->")
+    print("| dim | flops | temporaries | elements | bytes | largest | AST lines | hash |")
+    print("| ---: | ---: | ---: | ---: | ---: | ---: | ---: | --- |")
+    for dim in args.dims:
+        kernel = compile_target(dim)
+        source = str(kernel.ast)
+        digest = hashlib.sha256(source.encode()).hexdigest()[:12]
+        ntemp, nelem, largest = temporary_metrics(kernel)
+        print(f"| {dim} | {kernel.flop_count:.0f} | "
+              f"{ntemp} | {nelem} | {8 * nelem} | {largest} | "
+              f"{len(source.splitlines())} | {digest} |")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/pyproject.toml b/pyproject.toml
index 063133c4e1..49f8bffd25 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -23,7 +23,8 @@ dependencies = [
   # each Firedrake release to a specific UFL minor version (e.g. 2025.3.x)
   "fenics-ufl @ git+https://github.com/FEniCS/ufl.git@main",
   # TODO RELEASE
-  "firedrake-fiat @ git+https://github.com/firedrakeproject/fiat.git@main",
+  # DROP BEFORE MERGE: pinned to the paired FIAT branch for CI; revert to @main
+  "firedrake-fiat @ git+https://github.com/firedrakeproject/fiat.git@pbrubeck/simplex-sum-factor",
   "h5py>3.12.1",
   "firedrake-rtree>=2026.2.0",
   "immutabledict",
diff --git a/tests/firedrake/regression/test_quadrature.py b/tests/firedrake/regression/test_quadrature.py
index 225a4b244d..bb7189ddc4 100644
--- a/tests/firedrake/regression/test_quadrature.py
+++ b/tests/firedrake/regression/test_quadrature.py
@@ -52,3 +52,42 @@ def test_quadrature_element(mesh, family, mat_type, diagonal):
         a = inner(u, v) * dx
 
     assemble(a, mat_type=mat_type, diagonal=diagonal)
+
+
+@pytest.mark.parametrize("family", ["DG", "CG", "Bernstein"])
+@pytest.mark.parametrize("cell", ["triangle", "tetrahedron"])
+@pytest.mark.parametrize("degree", [1, 3])
+def test_collapsed_quadrature_sum_factorisation(cell, degree, family):
+    """Check sum-factorized residuals and matrices against dense tabulation."""
+    mesh = {"triangle": UnitSquareMesh(2, 2),
+            "tetrahedron": UnitCubeMesh(1, 1, 1)}[cell]
+    variant = None if family == "Bernstein" else "integral"
+    V = FunctionSpace(mesh, family, degree, variant=variant)
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    rg = RandomGenerator(PCG64(seed=0))
+    w = rg.uniform(V, 0, 1)
+
+    # translate_coefficient path (forward transform): residual with a
+    # derivative, mixing both the coefficient and argument sum-factorized
+    # tabulations.
+    L = inner(grad(w), grad(v)) * dx(scheme="canonical")
+    L_collapsed = inner(grad(w), grad(v)) * dx(scheme="collapsed")
+    b = assemble(L)
+    b_collapsed = assemble(L_collapsed)
+    assert np.allclose(b.dat.data, b_collapsed.dat.data, rtol=1e-10, atol=1e-10)
+
+    # translate_argument path (backward transform): mass matrix.
+    a = inner(u, v) * dx(scheme="canonical")
+    a_collapsed = inner(u, v) * dx(scheme="collapsed")
+    M = assemble(a).M.values
+    M_collapsed = assemble(a_collapsed).M.values
+    assert np.allclose(M, M_collapsed, rtol=1e-10, atol=1e-10)
+
+    # Bilinear derivatives exercise two independently transformed argument
+    # lattices.
+    a = inner(grad(u), grad(v)) * dx(scheme="canonical")
+    a_collapsed = inner(grad(u), grad(v)) * dx(scheme="collapsed")
+    K = assemble(a).M.values
+    K_collapsed = assemble(a_collapsed).M.values
+    assert np.allclose(K, K_collapsed, rtol=1e-10, atol=1e-10)
diff --git a/tests/tsfc/test_codegen.py b/tests/tsfc/test_codegen.py
index 8d0bc79655..7dd97e2ba1 100644
--- a/tests/tsfc/test_codegen.py
+++ b/tests/tsfc/test_codegen.py
@@ -1,6 +1,7 @@
+import numpy
 import pytest
 
-from gem import impero_utils
+from gem import gem, impero_utils
 from gem.gem import Index, Indexed, IndexSum, Product, Variable
 
 
@@ -24,6 +25,56 @@ def gencode(expr):
     assert len(gencode(e1).children) == len(gencode(e2).children)
 
 
+def test_jagged_index_codegen(monkeypatch):
+    import islpy as isl
+    import loopy as lp
+    import tsfc.loopy
+
+    # Execute the generated code so we check the numbers, not just the loop bounds
+    monkeypatch.setattr(tsfc.loopy, "target", lp.ExecutableCTarget())
+
+    n = 4
+    extent = n + 1
+    npts = 3
+    ndof = (n + 1) * (n + 2) // 2
+
+    rng = numpy.random.default_rng(7)
+    # Table zero-padded outside the simplex lattice p + q > n, and a
+    # clamped Morton index table for the coefficient gather
+    B = rng.random((extent, extent, npts))
+    morton = numpy.zeros((extent, extent), dtype=gem.uint_type)
+    for p_, q_ in numpy.ndindex(morton.shape):
+        if p_ + q_ > n:
+            B[p_, q_] = 0.0
+        else:
+            morton[p_, q_] = (p_ + q_) * (p_ + q_ + 1) // 2 + q_
+    c = rng.random(ndof)
+
+    i = Index(name="i", extent=npts)
+    p = Index(name="p", extent=extent)
+    q = gem.JaggedIndex(name="q", extent=extent, parents=(p,))
+
+    dof = gem.VariableIndex(Indexed(gem.Literal(morton, dtype=gem.uint_type), (p, q)))
+    integrand = Product(Indexed(Variable("c", (ndof,)), (dof,)),
+                        Indexed(gem.Literal(B), (p, q, i)))
+    expr = IndexSum(integrand, (p, q))
+
+    u = Variable("u", (npts,))
+    impero_c = impero_utils.compile_gem([(Indexed(u, (i,)), expr)], (i, p, q))
+    args = [lp.GlobalArg("u", dtype=numpy.float64, shape=(npts,)),
+            lp.GlobalArg("c", dtype=numpy.float64, shape=(ndof,))]
+    knl, _ = tsfc.loopy.generate(impero_c, args, numpy.float64)
+
+    # The jagged loop must have a domain parametrized by its parent iname
+    assert any(dom.get_var_names(isl.dim_type.param)
+               for dom in knl.default_entrypoint.domains)
+
+    u_out = numpy.zeros(npts)
+    knl(c=c, u=u_out)
+    u_ref = numpy.tensordot(c[morton], B, axes=((0, 1), (0, 1)))
+    assert numpy.allclose(u_out, u_ref, rtol=1e-14)
+
+
 if __name__ == "__main__":
     import os
     import sys
diff --git a/tests/tsfc/test_pickle_gem.py b/tests/tsfc/test_pickle_gem.py
index beb101f912..b68905cb0d 100644
--- a/tests/tsfc/test_pickle_gem.py
+++ b/tests/tsfc/test_pickle_gem.py
@@ -17,6 +17,20 @@ def test_pickle_gem(protocol):
     assert repr(expr) == repr(unpickled)
 
 
+@pytest.mark.parametrize('protocol', range(3))
+def test_pickle_jagged_index(protocol):
+    p = gem.Index(name='p', extent=4)
+    q = gem.JaggedIndex(name='q', extent=4, parents=(p,))
+    expr = gem.IndexSum(gem.Indexed(gem.Variable('A', (4, 4)), (p, q)), (p, q))
+
+    unpickled = pickle.loads(pickle.dumps(expr, protocol))
+    assert repr(expr) == repr(unpickled)
+    up, uq = unpickled.multiindex
+    assert isinstance(uq, gem.JaggedIndex)
+    assert uq.extent == 4
+    assert uq.parents == (up,)
+
+
 @pytest.mark.parametrize('protocol', range(3))
 def test_listtensor(protocol):
     expr = gem.ListTensor([gem.Variable('x', ()), gem.Zero()])
diff --git a/tests/tsfc/test_sum_factorisation.py b/tests/tsfc/test_sum_factorisation.py
index 891cf1c6cc..89488908af 100644
--- a/tests/tsfc/test_sum_factorisation.py
+++ b/tests/tsfc/test_sum_factorisation.py
@@ -3,7 +3,7 @@
 
 from ufl import (Mesh, FunctionSpace, TestFunction, TrialFunction,
                  TensorProductCell, dx, action, interval, triangle,
-                 quadrilateral, curl, dot, div, grad)
+                 tetrahedron, quadrilateral, curl, dot, div, grad, inner)
 from finat.ufl import (FiniteElement, VectorElement, EnrichedElement,
                        TensorProductElement, HCurlElement, HDivElement)
 
@@ -168,7 +168,110 @@ def test_vector_laplace_action(cell, order):
     assert (rates < order).all()
 
 
-if __name__ == "__main__":
-    import os
-    import sys
-    pytest.main(args=[os.path.abspath(__file__)] + sys.argv[1:])
+def simplex_mass(cell, family, degree, scheme='collapsed'):
+    m = Mesh(VectorElement('CG', cell, 1))
+    variant = None if family == "Bernstein" else "integral"
+    V = FunctionSpace(m, FiniteElement(family, cell, degree, variant=variant))
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    return inner(u, v) * dx(scheme=scheme)
+
+
+def simplex_laplacian(cell, family, degree, scheme='collapsed'):
+    m = Mesh(VectorElement('CG', cell, 1))
+    variant = None if family == "Bernstein" else "integral"
+    V = FunctionSpace(m, FiniteElement(family, cell, degree, variant=variant))
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    return inner(grad(u), grad(v)) * dx(scheme=scheme)
+
+
+@pytest.mark.parametrize('family', ["DG", "CG", "Bernstein"])
+@pytest.mark.parametrize(('cell', 'order'), [(triangle, 3), (tetrahedron, 4)])
+def test_simplex_mass_action(cell, family, order):
+    degrees = list(range(3, 9)) if cell is triangle else list(range(3, 8))
+    flops = [count_flops(action(simplex_mass(cell, family, degree)))
+             for degree in degrees]
+    rates = numpy.diff(numpy.log(flops)) / numpy.diff(numpy.log(degrees))
+    assert (rates < order).all()
+
+
+@pytest.mark.parametrize('family', ["DG", "CG", "Bernstein"])
+@pytest.mark.parametrize(('cell', 'order'), [(triangle, 3), (tetrahedron, 4.4)])
+def test_simplex_laplacian_action(cell, family, order):
+    degrees = list(range(3, 9)) if cell is triangle else list(range(3, 8))
+    flops = [count_flops(action(simplex_laplacian(cell, family, degree)))
+             for degree in degrees]
+    rates = numpy.diff(numpy.log(flops)) / numpy.diff(numpy.log(degrees))
+    assert (rates < order).all()
+
+
+@pytest.mark.parametrize('family', ["DG", "CG", "Bernstein"])
+def test_simplex_laplacian_action_compact_codegen(family):
+    form = action(simplex_laplacian(triangle, family, 3))
+    kernel, = compile_form(form, parameters=dict(mode='spectral'))
+    temporaries = kernel.ast.default_entrypoint.temporary_variables
+    assert len(temporaries) < 100
+
+
+def test_bernstein_laplacian_action_compact_literals():
+    degree = 5
+    form = action(simplex_laplacian(tetrahedron, "Bernstein", degree))
+    kernel, = compile_form(form, parameters=dict(mode='spectral'))
+    temporaries = kernel.ast.default_entrypoint.temporary_variables
+    literals = [numpy.asarray(temporary.initializer)
+                for temporary in temporaries.values()
+                if temporary.initializer is not None]
+    lattice_size = (degree + 1) ** 3
+    assert max(literal.size for literal in literals) <= lattice_size
+    assert sum(literal.size for literal in literals) < 10 * lattice_size
+
+
+def test_bernstein_laplacian_bilinear_compact_codegen():
+    import islpy as isl
+    import loopy as lp
+
+    degree = 10
+    collapsed = simplex_laplacian(
+        tetrahedron, "Bernstein", degree, scheme="collapsed")
+    canonical = simplex_laplacian(
+        tetrahedron, "Bernstein", degree, scheme="canonical")
+    collapsed_kernel, = compile_form(
+        collapsed, parameters=dict(mode="spectral"))
+    canonical_kernel, = compile_form(
+        canonical, parameters=dict(mode="spectral"))
+
+    # At this degree the lower asymptotic complexity of sum factorisation
+    # outweighs its setup cost.
+    assert collapsed_kernel.flop_count < canonical_kernel.flop_count
+
+    entrypoint = collapsed_kernel.ast.default_entrypoint
+    assert len(entrypoint.instructions) < 250
+    assert len(entrypoint.temporary_variables) < 250
+    code = lp.generate_code_v2(collapsed_kernel.ast).device_code()
+    assert sum(line.lstrip().startswith("for (")
+               for line in code.splitlines()) < 150
+
+    # A tetrahedral lattice has a loop whose bound depends on two parents.
+    assert max(domain.dim(isl.dim_type.param)
+               for domain in entrypoint.domains) >= 2
+
+
+@pytest.mark.parametrize('family', ["DG", "CG", "Bernstein"])
+@pytest.mark.parametrize(('cell', 'order'), [(triangle, 5), (tetrahedron, 7)])
+def test_simplex_mass_bilinear(cell, family, order):
+    degrees = list(range(3, 9)) if cell is triangle else list(range(3, 8))
+    flops = [count_flops(simplex_mass(cell, family, degree))
+             for degree in degrees]
+    rates = numpy.diff(numpy.log(flops)) / numpy.diff(numpy.log(degrees))
+    assert (rates < order).all()
+
+
+@pytest.mark.parametrize('family', ["DG", "CG", "Bernstein"])
+@pytest.mark.parametrize(('cell', 'order'), [(triangle, 5), (tetrahedron, 7)])
+def test_simplex_laplacian_bilinear(cell, family, order):
+    degrees = list(range(3, 9)) if cell is triangle else list(range(3, 8))
+    flops = [count_flops(simplex_laplacian(cell, family, degree))
+             for degree in degrees]
+    rates = numpy.diff(numpy.log(flops)) / numpy.diff(numpy.log(degrees))
+    assert (rates < order).all()
diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py
index 5d61a916aa..2aeddc5dc4 100644
--- a/tsfc/kernel_interface/common.py
+++ b/tsfc/kernel_interface/common.py
@@ -14,9 +14,10 @@
 import numpy
 from FIAT.reference_element import TensorProductCell
 from finat.cell_tools import max_complex
+from finat.duffy import DuffyElement
 from finat.quadrature import AbstractQuadratureRule
 from gem.node import traversal
-from gem.optimise import constant_fold_zero
+from gem.optimise import constant_fold_zero, unflatten_returns
 from gem.optimise import remove_componenttensors as prune
 from numpy import asarray
 from tsfc import fem
@@ -210,6 +211,8 @@ def compile_gem(self, ctx):
             assignments.extend(mode.flatten(var_reps.items(), ctx['index_cache']))
 
         if assignments:
+            # Rewrite flat FlattenedTensor scatters as jagged lattice loops
+            assignments = unflatten_returns(assignments)
             return_variables, expressions = zip(*assignments)
         else:
             return_variables = []
@@ -344,6 +347,11 @@ def set_quad_rule(params, cell, integral_type, functions):
         scheme = quad_rule
         fiat_cell = as_fiat_cell(cell)
         finat_elements = set(create_element(e) for e in elements if e.family() != "Real")
+        if (scheme == "default" and integral_type == "cell"
+                and any(isinstance(finat_el, DuffyElement)
+                        for finat_el in finat_elements)):
+            # Duffy tabulation requires a collapsed-coordinate point set.
+            scheme = "collapsed"
         fiat_cells = [fiat_cell] + [finat_el.complex for finat_el in finat_elements]
         if any(c.is_macrocell() for c in fiat_cells):
             if len(set(c.get_spatial_dimension() for c in fiat_cells)) > 1:
diff --git a/tsfc/loopy.py b/tsfc/loopy.py
index d4a31a36cb..23bd7e5eaa 100644
--- a/tsfc/loopy.py
+++ b/tsfc/loopy.py
@@ -123,6 +123,7 @@ def __init__(self, target=None):
         self.indices = {}  # indices for declarations and referencing values, from ImperoC
         self.active_indices = {}  # gem index -> pymbolic variable
         self.index_extent = OrderedDict()  # pymbolic variable for indices -> extent
+        self.index_parents = {}  # iname -> parent inames bounding a jagged index
         self.gem_to_pymbolic = {}  # gem node -> pymbolic variable
         self.name_gen = UniqueNameGenerator()
         self.target = target
@@ -257,7 +258,7 @@ def generate(impero_c, args, scalar_type, kernel_name="loopy_kernel", index_name
     instructions, event_name, preamble = profile_insns(kernel_name, instructions, log)
 
     # Create domains
-    domains = create_domains(ctx.index_extent.items())
+    domains = create_domains(ctx.index_extent.items(), ctx.index_parents)
 
     # Create loopy kernel
     knl = lp.make_kernel(
@@ -276,16 +277,23 @@ def generate(impero_c, args, scalar_type, kernel_name="loopy_kernel", index_name
     return knl, event_name
 
 
-def create_domains(indices):
+def create_domains(indices, index_parents=None):
     """ Create ISL domains from indices
 
     :arg indices: iterable of (index_name, extent) pairs
+    :arg index_parents: optional mapping from index_name to a tuple of parent
+        index names; the domain of a jagged index is parametrized by its
+        parents, with upper bound extent minus the sum of the parents.
     :returns: A list of ISL sets representing the iteration domain of the indices."""
 
     domains = []
     for idx, extent in indices:
-        inames = isl.make_zero_and_vars([idx])
-        domains.append(((inames[0].le_set(inames[idx])) & (inames[idx].lt_set(inames[0] + extent))))
+        parents = index_parents.get(idx, ()) if index_parents else ()
+        inames = isl.make_zero_and_vars([idx], parents)
+        bound = inames[0] + extent
+        for parent in parents:
+            bound = bound - inames[parent]
+        domains.append(((inames[0].le_set(inames[idx])) & (inames[idx].lt_set(bound))))
 
     if not domains:
         domains = [isl.BasicSet("[] -> {[]}")]
@@ -316,6 +324,13 @@ def statement_for(tree, ctx):
     assert extent
     idx = ctx.name_gen(ctx.index_names[tree.index])
     ctx.index_extent[idx] = extent
+    if isinstance(tree.index, gem.JaggedIndex) and \
+            all(parent in ctx.active_indices for parent in tree.index.parents):
+        # Tighten the loop bound of a jagged index nested inside its parents.
+        # If a parent loop is not in scope, the rectangular bound `extent`
+        # remains correct: jagged expressions are zero-padded.
+        ctx.index_parents[idx] = tuple(ctx.active_indices[parent].name
+                                       for parent in tree.index.parents)
     with active_indices({tree.index: p.Variable(idx)}, ctx) as ctx_active:
         return statement(tree.children[0], ctx_active)
 
@@ -427,6 +442,12 @@ def _expression_sum(expr, ctx):
     return p.Sum(tuple(expression(c, ctx) for c in expr.children))
 
 
+@_expression.register(gem.FactorisationAtom)
+def _expression_factorisation_atom(expr, ctx):
+    expression_, = expr.children
+    return expression(expression_, ctx)
+
+
 @_expression.register(gem.Division)
 def _expression_division(expr, ctx):
     return p.Quotient(*(expression(c, ctx) for c in expr.children))
diff --git a/tsfc/spectral.py b/tsfc/spectral.py
index 69e471104e..9b06827b55 100644
--- a/tsfc/spectral.py
+++ b/tsfc/spectral.py
@@ -2,14 +2,14 @@
 from functools import partial
 from itertools import chain, zip_longest
 
-from gem.gem import Delta, Indexed, Sum, index_sum, one
-from gem.node import Memoizer, MemoizerArg
+from gem.gem import Delta, FactorisationAtom, Indexed, Sum, index_sum, one
+from gem.node import Memoizer, MemoizerArg, traversal
 from gem.optimise import filtered_replace_indices
 from gem.optimise import delta_elimination as _delta_elimination
 from gem.optimise import replace_division, unroll_indexsum
 from gem.refactorise import ATOMIC, COMPOUND, OTHER, MonomialSum, collect_monomials
 from gem.unconcatenate import unconcatenate
-from gem.coffee import optimise_monomial_sum
+from gem.coffee import sum_factorise_monomial_sum
 from gem.utils import groupby
 
 
@@ -120,16 +120,22 @@ def group_key(pair):
         yield (variable, expression)
 
 
-finalise_options = dict(replace_delta=False)
+finalise_options = dict(replace_delta=True)
 
 
 def classify(argument_indices, expression, delta_inside):
     """Classifier for argument factorisation"""
+    if isinstance(expression, FactorisationAtom):
+        return ATOMIC
     n = len(argument_indices.intersection(expression.free_indices))
     if n == 0:
         return OTHER
     elif n == 1:
-        if isinstance(expression, (Delta, Indexed)) and not delta_inside(expression):
+        mapped = any(isinstance(node, FactorisationAtom)
+                     and node.linear_closure
+                     for node in traversal((expression,)))
+        if (mapped or isinstance(expression, (Delta, Indexed))) \
+                and not delta_inside(expression):
             return ATOMIC
         else:
             return COMPOUND
@@ -162,36 +168,13 @@ def prune(factors):
     variable = factors.pop()
     args = [f for f in factors if f != one]
 
-    assert set(var_indices) == set(variable.free_indices)
+    assert set(var_indices) <= set(variable.free_indices)
+    # A delta may replace a variable index by a contraction index.  That
+    # index now describes a scatter in the assignment, not a sum.
+    sum_indices = [i for i in sum_indices if i not in variable.free_indices]
     return variable, sum_indices, args, rest
 
 
 def sum_factorise(variable, tail_ordering, monomial_sum):
-    if tail_ordering:
-        key_ordering = OrderedDict()
-        sub_monosums = defaultdict(MonomialSum)
-        for sum_indices, atomics, rest in monomial_sum:
-            # Pull out those sum indices that are not contained in the
-            # tail ordering, together with those atomics which do not
-            # share free indices with the tail ordering.
-            #
-            # Based on this, split the monomial sum, then recursively
-            # optimise each sub monomial sum with the first tail index
-            # removed.
-            tail_indices = tuple(i for i in sum_indices if i in tail_ordering)
-            tail_atomics = tuple(a for a in atomics
-                                 if set(tail_indices) & set(a.free_indices))
-            head_indices = tuple(i for i in sum_indices if i not in tail_ordering)
-            head_atomics = tuple(a for a in atomics if a not in tail_atomics)
-            key = (head_indices, head_atomics)
-            key_ordering.setdefault(key)
-            sub_monosums[key].add(tail_indices, tail_atomics, rest)
-        sub_monosums = [(k, sub_monosums[k]) for k in key_ordering]
-
-        monomial_sum = MonomialSum()
-        for (sum_indices, atomics), monosum in sub_monosums:
-            new_rest = sum_factorise(variable, tail_ordering[1:], monosum)
-            monomial_sum.add(sum_indices, atomics, new_rest)
-
-    # Use COFFEE algorithm to optimise the monomial sum
-    return optimise_monomial_sum(monomial_sum, variable.index_ordering())
+    return sum_factorise_monomial_sum(
+        monomial_sum, tuple(tail_ordering), variable.index_ordering())