Sum the four gluon current into the cubic one carrying the same colour factor

[oliviermattelaer]

Each colour structure of the four gluon vertex carries exactly the same colour factor as the diagram obtained by splitting that vertex into two cubic ones. So the quartic current and the cubic current carrying that factor can be summed before the colour algebra is applied, and the amplitude reading the sum gets both contributions from a single call:

      CALL VVV1P0_1(W(3),W(6),GC_10,ZERO,ZERO,W(19))   ! cubic pair
      CALL SUMW_1(W(19),W(12),W(20))                   ! + quartic current W(12)
      CALL VVV1_0(W(4),W(5),W(20),GC_10,AMP(33))       ! one call, both pieces

The option, and the default to decide

set merge_quartic_vertices False   # default today — off
set merge_quartic_vertices speed   # current sums, and the diagram order that allows them
set merge_quartic_vertices slots   # the order that keeps fewest currents alive; no current sums
set merge_quartic_vertices auto    # slots for a gpu backend, speed otherwise, decided per output

True is accepted as a synonym for speed. With no set at all the generated matrix.f and CPPProcess.cc are byte-identical to main for g g > N g, N=2..5, in both backends.

speed and slots have to be set and not output options: they fix the diagram order while the diagrams are generated. auto is the exception, and only because slots is the speed order reversed — see below.

Why two modes

The two objectives pull in opposite directions, and the constraint is sharp: a current sum needs its quartic current emitted before the target amplitude, and that ordering constraint is exactly what costs slots.

g g > 5 g	amplitude calls	slots	wavefunction store
off	7245	268	25.1 kB
speed	6813 (−6%)	259	24.3 kB
slots	7245	199 (−26%)	18.7 kB

On CPU that store is a stack frame and stays cache-resident, so the amplitude calls decide — speed wins, measured. On GPU it is per thread, ~24 kB of local memory at seven gluons, and it is what caps occupancy — so 6% more arithmetic for 23% less memory is plausibly the better trade.

I have not measured that. There is no CUDA device on the machine I ran on. Everything below is CPU/SIMD; the case for slots is an argument, not a measurement.

auto — letting each output choose

The wrinkle behind the default question is that the option is read at generate time while the backend is chosen at output time, and one generation can feed both. auto closes that: it generates the speed order and reverses it at output time, once the backend about to receive the matrix elements is known.

The choice comes from the matrix element exporter, not the output format. output madevent --me_exporter=<gpu backend> hands one _curr_matrix_elements to both exporters in a single export_processes call, so its fortran driver and its gpu matrix elements cannot carry different orders. That is the cudacpp production path — the one case the option exists for. Keying on the format would hand it the cpu order; keying on the me exporter gives it slots (NWAVEFUNCS=54 rather than 78 at six gluons).

Verified byte for byte: auto + standalone reproduces a native speed generation, auto + standalone_mg7 a native slots one, output madevent --me_exporter=standalone_mg7 a native slots madevent, and a session going standalone → standalone_mg7 → standalone reproduces its first output for the third. Processes the seed rule never applied to are untouched.

Two pre-existing behaviours had to be worked around, neither specific to this option:

• An export mutates the diagrams it is given — 345 of 757 vertex leg records on g g > g g g g. Reversing mutated diagrams gives an equivalent but differently numbered result, so the reordering starts from a copy taken before the first export. copy.deepcopy of the amplitudes is not available (it drags the model along and trips assert type(col_obj) != array.array in color_algebra.create_copy); the diagrams alone copy cleanly, 0.011 s at seven gluons.
• An export also drops the marks saying the diagrams came from a seed — afterwards quartic_unroll_tags has 0 entries instead of 405 — so whether an amplitude may be reordered is read before the first export, not at the output that wants to reorder it.

So what should the default be?

Still yours to pick, but auto removes the reason to compromise:

1. Keep False (what this PR does). Opt-in, no behaviour change, nothing to decide until someone runs slots on a device.
2. Default auto. Each backend gets the order it wants, with no per-process choice to make. This is what I would move to once there is a gpu number.
3. Default speed. Best for cpu, but makes a choice on behalf of the gpu backend.

I still recommend (1) for this PR: auto is only worth defaulting to once the slots half of it is backed by a measurement, and there is no CUDA device on the machine I ran on — everything below is CPU/SIMD, and the case for slots remains an argument.

Results (speed against off)

Two runs each, reproducible to ~0.1%. Standalone is the shipped check driver looping SMATRIX; madmatrix is check_sa.exe perf built FPTYPE=d on cppsse4 — the default mixed-precision build rounds both to the same value and would hide any difference.

Speed

	standalone (per call)		madmatrix (evt/s)
g g > g g	11.00 → 11.04 s	−0.4%	875150 → 878724	+0.4%
g g > g g g	34.90 → 34.99 s	−0.3%	74229 → 67566	−9.0%
g g > g g g g	47.77 → 45.25 s	+5.3%	2651 → 2880	+8.6%
g g > 5 g	42.19 → 40.10 s	+4.9%	41.75 → 43.33	+3.8%

Four gluons is a wash — there is nothing to sum, the only quartic vertex is the whole amplitude. Five gluons loses on madmatrix: 7 sums against 7 saved amplitude calls does not cover the extra slots, and madmatrix has no amplitude folds to make up the difference. Six and seven gluons win on both.

Timings were taken in batches with ~1% of drift between them, so each row should be read against the off measured alongside it.

Memory and work — slots / wavefunction calls (+ sums) / amplitude calls

standalone Fortran:

	off	speed
g g > g g g	12 / 33 / 45	19 / 39+7 / 38
g g > g g g g	51 / 111 / 510	78 / 126+30 / 450
g g > 5 g	268 / 898 / 7245	259 / 925+60 / 6813

madmatrix — the slot counts are now the same as Fortran's:

	off	speed
g g > g g g	12 / 45	19 / 38
g g > g g g g	51 / 510	78 / 450
g g > 5 g	268 / 7245	259 / 6813

How it works

1. Seed rule in generation. Unrolling a quartic vertex always yields two cubic vertices sharing the line that replaced it, so a diagram is reachable that way exactly when two of its cubic vertices already share a line. reduce_leglist drops any combination putting two of them on the same line: 1 / 10 / 55 / 385 diagrams instead of 4 / 25 / 220 / 2485.
2. Expansion. The seed is unrolled back to the full diagram set — same diagrams, same count, nothing user-visible moves — but each unrolled diagram keeps the decomposition of the seed it came from, so it is rooted like the diagram it must be summed with.
3. The current sum, where the vertex reading it is an amplitude. Nothing sits above an amplitude, so no consumer can be handed a term it must not have.
4. Canonical leg order, so the same current is never built two ways (see below).
5. Slot recycling and ordering, so the sums come out of the same pool as the wavefunctions and each diagram lands where it keeps lifetimes short.

Bugs fixed on the way

Four, all pre-existing and none of them about this optimisation:

• get_color_amplitudes dropped every merge source from the JAMPs, assuming the caller writes the amplitude sums to put them back. Only the Fortran writer does — so C++/madmatrix/python output with the flag set was silently losing four fifths of the amplitude. It now takes merge_quartic_amplitudes.
• hel_recycle.add_indices could not index a statement-initial AMP( — the pattern ate the character in front of it, and there is none at the start of a line. Latent until something emitted AMP(31) = AMP(31) + AMP(1).
• reuse_outdated_wavefunctions handed a second slot to a wavefunction listed by more than one diagram, leaking the first.
• An antisymmetric current merged by equality. VVV1P0_1(a,b) = −VVV1P0_1(b,a) (measured), but HelasWavefunction.__eq__ compares mothers by sorted number and calls them equal, and export_cpp renumbers wavefunctions by that equality — so the two landed on one slot inside one matrix element and one silently carried the wrong sign. Sidestepped here by building the unrolled vertices with their legs in a canonical order, so the pair never arises; the underlying __eq__/export_cpp mismatch is still there for anything else that produces both orders.

Validation

• |M|² bit-identical at four gluons, 1e-16 at five, 1e-14 at six and seven, in both backends and for all three option values; unchanged for g g > t t~ g g and u u~ > g g g.
• With the option unset, generated code byte-identical to main for N=2..5 in both backends.
• madevent end to end: g g > g g g, 10000 events, three seeds — cross sections agree. The channel weights now see all 510 amplitude pieces instead of 105, because get_amp2_lines skips any diagram with a four point vertex. Measured effect on the integration: none either way — every difference is under two sigma once the seed-to-seed scatter is accounted for, and generate_events wall time swings by a factor four between seeds for the same configuration, so a fixed-work survey comparison over six seeds was added.
• 916 unit tests including 16 new ones; the two failures are pre-existing on main (a test-ordering artefact and a missing scipy).

Not done

• The sum is only taken where the vertex reading it is an amplitude. Deeper it is blocked by counting: expanding the seed makes 340 (seed, choice) instances for the 220 diagrams at six gluons, so a diagram is reached from several seeds, and carrying one decomposition it can be rooted to match at most one of its quartic partners.
• madmatrix has no amplitude-level folds — there is no AMP array to fold into — which is why five gluons still loses there.
• No GPU measurement, so slots -- and therefore the auto default -- is unproven.

Full write-up, including the pitfalls and the measurements that closed off the alternatives, in docs/gluon-quartic-plan.md.

🤖 Generated with Claude Code

[oliviermattelaer]

Follow-up on the diagram order (thanks @oliviermattelaer — the premise of my earlier answer was half wrong).

I claimed the wavefunction count was order-independent. That is true of the wavefunction set (content addressed, wavefunctions[wavefunctions.index(new_wf)]) but not of NWAVEFUNCS: reuse_outdated_wavefunctions is a linear-scan allocator over lifetimes taken in emission order. Emitting every seed first, as expand_seed_diagrams does, is close to the worst case for it.

Four orders built and run. A sum can only be formed when its quartic current is emitted before the target amplitude, which makes this a trade rather than a free win:

NWAVEFUNCS / sums	off	seeds first (shipped)	by quartic count	seed + own unrollings	last discovery
g g > g g g	12	19 / 7	19 / 7	11 / 3	15 / 7
g g > g g g g	51	66 / 30	64 / 30	33 / 0	27 / 30
g g > 5 g	268	290 / 60	314 / 60	245 / 60	219 / 60

Emitting a seed followed by its own unrollings gives the locality but loses the sums — the fully-cubic target is usually claimed by an earlier seed, so the quartic current arrives after it. Placing each diagram at its last discovery instead of its first fixes that and takes both: all the sums, and a slot count below the flag-off baseline (27 vs 51 at six gluons, 219 vs 268 at seven). Fortran stays exact at N=2..5 and is marginally faster (+3.7% / +7.4%).

Not in this PR: it makes g g > g g g g come out wrong in madmatrix (2.43e-04 against 1.59e-04). The generated code has no read-before-write, and a matrix element built directly has no duplicate wavefunction listings — those appear only through the madmatrix exporter path, which is also what forced the allocated guard in the slot-recycling commit. So the reordering looks to be exposing something latent there rather than being wrong itself, but that needs understanding first.

Worth chasing: at six gluons it takes the wavefunction store from 6600 B to 2700 B, under the 5100 B of the unoptimised code, with every current sum kept. Written up in docs/gluon-quartic-plan.md.

[oliviermattelaer]

The missing row is in. g g > 5 g in madmatrix finished compiling (7245 amplitudes in one translation unit takes a while), so the speed table is now complete in both backends:

	standalone (per call)		madmatrix (evt/s)
g g > g g	11.00 → 11.04 s	−0.4%	875150 → 878724	+0.4%
g g > g g g	34.90 → 34.99 s	−0.3%	72359 → 66757	−7.7%
g g > g g g g	47.35 → 45.61 s	+3.7%	2699 → 2784	+3.1%
g g > 5 g	43.05 → 39.98 s	+7.1%	41.75 → 43.33	+3.8%

|M|² = 6.6739867626784836e-07 against 6.6739867626784825e-07, agreeing to 1.6e-15. Wavefunctions 268 → 320, amplitude calls 7245 → 6813, jamp lines 231850 → 218026.

So the shape is the same in both backends: a wash at four gluons, a loss at five in madmatrix, and a growing win from six on. Fortran gains roughly twice what madmatrix does at the same multiplicity, which is the amplitude-fold difference (madmatrix has no AMP array to fold into).

One caveat I checked rather than assumed: this build had been launched before the diagram-ordering experiments, so it could have been generated with an experimental order. I regenerated with the committed code and confirmed CPPProcess.cc is byte-identical, so the numbers do belong to what is in the PR.

[oliviermattelaer]

Root cause of the madmatrix breakage — found, and it is a latent MG5 bug rather than one of this optimisation.

I could not do the two steps in the order asked (fix the duplication, then apply the reordering), because the "duplication" turned out not to be a duplication at all. The chain:

1. The reconstruction can build the same cubic current with its two mothers in either order — the colliding pairs are interaction 3, colour_key 0, mothers [7,3] against [3,7]. sorted_mothers leaves them alone: for two identical gluons its key ties and the sort is stable.
2. VVV1P0_1 is antisymmetric under exchanging its two inputs. Measured, not inferred — with a fixed pair of wavefunctions, VVV1P0_1(a,b) + VVV1P0_1(b,a) = 0 exactly. So the two objects are negatives of each other and write out calls differing by a sign.
3. HelasWavefunction.__eq__ compares mothers by sorted number and so calls them equal — its docstring says "the number for this wavefunction, the pdg code, and the interaction id are irrelevant".
4. export_cpp.generate_process_files renumbers wavefunctions by that equality, to share them between matrix elements. The two therefore land on one number and one slot inside a single matrix element, and whichever is written last wins — with the wrong sign for the other.

The Fortran writer never hits this because it does not renumber by equality. With the committed diagram order the collisions happen not to matter; the last-discovery order moves them somewhere they do.

Ruled out along the way, each by measurement rather than argument: the reordering itself (exact with the sums off), slot corruption (no read-before-write, no sum aliasing its own input, nwf sufficient), coloramps.h, the merge_quartic_amplitudes=False path (the Fortran writer forced into madmatrix's semantics is exact), and the duplicate listings on their own.

I have not applied the reordering, because the fix is upstream of this PR and the direction is yours to choose:

• make __eq__ compare mothers in order rather than sorted — arguably strictly safer, since for any vertex whose particles differ sorted_mothers fixes the order anyway, so only the identical-particle case changes. But it alters the wavefunction CSE everywhere and needs validating with the flag off across the suite;
• or make the reconstruction produce only one of the two orders. Canonicalising the pair inside split_quartic_vertex is not sufficient — I tried it and the collision count stayed at 20, because the other copy can come from a vertex the reconstruction did not build.

Shipped in the meantime (65c457fd6): is_unrolled_pair now checks that the quartic's colour structure is the one separating the cubic's inner pair. It matched on the lines alone before, and all three colour structures of a quartic vertex share the same lines and the same outgoing line — that was only not biting by luck. No change to any result on the committed ordering.

[oliviermattelaer]

You were right that the wavefunction count should not have gone up — chasing that fixed the madmatrix blocker as well.

At six gluons the flag took the wavefunctions from 111 to 146. Breaking that down: 20 of the 35 extra are the same current built with its two mothers in the opposite order. sorted_mothers leaves them alone (for two identical gluons its key ties and the sort is stable), so they stay two objects — and since VVV1P0_1 is antisymmetric they are negatives of each other and cannot be shared.

Taking the legs of both unrolled vertices in a canonical order removes every one:

	wavefunctions	order-flipped twins
g g > g g g	33 → 39	0
g g > g g g g	111 → 126 (was 146)	20 → 0
g g > 5 g	898 → 925 (was 955)	30 → 0

And those twins were exactly what broke madmatrix. They are what __eq__ calls equal while they differ by a sign, and export_cpp renumbers wavefunctions by that equality — so they landed on one slot inside one matrix element and one silently carried the wrong sign. With none left there is nothing to collide, so the last-discovery diagram order is now in, with every current sum surviving.

Where that leaves the numbers

standalone Fortran, slots / wavefunction calls + sums / amplitude calls:

	flag off	before this work	at PR open	now
g g > g g g	12 / 33 / 45	12 / 33 / 45	19 / 39+7 / 38	19 / 39+7 / 38
g g > g g g g	51 / 111 / 510	51 / 111 / 510	66 / 146+30 / 450	78 / 126+30 / 450
g g > 5 g	268 / 898 / 7245	268 / 898 / 7245	290 / 955+60 / 6813	259 / 925+60 / 6813

madmatrix, slots / amplitude calls — no "before" numbers, because with the flag set the C++ output was silently losing 405 of its 510 amplitudes at six gluons (get_color_amplitudes dropped every merge source and only the Fortran writer wrote them back; fixed in this PR):

	flag off	at PR open	now
g g > g g g	12 / 45	19 / 38	19 / 38
g g > g g g g	51 / 510	86 / 450	78 / 450
g g > 5 g	268 / 7245	320 / 6813	259 / 6813

At seven gluons the slot count is now below the unoptimised build in both backends: 259 against 268.

Speed:

	off	at PR open	now
fortran g g > g g g g	47.77 s	45.61 s	45.25 s (+5.3%)
fortran g g > 5 g	42.19 s	39.98 s	40.10 s (+4.9%)
madmatrix g g > g g g	74229	66757	67566 (−9.0%)
madmatrix g g > g g g g	2651	2784	2880 (+8.6%)

Six gluons roughly doubles in madmatrix (+3.1% → +8.6%) and improves in Fortran (+3.7% → +5.3%). Five gluons still loses in madmatrix — 7 sums against 7 saved amplitude calls does not cover the extra slots, and madmatrix has no amplitude folds to make up the difference.

|M|² stays exact: bit-identical at four gluons, 1e-16 at five, 1e-14 at six and seven, in both backends. With the flag off matrix.f and CPPProcess.cc are byte-identical to before. 913 unit tests, same two pre-existing failures.

Note the timings are taken in batches with about 1% of drift between them, so each column should be read against the off measured alongside it — that is why g g > 5 g reads +4.9% now against +7.1% earlier despite the faster absolute time.

[oliviermattelaer]

Full sweep of both series, all three modes, on speed / memory / generation time. This supersedes the standalone timings in the description, which were single runs on g g > N g only.

Re-measured after the AMP array was made recyclable (see the follow-up comment below) — the AMP column is now the entry count, not the amplitude count, and the two differ: at seven gluons slots runs 7245 amplitude calls through 946 entries.

Timings are the minimum of five runs of the shipped check driver looping SMATRIX. The minimum matters: a single run carries ~3% of noise and most of the effects here are smaller than that. The floor of the method is 0.6%, measured on g g > t t~, where off and speed produce a byte-identical matrix.f and still time 4.000 against 3.988 µs. |M|² agrees to 8.5e-15 or better on every row.

W slots is NWAVEFUNCS, the length of TYPE(ALOHA) W(NWAVEFUNCS) — not the number of wavefunctions computed, since reuse_outdated_wavefunctions frees an entry as soon as its last reader has run. One entry is 104 bytes, measured with storage_size; one AMP entry is 16 bytes.

g g > N g

process	mode	generate	matrix.f	matrix.o	W slots	W array	AMP entries	amps computed	AMP array	per call	speed
g g > 2g	off	1.8 s	27 kB	16 kB	5	0.5 kB	6	6	0.1 kB	5.48 us	-
	speed	1.5 s	27 kB	17 kB	5	0.5 kB	4	6	0.1 kB	5.52 us	-1%
	slots	1.6 s	27 kB	17 kB	5	0.5 kB	4	6	0.1 kB	5.55 us	-1%
g g > 3g	off	2.1 s	40 kB	28 kB	12	1.2 kB	45	45	0.7 kB	87.50 us	-
	speed	1.8 s	39 kB	29 kB	19	1.9 kB	24	38	0.4 kB	90.25 us	-3%
	slots	2.4 s	40 kB	29 kB	12	1.2 kB	16	45	0.2 kB	93.00 us	-6%
g g > 4g	off	2.8 s	181 kB	141 kB	51	5.2 kB	510	510	8.0 kB	2.37 ms	-
	speed	2.4 s	165 kB	146 kB	78	7.9 kB	316	450	4.9 kB	2.29 ms	+4%
	slots	2.4 s	168 kB	152 kB	54	5.5 kB	106	510	1.7 kB	2.41 ms	-2%
g g > 5g	off	35.2 s	3.3 MB	5.8 MB	268	27.2 kB	7245	7245	113.2 kB	141.00 ms	-
	speed	16.6 s	2.5 MB	3.4 MB	259	26.3 kB	5869	6813	91.7 kB	132.75 ms	+6%
	slots	18.4 s	2.5 MB	3.4 MB	199	20.2 kB	946	7245	14.8 kB	134.25 ms	+5%

g g > t t~ N g

process	mode	generate	matrix.f	matrix.o	W slots	W array	AMP entries	amps computed	AMP array	per call	speed
g g > t t~	off	1.5 s	26 kB	16 kB	5	0.5 kB	3	3	0.0 kB	4.00 us	-
	speed	1.5 s	26 kB	16 kB	5	0.5 kB	3	3	0.0 kB	3.99 us	+0%
	slots	2.0 s	26 kB	16 kB	5	0.5 kB	3	3	0.0 kB	4.09 us	-2%
g g > t t~ g	off	1.8 s	31 kB	20 kB	12	1.2 kB	18	18	0.3 kB	30.00 us	-
	speed	1.6 s	31 kB	20 kB	12	1.2 kB	15	15	0.2 kB	29.58 us	+1%
	slots	1.7 s	31 kB	21 kB	12	1.2 kB	15	18	0.2 kB	30.33 us	-1%
g g > t t~ 2g	off	2.6 s	68 kB	51 kB	26	2.6 kB	159	159	2.5 kB	377.00 us	-
	speed	2.4 s	65 kB	49 kB	35	3.6 kB	109	126	1.7 kB	343.00 us	+9%
	slots	2.4 s	66 kB	53 kB	29	2.9 kB	106	159	1.7 kB	391.00 us	-4%
g g > t t~ 3g	off	6.2 s	576 kB	479 kB	121	12.3 kB	1890	1890	29.5 kB	9.63 ms	-
	speed	4.9 s	463 kB	408 kB	213	21.6 kB	1159	1551	18.1 kB	8.97 ms	+7%
	slots	5.0 s	493 kB	466 kB	141	14.3 kB	946	1890	14.8 kB	9.93 ms	-3%

What the columns show

The AMP array is recycled, and it is where slots wins. It used to be declared at the full diagram count in every mode — 113 kB at seven gluons, with speed leaving 432 entries written by nobody. Now slots runs 7245 amplitude calls through 946 entries at seven gluons, 14.8 kB rather than 113.2, while speed only reaches 5869 because of the order it emits them in.

slots mode computes every amplitude — amps computed equals the off-mode amplitude count on all its rows. With no current sums nothing is skipped, so it does the baseline's amplitude work plus the folds, buying a shorter JAMP block and much smaller arrays rather than less arithmetic. That is why it is slower than off on CPU almost everywhere. speed is the opposite — it skips amplitudes outright (450 of 510, 6813 of 7245, 1551 of 1890), which is where its 4–9% comes from, and pays in W slots (121 → 213 at t t~ 3g).

What it is actually good for

Generation time and code size, more than speed. g g > 5 g generates in 16.6 s rather than 35.2 s (−53%, reproducible: 385 seed diagrams unrolled beats 2485 generated), its matrix.o goes 5.8 → 3.4 MB and its matrix.f 3.3 → 2.5 MB. Runtime is 4–9% above six particles and nothing at all below. The t t~ series gains more than the pure gluon one at equal particle count: +9% at t t~ 2g against +4% at 4g. Peak RSS is flat except at seven gluons, because the wavefunction store is a stack frame and the code image dominates.

Full write-up in docs/gluon-quartic-plan.md, sections "Full sweep" and "Recycling the AMP array".

How far it generalises

Both series above are gluon-rich, so a third process was measured as a control: u u~ > z g g g g — seven legs like g g > t t~ 3g, but with a Z and a quark line, so most of its diagrams have no four-gluon vertex at all. (q is not a defined multiparticle, hence u u~.)

process	mode	generate	matrix.f	matrix.o	W slots	AMP entries	amps computed	per call	speed	total/call
g g > t t~ 3g	off	6.2 s	576 kB	479 kB	121	1890	1890	9.63 ms	–	41.8 kB
	speed	4.9 s	463 kB	408 kB	213	1159	1551	8.97 ms	+7%	39.7 kB (−5%)
	slots	5.0 s	493 kB	466 kB	141	946	1890	9.93 ms	−3%	29.1 kB (−30%)
u u~ > z 4g	off	2.9 s	135 kB	115 kB	76	516	516	1.58 ms	–	15.8 kB
	speed	2.3 s	127 kB	110 kB	85	391	450	1.51 ms	+4%	14.7 kB (−7%)
	slots	2.3 s	132 kB	117 kB	84	384	516	1.65 ms	−5%	14.5 kB (−8%)

The payoff tracks the quartic fraction, exactly. In slots mode the AMP entry count is total amplitudes − merge sources, to within one, on every process measured:

	amplitudes	merge sources	AMP in slots	saving
g g > 5 g	7245	6300 (87%)	946	−87%
g g > g g g g	510	405 (79%)	106	−79%
g g > t t~ 3g	1890	945 (50%)	946	−50%
u u~ > z 4g	516	132 (26%)	384	−26%

So u u~ > z 4g is the weakest case measured, and predictably — there is simply little to merge. Runtime follows at +4% rather than +7 to +9%, and the working set at −8% rather than −30 or −75%.

It is also the one process where slots does not reduce the wavefunctions either — 84 against off's 76, worse — breaking the pattern from the gluon-rich processes, where reversing the order always recovered them. There slots is the worst of the three: 5% slower than off and larger in W, for an AMP count it barely wins over speed (384 vs 391). speed still behaves: +4%, generation down 21%, smaller source and object — which is what the auto gate picks at seven legs anyway.

[oliviermattelaer]

Two follow-ups from the sweep: a multiplicity gate on auto, and the AMP array is now recycled.

auto gates on the multiplicity

The sweep puts the turnover of the full merging at six external legs, so auto now only takes the seed rule at merge_quartic_min_legs (6) or above. speed and slots asked for by name stay unconditional.

What is left below the threshold is not nothing, and it was worth measuring rather than assuming: the amplitude merges do not need the seed rule — they come from the colour algebra through unroll_quartic_vertices — so they still apply, and they shrink the JAMP block without touching the wavefunctions.

g g > g g g	slots	JAMP temporaries	per call
off	12	72	87.75 µs
auto (merges only)	12	42	84.25 µs, +4.0%
speed (full)	19	42	87.25 µs, −1%

So below the threshold auto beats both — it keeps the JAMP fold, which is free, and drops the reordering, which is what costs the seven extra slots. At four and five legs elsewhere it is neutral and never negative, at an unchanged slot count. At six legs and above auto generates a byte-identical matrix.f to speed.

The AMP array is recycled

AMP was declared COMPLEX*16 AMP(NGRAPHS) at the full diagram count in every mode — 113 kB at seven gluons, with speed leaving 432 entries written by nobody. get_amplitude_slots now does for AMP what reuse_outdated_wavefunctions does for W. The enabling change is where the merges are written: they used to be emitted in one block at the very end, keeping every source alive to the end, and each is now written as soon as both of its amplitudes exist.

process	mode	AMP entries	AMP	W slots	W	total per call
g g > g g g g	off	510	8.0 kB	51	5.2 kB	13.1 kB
	speed	316	4.9 kB	78	7.9 kB	12.9 kB (−2%)
	slots	106	1.7 kB	54	5.5 kB	7.1 kB (−46%)
g g > t t~ g g	off	159	2.5 kB	26	2.6 kB	5.1 kB
	speed	109	1.7 kB	35	3.6 kB	5.3 kB (+3%)
	slots	106	1.7 kB	29	2.9 kB	4.6 kB (−10%)
g g > 5 g	off	7245	113.2 kB	268	27.2 kB	140.4 kB
	speed	5869	91.7 kB	259	26.3 kB	118.0 kB (−16%)
	slots	946	14.8 kB	199	20.2 kB	35.0 kB (−75%)

slots reaches the floor and speed does not, and the reason is the diagram order rather than the allocator. Only (2n−5)!! amplitudes are read by the JAMPs — 105 at six gluons, 945 at seven — and the rest are merge sources that could share a handful of entries. speed emits every seed before its unrollings, so a source is born early and its target arrives late and nothing can be reclaimed in between: 5869 rather than 945. Reversing that order puts each source beside its target, so slots lands one above the floor.

That changes the slots case substantially. It used to buy 23% of the wavefunction store for 6% more arithmetic; it now buys 75% of the whole per-call working set at seven gluons — which is the number that matters on a GPU, where this is per thread.

It buys no time, and I want to be clear about that: speed 2247 → 2260 µs and slots 2347 → 2373 µs across the change at six gluons, both inside the 0.6% floor and if anything marginally the wrong way. The arrays were already cache-resident at this size. This is a memory optimisation, not a speed one.

Correctness

NGRAPHS only ever dimensioned AMP inside matrix.f, so it becomes the entry count while ngraphs.inc keeps the diagram count. Everything reading AMP afterwards goes through the same map — the JAMPs via map_color_amplitudes, AMP2 via get_amplitude_slot_map. AMP2 was the one worth checking, since multichannel reads individual amplitudes: it reads only merge targets, which are exactly the entries that stay put. Verified on a madevent output at six gluons — 316 written, 316 read, nothing read that is never written.

|M|² unchanged on every row; 918 unit tests, two new, same two pre-existing failures.