Preserve finite element pullbacks when collecting monomials

[pbrubeck]

TLDR

This is the TSFC half of the FIAT stack that ends at firedrakeproject/fiat#286.

It does two things. TSFC now costs two factorisation plans and keeps the
cheaper one, instead of choosing by a rule. And Loopy now puts adjacent
shared tabulations in one loop, instead of one loop each.

Raviart--Thomas in 3D loses a third of its operations. Q and NCE do not change.

What this does

Cost the two plans (tsfc/spectral.py).

A pullback reaches TSFC as a sum over one argument axis. There are two ways to
handle it.

• Expand it. This exposes scalar factorisation across its entries.
• Keep it. This exposes one tabulation that several argument axes can share.

Neither wins everywhere. So flatten builds both plans and keeps the cheaper
one by estimate_cost. It skips the second plan when has_linear_maps says
there is no map to keep, so the extra work is only paid where it can help.

On a CG degree 3 Laplacian in 2D, keeping the map costs 8,305 operations against
8,125 for expanding it. Only a cost model rejects it there.

Fuse adjacent tabulations (tsfc/loopy.py).

One kept map per ComponentTensor gave one loop each, so a single fused nest
became three loops over the same extent. TSFC now reuses one iname per bound
index and fuses them back.

The iname is reused only between adjacent tabulations. Impero can place
other statements between two tabulations that depend on each other, and one
iname cannot be both inside and outside such a statement. statement_block
therefore resets the memo on any statement that is not a ComponentTensor
Evaluate.

Follow the GEM module split. FIAT moved the cost model to gem.cost and
renamed eliminate_deltas to cancel_nested_deltas.

Effect

Raviart--Thomas is where keeping the pullback pays. On
inner(u, v)*dx + inner(curl(u), curl(v))*dx, main against this PR:

family	dim	degree	flops main	flops PR	AST lines main	AST lines PR
RT	2	3	23,318	18,038	90	81
RT	2	5	253,051	192,564	90	81
RT	3	3	374,589	255,933	136	109
RT	3	5	12,965,694	8,711,463	136	109

Arithmetic falls 32% at degree 3 and 33% at degree 5 in 3D. There are 31% fewer
scalar temporaries and 20% fewer lines of C.

Q and NCE do not change. On a tensor-product cell no sum depends on exactly one
argument index, so the second plan is never built.

On time this reaches parity, not a speedup. Best of four interleaved
repetitions:

case	main (s)	PR (s)	change	main -fno-fast-math	PR	change
RT3 3D	0.014052	0.013975	-0.5%	0.015036	0.013638	-9.3%
RT5 3D	0.068845	0.068857	+0.0%	0.075158	0.068414	-9.0%

PyOP2 compiles with -ffast-math. gcc may then turn a*x + a*y into
a*(x + y) and hoist the sum out of the loop. That is the same saving that
keeping the map makes structurally. Forbid the reassociation and the two
separate: main loses 7.0% and 9.2%, and this branch loses 2.4% and 0.6%.

So the counted 32-33% is real, and the C compiler was already collecting most
of it. What this PR adds is a third fewer operations to schedule, and a
representation that does not depend on gcc choosing to reassociate.

Validation

• tests/tsfc and the zany helmholtz, interpolate and projection regression
  suites: 437 passed
• test_impero_loopy_flop_counts.py covers the ComponentTensor case that
  count_flops used to turn into a silent zero
• test_sum_factorisation.py covers the plan selection

Notes

A DROP BEFORE MERGE commit adds a step to
.github/actions/install/action.yml. It installs the head of the FIAT stack
over the one pyproject.toml resolves from main, so CI exercises both halves
together. Revert it once they land.

AI assistance

Claude Code was used for implementation, benchmarking, and drafting this
section. The human contributor remains responsible for understanding,
validating, and maintaining the changes.