SecRed: A program for sector-wise IBP reduction of Feynman integrals

Program Summary

Program Title: SecRed 1.0

CPC Library link to program files: to be added by Technical Editor

Developer’s repository link: https://github.com/chiyutuci/secred

Licensing provisions: GPL-3.0

Programming language: C++

External routines/libraries: FLINT [1] (arbitrary-precision arithmetic), GiNaC [2] (symbolic manipulation), yaml-cpp [3] (input parsing), FiniteFlow [4] (finite-field linear solving and rational-function reconstruction), FireFly [5] (reconstruction of symbolic systems), and a C++ thread pool [6]

Nature of problem: Integration-by-parts (IBP) reduction expresses a large number of Feynman integrals of a given family as linear combinations of a small set of master integrals. The Laporta algorithm generates IBP and Lorentz-invariance (LI) identities for all seed integrals below prescribed depth and rank bounds and solves the resulting linear system. For high-loop, many-leg and multi-scale families the system grows as a power of the seed bounds. Two problems follow. First, the elimination has to be repeated at every prime of the coefficient reconstruction, and each pass becomes slow and time-consuming as the system grows. Second, the complete set of raw equations can exceed the memory of a single machine, so the system cannot even be assembled as a whole.

Solution method: The reduction is carried out sector by sector. Each sector receives the targets that higher sectors and the user require, derives an adaptive seed box from them, cuts the upper-right corner of the depth–rank rectangle, can truncate the seeds that are expected to be intermediate integrals, and, when enabled, widens the truncation if the settled basis shows that it cut too deeply. It also filters subsector columns far below the box top. The retained system is solved over a finite field. Only the linearly independent equations and the local masters are kept. The global stage assembles these equations and reconstructs the coefficients with FiniteFlow for numeric systems or FireFly for symbolic ones. Sectors run in parallel under a memory-aware scheduler that keeps the process’s resident set size below a user-defined budget.

Introduction

Precision calculations for the LHC and for rare decays rest on multi-loop Feynman integrals. Once a process is computed beyond leading order, its diagrams are assembled into families of scalar integrals, and every family is reduced to a small set of master integrals with integration-by-parts (IBP) identities [7, 8, 9]. The masters are then evaluated analytically or numerically. Their number is finite [10], yet choosing a basis that keeps the coefficients small is a problem of its own [11, 12]. Typical applications include Higgs production at three loops [13], charged-current Drell–Yan production at N$^{3}$LO [14], two-loop five-point amplitudes [15], and the $2 \to3$ reactions on the Les Houches precision wishlist [16]. The reduction is the common step, and recent reviews describe its algorithms and its cost [17, 18, 19, 20, 21].

The standard Laporta algorithm [9] turns the reduction into the elimination of a finite system, and the size of that system sets the cost. An integral is characterized by its depth, the total power of the propagators in the denominator, and its rank, the total power of the numerator. The algorithm chooses two bounds, $D_{\max}$ and $R_{\max}$, one layer above the largest depth and rank among the integrals to be reduced, and applies the identities to every integral with depth less than $D_{\max}$ and rank less than $R_{\max}$. These integrals are the seeds. Integrals with the same set of propagators in the denominator form a sector. For a family with $L$ loops, $N$ propagators and $E$ independent external momenta, a sector with $k$ active propagators contains $\binom{D_{\max}-1}{k}\binom{R_{\max}+N-k-1}{N-k}$ seeds, polynomial in the two bounds but exponential in $k$. Every seed carries $L(L + E) + E(E - 1)/2$ identities. Each identity connects an integral to its nearest neighbours, so the generated system is sparse, but the number of rows and the total number of terms grow with the same combinatorial factor. The three-loop five-point family of Ref. [22] already lies beyond the reach of the standard programs and required a dedicated syzygy-based generator. At four loops the equation count reaches the ten-billion scale: for the four-loop cut family F1 studied in § 7, a global Laporta run has to generate $1.02 \times10^{10}$ equations before the first elimination step. Holding such a system exceeds the memory of a single node, and eliminating it at this scale takes days or months even after it is partitioned to fit. Memory and run time therefore limit the reduction.

Three quantities set this limit: the size of the equation system that has to be generated, the time its elimination takes, and the memory needed to hold it. The $2 \to3$ frontier is one reason these quantities grow [16], since an extra external leg adds a scale and a loop-momentum routing. Higher powers of propagators in four-loop corrections, additional mass scales and larger master bases have the same effect. Programs have automated the Laporta algorithm [23, 24, 25, 26, 27, 28] and computer-algebra systems make the generation feasible [29], but in demanding calculations the reduction, together with the reconstruction of its rational coefficients, dominates the computing cost [20].

Part of this cost is arithmetic, and it now has a standard treatment. Instead of eliminating over the rational numbers, where intermediate coefficients grow without bound, current programs eliminate modulo a large prime and recover the exact result by Chinese remaindering and rational reconstruction. The dataflow framework of FiniteFlow [4, 30] and the interpolation library FireFly [5, 31] supply this machinery, and Kira 2 [32], Kira 3 [33] and FIRE 6/7 [28, 34] use it by default. This machinery is now common infrastructure rather than a distinguishing feature of any single program. Recent work extends it: balanced reconstruction of sparse rational functions [35, 36] and factorized storage of denominators during the reduction [37].

Another line of work compresses the equation system through algebraic structure before the elimination. Unitarity-cut and syzygy methods construct IBP relations with a controlled propagator degree [38, 39, 40, 41, 42, 43, 44], the group structure exposes relations among the identities [45], NeatIBP generates a small set from module intersections [46, 47], and Blade builds a block-triangular form and reduces spanning sector by spanning sector [48]. More recent constructions order the reduction by a geometric relation and produce $\varepsilon$-factorized differential equations [49], use the reduction operators of the GKZ system associated with a graph [50], or fix the coefficients from analyticity constraints with a few IBP samples [51]. Each of them reduces the size of the system that reaches the elimination independently of the arithmetic used to solve it. Their constructions, however, are expensive in their own right.

A further line of work controls the seed set. Kira 3 introduces a truncation of the seed set, truncate_sp, which keeps only the seeds expected to generate useful equations and leaves the others as intermediate integrals [33]. Blade restricts the seed set through its block and spanning-sector structure [48]. Tube seeding restricts it to a thin path between a target and the masters, so that the number of seeds grows linearly rather than polynomially with the numerator power [52]. Machine-learning approaches search for the priority function that selects the seeds to retain, by genetic programming or language-model search [53, 54], or by reinforcement learning over the reduction steps [55]. Seedless formulations avoid the seed set altogether [56].

SecRed, the program presented in this article, reduces an integral family sector by sector and, within this organization, improves seed selection, equation-size control, elimination efficiency and memory management. The sector-wise organization is also used by Reduze [57, 58], FIRE [27, 28, 34] and NeatIBP [46, 47]. It is adopted for two practical reasons: only the system of the sector currently being processed has to reside in memory, and the sectors form independent tasks whose run times differ by orders of magnitude, which makes them suitable for dynamic scheduling. Sectors run under a memory-aware scheduler that admits tasks according to the measured resident set size and keeps the process within a user-defined budget (§ 4.6). SecRed uses an adaptive seed box with a corner cut (§ 4.2), can truncate the seed set with self-correcting caps (§ 4.3), applies a subsector level filter (§ 4.4), and expands symmetry relations only from the master integrals of a sector (§ 4.5). On the elimination side, it records the elimination plan once and replays it at every prime, and orders the equations by a head score (§ 5).

The benchmarks (§ 7) cover F234, a two-loop five-point non-planar family, and F1, a four-loop four-point cut family whose global system exceeds the memory of a single node. They test these components and compare SecRed with other reducers. On F234 the global system remains compact, the memory footprint modest and the wall-clock time short. On F1, SecRed completes a reduction whose global system cannot be assembled in memory. The reductions have been cross-checked against Kira 3 [33] and LiteRed [25, 26].

SecRed also serves as the IBP engine of AMFlow [59] for the differential-equation method [60, 61, 62, 63, 64], and the numerical results for the target integrals agree with those obtained with other programs. These benchmarks show that SecRed is competitive with existing reducers and extends the range of reductions that can be completed on a single shared-memory node. Section 2 recalls the definitions and identities used throughout. Section 3 describes the architecture and the sector-level algorithms, and Sec. 4 describes the elimination engine. Sec. 5 documents installation and usage, and Sec. 6 reports the benchmarks. Sec. 7 summarizes the results.

Preliminaries: integration-by-parts reduction

We consider a family of Feynman integrals generated by a graph with $L$ loops, $N$ internal lines and $E$ independent external momenta. Every integral is labeled by a vector of integer indices $\mathbf{a}=(a_1,\ldots,a_N)$,

where $l_i$ are the loop momenta, $d$ is the space–time dimension in dimensional regularization, and the propagators $D_k=q_k^2-m_k^2$ are quadratic in the loop momenta, with $q_k$ a linear combination of loop and external momenta. A negative index moves the propagator to the numerator, $1/D_k^{a_k}=D_k^{-a_k}$ with $-a_k>0$, and $a_k=0$ means that the propagator is absent.

Depth, rank, level and sectors

Two integers characterize the complexity of an integral,

the depth $\tau$ (the total power of the propagators in the denominator) and the rank $\rho$ (the total power of the numerator factors). Their sum is the level,

Integrals are classified by the sector

encoded in a bitmask as the sector ID $\chi(S)=\sum_i s_i2^{i-1}$. We identify a sector with the set of its active propagators, so $k\in S$ means $s_k=1$. We write $\omega(S)=\sum_i s_i$ for the number of active propagators (the popcount) and call

the dot number $\delta$ of the integral: the number of powers by which its positive indices exceed one, i.e. the number of dots on the active propagators. For an integral of sector $S=S(\mathbf{a})$ one has $\ell(\mathbf{a})=\omega(S)+\delta(\mathbf{a})+\rho(\mathbf{a})$. A sector $S'$ is a subsector of $S$ if $S'\subset S$ in the bitmask sense. Sector membership depends only on the signs of the indices: numerator factors contribute to the rank but cannot make a propagator active.

The order in which integrals are enumerated is itself an algorithmic choice, and SecRed uses a single key, the balance key, in three places: the enumeration of the seeds (the integrals to which the identities are applied), the ordering of the subsector columns (§ 4.2), and the global integral numbering. For an integral $\mathbf{a}$,

compared lexicographically. The level comes first so that the key is compatible with the level cap of § 4.2 and the subsector level filter of § 4.4. Depth and rank then separate deep-and-narrow from shallow-and-wide integrals, and the largest denominator power and the largest numerator power favor one high power over several moderate ones. The alternative normal order sorts by depth, then rank, with the index vector as the final tie-break, and is used only for comparison (§ 7.2).

Integration-by-parts identities and the Laporta algorithm

Integrals of total derivatives vanish in dimensional regularization [7, 8]. For every pair $(l_i,q_j)$, where $l_i$ runs over the loop momenta and $q_j$ over all loop and external momenta,

Expanding the derivative and decomposing the quadratic factor on the propagators, $q_j^\mu\partial D_k/\partial l_i^\mu=\sum_n\lambda_{kn}D_n+\beta_k$, gives the index-shift form

with $d_{ij}=\partial q_j^\mu/\partial l_i^\mu$, which equals $d$ when $q_j=l_i$ and zero otherwise. There is one identity for each of the $L(L+E)$ pairs, and each identity connects an integral only to its nearest neighbours in index space. In addition, the invariance of the integral under Lorentz transformations gives the Lorentz-invariance (LI) identities [65]. They add further relations among the integrals of the family and can reduce the number of master integrals. Contracting the generator of the Lorentz group with two external momenta $p_r$ and $p_s$ yields, for $1\le r<s\le E$,

The contraction is antisymmetric in $r$ and $s$, so each unordered pair gives one independent identity and there are $E(E-1)/2$ LI identities. The derivatives in Eq. (9) act on the propagators as in Eq. (8), so the LI identities have the same shift structure.

The shift structure determines how the sector changes under an identity. In Eq. (8) the index $k$ is raised and $n$ is lowered, and every term that raises $k$ carries a factor $a_k$. If $k\notin S$ then $a_k\le0$: the terms vanish when $a_k=0$, and when $a_k<0$ the raised index $a_k+1$ is still non-positive. Raising an index therefore never activates a propagator outside $S$. Lowering $n$ either keeps it active ($a_n\ge2$), removes it from the denominator ($a_n=1$, so the integral passes to the subsector $S\setminus\{n\}$), or only increases a numerator power ($a_n\le0$). No operation activates a propagator outside $S$, and at most one active propagator is removed at a time. The LI identities produce the same two shift types, so the same reasoning applies to them. Consequently every integral generated from an integral of $S$ by an IBP or LI identity belongs to $S$ or to a strict subsector of $S$. The equations generated from the seeds of $S$, together with the master integrals they determine, therefore depend only on $S$ and its subsectors.

The IBP and LI identities express every integral of a family as a linear combination of a finite set of master integrals, with coefficients that are rational functions of $d$ and the kinematic invariants [10]. The Laporta algorithm [9] organizes the reduction as follows. The integrals to be reduced are the targets. Two bounds $D_{\max}$ and $R_{\max}$ are set from the targets, one layer above the largest depth and rank among the targets of the sector, and every integral with $\tau<D_{\max}$ and $\rho<R_{\max}$ becomes a seed. Every identity is applied to every seed, producing a linear equation. Neighbours outside the bounds become new unknowns, while neighbours in lower sectors are moved to the right-hand side and reduced separately. The system is ordered by a complexity measure and solved by Gaussian elimination. The integrals that cannot be eliminated are the masters. This is the standard construction. Secs. 3.2 and 3.3 describe how SecRed selects the seeds and can truncate them, and Sec. 3.4 how it filters the subsector columns. Which masters are reported depends on the targets: an irreducible integral on which no target depends is not counted.

Trivial sectors and symmetry relations

The Schwinger parametrization of Eq. (1) reads

where the first and second Symanzik polynomials [66, 67] are sums over spanning trees and spanning two-forests,

with $p_{T_1}$ the momentum flowing into $T_1$ through the cut lines. $\mathcal{U}$ and $\mathcal{F}$ are homogeneous of degree $L$ and $L+1$ in the $\alpha_i$. A sector is trivial if all its integrals are scale-free, i.e. if they carry no dimensionful scale, so that they vanish identically in dimensional regularization. Triviality can be decided from the Symanzik polynomials: writing $G = \mathcal{U} + \mathcal{F}$, a sector is trivial if, after all $\alpha_i$ outside the sector are set to zero, the linear system

has a solution [25, 68]. Trivial sectors can therefore be discarded before any equation is generated.

Sectors related by a relabeling of the propagators belong to the same integral family. The relabeling is detected from the Symanzik polynomial $G = \mathcal{U} + \mathcal{F}$ of Eq. (12): setting to zero all $\alpha_i$ outside a sector gives its restricted polynomial, and two sectors are symmetric if their restricted polynomials become identical under a relabeling of the active propagators. If the relabeling is a pure permutation that leaves the family invariant, the integral is unchanged,

where $\sigma$ denotes the induced permutation of the index vector. A general relabeling is accompanied by a transformation of the loop momenta. The relabeling is promoted to a momentum map by requiring

where $q_k$ is the momentum of propagator $k$, $\sigma$ is the relabeling, and $M$ and $K$ are constant matrices acting on the loop and external momenta. Rewriting the propagators in the transformed basis gives a linear relation among integrals with shifted indices,

where $\mathcal{N}(\mathbf{a}) \subset\mathbb{Z}^{N}$ is a finite set of shift vectors and the coefficients $c_{\mathbf{n}}(\mathbf{a})$ are rational functions of $d$ and the kinematic invariants. For a pure permutation with $\sigma(\mathbf{a}) \ne\mathbf{a}$, $\mathcal{N}(\mathbf{a})$ contains at most the two shifts $\mathbf{n} = 0$ and $\mathbf{n} = \sigma(\mathbf{a}) - \mathbf{a}$, and Eq. (16) reduces to Eq. (14). Equations (14) and (16) are the symmetry relations [57, 69]. A symmetry relation that contains a master $I_M$ can be solved for $I_M$,

where the sum runs over the other integrals in the relation, so $I_M$ is not an independent master. When the relabeling connects different sectors, the relation links integrals across the sector boundary. The relations are added to the equation system before the elimination. § 4.5 describes how SecRed constructs them.

Finite-field methods and reconstruction

Elimination over the rational numbers, or with symbolic polynomials, suffers from intermediate expression swell. Finite-field methods [30, 70] remove this limitation by carrying out the elimination modulo a large prime $p$, repeating it for several primes, and reconstructing the exact rational function from the modular images [71]. Given residues $r_i \equiv a \pmod{p_i}$ for pairwise coprime moduli, the Chinese remainder theorem gives

A rational number $a = u/v$ with $\gcd(u,v) = 1$ and $|u|, |v| < \sqrt{P/2}$ is uniquely recovered from its residue modulo $P$ by the extended Euclidean algorithm [72]. Polynomial dependence is recovered by Newton interpolation in one variable, by the Zippel algorithm [73, 74] in several variables, or by balanced reconstruction of numerator and denominator [31, 35, 36]. FireFly and FiniteFlow provide these algorithms, and SecRed uses them for the reconstruction in the global stage (§ 4). The number of primes is determined by the size of the largest numerator and denominator and ranges from a few dozen to a few hundred.

Architecture and algorithms for sector-wise reduction

SecRed reduces an integral family in two stages. The first stage processes the non-trivial sectors independently and in parallel, determines their local master integrals, and retains only the linearly independent equations. The second stage assembles those equations into one global system, eliminates it, and reconstructs the reduction coefficients. This section describes the architecture and the four sector-level algorithms. The elimination engine is the subject of § 5.

Overall architecture and data flow

Figure 1 shows the data flow. Stage 1 processes the non-trivial sectors in parallel, ordered by decreasing popcount. Each sector is reduced locally, and the subsector integrals that its retained equations involve are forwarded downwards as targets.

Stage 1 does not compute reduction coefficients. It determines which equations are independent, which integrals are masters of a sector, and which subsector integrals a higher sector needs. The arithmetic is therefore finite-field arithmetic only, and the symbolic coefficients are regenerated for the retained equations before they are written to the global container.

As shown in § 3.2, the equations generated from the seeds of a sector $S$ involve only $S$ and its subsectors. Ordering the unknowns by popcount, then by depth and rank, makes the global matrix block triangular, with each block referring only to itself and to blocks below it. This permits the top-down traversal: a sector can be reduced as soon as its parents have finished, and sectors of the same layer that are not in a subsector relation are independent.

The information that travels downwards is the list of upstream targets. When a sector $S_{\mathrm{p}}$ finishes, every subsector integral that appears in its retained equations and is not reducible inside $S_{\mathrm{p}}$ becomes a target of the sector to which it belongs. The receiving sector then treats it like a target supplied by the user: its reduction must be completed there. The shift argument of § 3.2 shows that such an integral has depth at most the largest depth among the targets of $S_{\mathrm{p}}$: reaching a proper subsector lowers one active denominator, which offsets the at most one unit that the raised index can add. In practice it lies one level below the box top of $S_{\mathrm{p}}$. The initial depth bound of the child therefore does not exceed that of the parent. The propagation is exact: every integral that the global stage may need is either reduced in the sector where it was generated or forwarded to a sector that can reduce it. A sector-level reduction is nevertheless not complete in the global sense. The equations of $S$ may connect integrals of several different subsectors: after the integrals of $S$ itself have been eliminated, what remains can be a relation among integrals of two or more subsectors. Such a relation cannot be obtained by reducing either subsector alone, because the equations of a subsector never contain integrals of a sibling. Relations of this type are called magic relations [32]. They become available only in stage 2, when the equations retained from all sectors are assembled and the subsector integrals become common unknowns. The sector-level master set is therefore a superset of the final one, and the global stage removes the redundancy. This division allows the sector stage to remain local and inexpensive without sacrificing completeness.

Within a sector, the algorithms of § 4.2–§ 4.5 form a pipeline. Each stage can only shrink the set of objects passed to the next:

Here $\mathcal{B}$ is the seed box, $\mathcal{B}_{\mathrm{true}}$ the true seeds, $\mathcal{C}$ the columns that enter the elimination, and $\mathcal{E}$ the independent equations that leave it. The seed box is the conservative starting point, determined from the requests of the sector. The optional truncation removes seeds that are expected to be intermediate integrals, without removing their columns. The level filter removes subsector columns far below the box top. The elimination then acts on the remaining system and returns the local masters and the independent equations. Symmetry relations are added after the local masters are known, since they are expanded from those masters only, and any integral they involve in another sector is forwarded like an upstream target. Each arrow in Eq. (19) is a place where the system can be made smaller without changing the final answer.

Each sector is solved once to determine its linear structure. The coefficients are evaluated at a single point modulo a learn prime $p_0$, the first of a fixed list of primes just below $2^{31}$, and the solve returns three results: the local masters, the independent equations, and the dependencies between them from forward selection [33]. All three are properties of the sparsity pattern rather than of the coefficient values. A coefficient that is nonzero generically vanishes at the learn point with probability at most $D/p_0$, where $D$ is the degree of its numerator, and a union bound over the coefficients of the local system bounds the probability that a structural decision differs from the generic one. The decision is checked as well. Forward selection re-solves the selected equations at $p_0$ and requires every target to reduce to masters, repairing a missing pivot or falling back to the full closure when it cannot, and a stored coefficient that vanishes at a later prime sends the replay to a fresh elimination order (§ 5.1). A wrong structural decision therefore triggers a fallback or a rebuild rather than propagating unchecked. The plan that records the decisions is reused at every prime.

The coefficients of the selected equations are regenerated afterwards. Only the equations that survived forward selection are regenerated, not the raw system. Writing $E_{\mathrm{fwd}}(S)$ for the number of equations forwarded by sector $S$ and $R_{\mathrm{gen}}(S)$ for the number generated, the ratio is

The relation is approximate because the global count also contains the symmetry relations, which are expanded from the local masters rather than generated from the seeds. On F234 they are 204 of the 1,162,358 forwarded equations, so the difference is negligible. The ratio is 8.0% on F234 (1,162,358 of 14,478,997) and 3.1% on F1 at $n = 0$ (34,557,830 of $1.11 \times10^9$), with both families defined in § 7. More than 90% of the generated equations never receive a symbolic coefficient.

For a symbolic system the regeneration uses a row table shared by the whole family. For every IBP identity it stores the coefficient of each index shift as a polynomial in the seed indices, so the coefficient of a particular row is obtained by evaluating that polynomial at the row’s index vector rather than by expanding the identity again. Identities whose coefficients are not polynomial in this form fall back to GiNaC and are counted separately. For a numeric system the coefficients are not regenerated at all: the selected rows are replayed at each prime by the elimination plan. In both cases the structural work of stage 1 is done once, and only the retained equations carry coefficient data into stage 2.

Adaptive seed box and corner cutting

For a sector $S$ with bounds $D_{\max}$ and $R_{\max}$, the system has one column per integral with $\tau(I) \le D_{\max}$ and $\rho(I) \le R_{\max}$. An integral generates equations only when

where $\ell_{\max}$ is a cap on the level. The first two conditions define a rectangle in the depth–rank plane, and the third cuts off its upper-right corner. An integral on the far face, with $\tau= D_{\max}$ or $\rho= R_{\max}$, is kept as a column because it can appear as a term in an equation, but it seeds no equation of its own. The starting box is determined by the requests of the sector, $\mathcal{Q}(S)$ for user targets and $\mathcal{P}(S)$ for upstream targets,

and the level cap is placed one layer below the corner,

Two request sets enter Eq. (22): the user targets $\mathcal{Q}(S)$ and the upstream targets $\mathcal{P}(S)$ registered by the parents. A box built from $\mathcal{Q}(S)$ alone is the target box, and the union gives the request box, which is the one used. The largest level admitted by the box is the box top

Initially, Eq. (23) gives $\ell_{\mathrm{top}} = D_{\max} + R_{\max} - 3$. The corner integral, at depth $D_{\max} - 1$ and rank $R_{\max} - 1$, is outside the box.

The corner cut is exact and inexpensive to evaluate. For a sector with $t = \omega(S)$ active and $m = N - t$ inactive propagators, the number of seeds of depth $\tau$ and rank $\rho$ is $\binom{\tau-1}{t-1}\binom{\rho+m-1}{m-1}$. Summing over the rectangle gives

so that the cut removes the fraction

For a representative F234 sector with $t = 5$, $D_{\max} = 9$, $R_{\max} = 4$ and $m = 4$, this gives $f_{\mathrm{cut}} = 20/56 \approx35.7\%$: more than a third of the seeds are not generated, without changing which integrals are reachable.

The initial box is a starting point, not a guarantee. In a round the sector system for the current box is built and solved, and the local masters that the requests depend on are then inspected. Because the elimination orders integrals by the balance key of Eq. (6), whose first criterion is the level, a master needed by a target has no larger level than the target. The corner cut therefore has three cases. If the corner holds no target, no needed master can appear there, so the cut is safe and the first round settles. If the corner holds a target but no master appears there, the target is expressed through integrals below the corner, so the cut is safe and the first round settles. If a master appears at the corner, the corner integral has to generate an equation. Only this third case produces a master above the top,

so the round is discarded and $\ell_{\max}$ is raised by one level. The raised cap admits the corner, and the second round settles. Whichever round settles, its system is the one retained. The discarded first round is cheaper because its box is smaller. Algorithm 1 states the procedure.

Algorithm 1 Adaptive seed box (one sector $S$)

Require: $S$, requests $Q(S) \cup P(S)$

1. $D_{\max} \leftarrow\max_t \tau(t) + 1$; $R_{\max} \leftarrow\max(1, \max_t \rho(t) + 1)$

2. $\ell_{\max} \leftarrow D_{\max} + R_{\max} - 2$; $\ell_{\mathrm{top}} \leftarrow\min(D_{\max} + R_{\max} - 2, \ell_{\max} - 1)$

3. for round $= 1, 2$ do

4. enumerate the columns, select the seeds of Eq. (21)

5. apply truncation (§ 4.3) and level filter (§ 4.4), and solve the sector system

6. if no needed master $m$ has $\ell(m) > \ell_{\mathrm{top}}$ then

7. settle: return local masters and independent equations

8. end if

9. raise $\ell_{\max}$ by one level to admit the corner

10. discard the round

11. end for

The measurements on F234 at $n = 1$ with the truncation off (§ 4.3) show how the corner cut performs. It succeeds at the first round in 208 of the 219 non-trivial sectors, a rate of 95.0%. The remaining 11 need one extension. Across the family, the cut removes 24.5% of the seeds before truncation, from 1,421,510 to 1,073,519.

Adaptive seed truncation

The adaptive seed box of § 4.2 follows a conservative rule: every unreduced subsector integral that appears in a selected independent equation is registered as an upstream target and must be expressed in its own sector. The rule leaves no intermediate integral unreduced, but the sector stage has no global view. An intermediate integral that a sector is required to reduce need not lie on a reduction path that the global stage uses, and the box covers requests that the sector stage may not need.

The practice in Kira 3 [33] is more aggressive. Its seed truncation (truncate_sp) generates equations only from seeds inside depth and rank caps, and keeps the integrals outside them as intermediate integrals that remain columns without generating equations. The caps can be tightened so that only the user targets are guaranteed to reduce. The upstream targets are then ignored for the moment and forwarded to the global stage as irreducible integrals. The expectation is not that the global stage reduces these integrals, but that it reduces the user targets along paths that avoid them. The risk is that a truncation that is too tight removes a reduction path that the global stage cannot reconstruct, and Kira 3 relies on manual tuning of the truncation to reach a minimal master basis.

The seed box fixes which integrals must be present in the system, not which of them must generate equations. The truncation makes that second choice, separating the true seeds, which generate equations, from the intermediate integrals, which remain columns without equations. SecRed implements the truncation as a switch, off by default, and makes the choice adaptive rather than manual when it is on. Within a sector, a refused seed that the round needs is detected from the returned masters, and the cap on the corresponding axis is lifted for the next round. Each sector also records the requests that its truncation refuses. If a final master in the settled global basis is one of these requests, it can be traced back to the sector whose cut was too tight, and the truncation rule is relaxed before the family is reduced again. The caps are computed from the input targets as follows.

Let $Q_{\mathrm{in}}$ be the targets read from the input file, and let $t_{\mathrm{top}}$ be the number of propagators of the top sector, whose active propagators are exactly those that carry a positive power in at least one input target. Define

with $\rho_{\mathrm{exc}}(S)=0$ if no input target belongs to $S$. The level offset and the rank cap of the sector are

where $t_S=\omega(S)$ and $o$ is the offset seed_level_offset. The offset is 0 if all input targets have $D_\delta=0$ and 1 otherwise, unless the input sets it explicitly. The two truncation caps are

and a seed generates an equation only if, in addition to Eq. (21),

A truncated seed is not removed from the system. It remains a column and receives a column weight, so a target can still be expressed through it. It does not generate equations. The truncation is therefore self-correcting. If a truncated intermediate integral proves to be needed, it appears among the masters returned by the round, and the program can detect that the truncation refused a needed seed and lift the cap on the corresponding axis for the next round.

Which refused masters trigger a correction is controlled by probe_mode, which is consulted only while the truncation is on. In targets mode, only masters on the reduction path of an input target are inspected. A master needed only by an upstream target is accepted as a local master and left to the global stage. In all mode, every refused master triggers a correction. Targets mode keeps the truncation adaptive rather than recursive: an upstream target cannot force a sector to enlarge its box, and then its parents to do the same, along the whole request chain. The price is that the local basis may contain a few integrals that a larger box would have eliminated. The global stage can still remove them if the equations suffice.

Two feedback checks act after the global basis has settled. Let

be the largest dot number in the final basis $\mathcal{M}$. If $D_{\delta}^{\mathrm{settled}}>D_\delta$, the truncation has removed a dot layer that the basis needs. $D_\delta$ is then raised to the observed value and the whole run is repeated. If the settled basis contains a request that the truncation refused, the request truncation is lifted as a whole and the run is repeated. Lifting single requests does not converge, because releasing one request changes the targets propagated to its subsectors. Both checks share a limit of four passes, and the final pass is the reported result. Both feedback reruns occur in the samples distributed with the program: tennis raises the dots limit and L04 lifts the request truncation.

On F234, the input targets give $t_{\mathrm{top}} = 8$, $\rho_{\max} = 4$ and $D_{\delta} = 2$, so the automatic offset $o = 1$ gives $\ell_{\mathrm{off}} = 8 - 4 - 1 = 3$. For a five-line subsector that carries no input target, $\rho_{\mathrm{exc}} = 0$ and Eq. (30) gives

so the true seeds of that sector are the box seeds with $\tau\le7$ and $\rho\le3$, while the rank-4 seeds remain columns without generating equations. A sector that does carry a rank-4 input target has $\rho_{\mathrm{exc}} = 4$, and the rank cap is lifted back to $R_{\max}$. This is the exception mechanism in its simplest form. With the truncation switched off, the run admits the whole seed box and generates 14,478,997 raw equations. With the switch on in targets mode, 182 of the 219 non-trivial sectors are truncated, no cap is lifted by the truncation probe, and the run generates 1,113,377 raw equations, 92.3% fewer. Both runs settle on the same 309 masters.

When the truncation is on, it also masks the corner cut. In a truncated sector at least one truncation cap lies below the corresponding box bound, so every seed inside the caps has $\ell\le D_{\max} + R_{\max} - 3$, at or below the box top. The corner cut removes only seeds above that top, so it changes nothing in the generated system of the 182 truncated sectors. Its effect is confined to the 37 sectors without truncation. The cut still contributes: on F234 with the truncation on, it removes 9.6% of the seeds, from 94,781 to 85,681, and acts in 36 of those 37 sectors.

Subsector level filtering

The subsector filter rests on a linear-algebra fact. In Gaussian elimination performed from left to right, the pivot of a column depends only on that column and the columns before it, so deleting columns to the right cannot change which of the earlier columns become masters.

The filter is motivated by the cost of the local elimination. The local system is ordered so that the sector’s own integrals come first and the subsector columns after them. Most of the columns are subsector integrals generated by the identities: a single equation can refer to several different subsectors, and in a sector with many active propagators they can exceed 80% of the columns. Eliminating the complete local system of one sector can be even more demanding than eliminating the whole global system, because the subsector integrals at the back of the equations bring a large number of terms and much complexity into the local elimination, while they do not affect which of the sector’s own integrals at the front become masters. A natural approach is to filter all subsector integrals. This removes the subsector complexity without changing the mathematical structure of the sector’s own reduction (FIRE, for example [34]). It also removes information, however. If all subsector integrals are masked, the magic relations of § 4.1, which appear when the sector’s own integrals cancel among themselves, collapse to zero, and the global stage loses the information it uses to remove redundant masters. The problem is therefore a balance between efficiency and information: the filter should remove enough subsector integrals to make the local elimination faster while losing as few potential magic relations as possible.

Filtering a column means omitting its term from every equation and keeping the rest, i.e. projecting the system onto the retained columns. We work in the balance ordering of § 3.1 and keep only the subsector integrals in the top $n$ levels, counted by their distance to the box top.

A subsector integral $J$ is retained when its level satisfies

and discarded otherwise. The offset $n$ is subsector_filter_n, default 1. At $n = 0$ the threshold is $\ell_{\mathrm{top}} + 1$, so every subsector column is discarded and the system contains only the sector’s own integrals. At $n = 1$ only the top level is kept. As $n$ grows, the retained band widens until the unfiltered Laporta system is recovered. We denote this limit by $n = \infty$, although the value is finite in practice. Figure 2 shows the retained band and the corner cut in the depth–rank plane.

The effect is measured by scanning $n$ on F234 (Table 3). At $n = 0$ the basis has 314 masters. At $n = 1$ it drops to 309 and stays there for $n = 2, 3, \infty$, while the wall-clock time grows from 1m11s at $n = 0$ to 12m28s in the unfiltered limit. The five extra masters at $n = 0$ are the price of discarding all subsector columns: they are reducible integrals that the cross-subsector relations of § 4.1 would have removed. The default $n = 1$ is the smallest offset that closes the reduction on this family. On F1 the trade-off reverses: $n = 0$ gives 566 masters in 1h14m43s, while $n = 1$ gives 559 masters but takes 3h24m30s, so the larger family is run at $n = 0$ (§ 7.3).

$n$	masters	raw eq.	to global	to recon.	wall	peak
0	314	14,221,991	722,748	531,490	1m11s	6.68 GB
1	309	14,478,997	1,162,358	555,833	1m30s	6.97 GB
2	309	14,479,555	1,953,326	531,560	1m40s	7.56 GB
3	309	14,480,135	2,940,516	528,117	2m51s	8.29 GB
$\infty$	309	14,481,498	3,464,762	519,090	12m28s	11.60 GB

Symmetry relations from local masters

Besides the magic relations of § 4.1, symmetry relations form a second class of supplementary relations that remove redundant masters. In the standard construction one relation is generated from every seed, so the number of relations is proportional to the size of the seed set.

A relation is obtained by rewriting each propagator of the seed in terms of the propagators $D'_j$ of the target sector. The map of Eq. (15) gives the loop-momentum transformation $l = M^{-1}(l' - Kp)$ and the correspondence $k \leftrightarrow\sigma(k)$ between the propagators. Substituting the transformation into $D_k = q_k^2 - m_k^2$ expresses the propagator as

where $f(k,j)$ is the coefficient of $D'_j$ and $c_k$ is the constant remainder, which contains the mass shift of the mapped propagator. A seed that carries numerator powers produces a relation whose number of terms grows with the rank of the seed. If propagator $k$ is replaced by $f_k$ terms, where $f_k$ is the number of nonzero coefficients $f(k,j)$ of Eq. (35), a numerator power $|a_k|$ multiplies the number of terms of the relation by $\binom{|a_k|+f_k-1}{f_k-1}$, so

which is smallest for the low-rank objects.

In the sector-wise architecture, we observe that generating symmetry relations only from the local masters suffices to remove the redundancy. This observation greatly simplifies the symmetry relations added to the global reduction. Symmetry relations are generated in the settled round, after the local masters of the sector are known. On F234 this produces 204 relations. Generating from all seeds of the sector instead produces 703,005, while the final basis is the same 309 masters in both cases. The local masters are the lowest-complexity integrals of the sector, so the 204 relations are both fewer and individually cheaper than the 703,005, and the total saving is larger than the ratio of the counts suggests. The sector-wise organization enables this saving, because the local masters of a sector are known before the global reduction.

When a relation involves an integral of another sector, that integral is registered as an upstream target of that sector and the dependency is added to the scheduling graph.

Memory-aware scheduling

Sector tasks are submitted to a thread pool of size sthreads. Two scheduling modes are available. In layer mode the sectors are grouped by popcount and a barrier separates consecutive layers, which guarantees the dependency order and allows the upstream targets registered by a layer to be merged before they are propagated, which removes duplicates. In fine mode, the default, a sector waits only for its direct parents through an atomic counter. The root sectors start immediately, and a finishing sector decrements the counters of its children. The fine mode keeps more threads occupied when sector run times differ widely.

Concurrency is limited by the measured memory of the process, not by a per-sector estimate. Let $M_{\max}$ be the budget (max_memory) and $M_{\mathrm{res}}$ the reserve (memory_reserve). The program reads the resident set size $M_{\mathrm{used}}$ from /proc/self/status and computes

A task is admitted by a three-tier gate. It is admitted unconditionally if no sector is running, and directly if $M_{\mathrm{avail}} > 2M_{\mathrm{res}}$. In the intermediate range it is admitted only if a binary flag $B$ is clear,

The flag records that a task has been admitted in this range, and a single admission sets $B = 1$, so that the free space released by one sector is taken by one sector rather than by all waiting tasks at once. The flag is cleared when $M_{\mathrm{avail}} \le M_{\mathrm{res}}$ or when $M_{\mathrm{avail}}$ has dropped by more than $\Delta= 10$ GB from its value at admission. The first tier excludes deadlock: if nothing runs, one task is always admitted.

The reconstruction stage has its own budget, because each worker holds a complete copy of the reduction store for the whole stage. After measuring the store size $S_{\mathrm{store}}$ on a single-prime trial solve, the pool is sized as

where $M_{\mathrm{mach}}$ is the total physical memory of the machine. The pool is bounded by the number of cores and by one less than the number of primes. On F1 at $n = 1$, Eq. (39) gives 48 workers out of 64, with a measured 7.68 GB per worker inside a 420 GB budget. The run peaks at 425.30 GB.

Optimizations of the elimination engine

The sector stage decides which equations are retained and which integrals are local masters. The global stage then solves the assembled linear system many times, once for each pair of a prime and a sample point used in the reconstruction. Generically all these solves have the same rows and the same positions of nonzero coefficients, and only the coefficient values change. Sec. 4.1 uses this repetition: the order of the elimination is computed once and stored. Sec. 4.2 uses the freedom that remains inside that order: rows with the same first unknown can be processed in any order, and the choice changes how many terms the elimination moves.

Elimination plan and replay

Let the global system be $Ax = 0$ with $m$ rows. The unknowns are the integrals that occur in the retained equations, numbered as in Sec. 2.1. A row is a list of terms, each term a pair of an unknown number and a coefficient, and the row length is the number of terms. The coefficients are rational functions of $d$ and the kinematic invariants. Because they are reconstructed from values in finite fields, the same system is solved at every pair of a prime and a sample point.

An elimination processes the rows in some order and solves one unknown per row. The terms of a row are ordered by decreasing unknown number, and the first unknown in this order is called the head of the row. When a row is reached, its first unknown that is not yet solved becomes its pivot. The row is multiplied by the inverse of the pivot coefficient, so that the pivot has coefficient one, and the result is stored as the solution of the pivot. A later row that contains an already solved unknown is reduced by adding a multiple of the stored solution, which cancels that unknown. This addition is a merge. The order of the rows depends only on the positions of the nonzero coefficients, not on their values.

The elimination plan is the stored record of the structural part. For each row $r$ it stores

the head $h(r)$, the row length $n(r)$, and the equation number $e(r)$, which is the position of the row in the generated system and makes the order deterministic when all other entries are equal. The plan also stores the order in which the rows are processed. It does not copy the terms: all rows are kept in one contiguous array, the flat store. Each row occupies a segment of the array, and its terms are stored as consecutive pairs (unknown number, coefficient) in the order stated above. The plan refers to a row by its location in the array. Keeping all rows in one array rather than in one object per row makes the memory access sequential.

The order is obtained by sorting the rows with a key that groups them by head, so that rows with the same head are adjacent. The full key is given in § 5.2. Given the order, the pivot of every row is fixed by the rows before it, and so is the sequence of merges. The solved unknowns in a row are the pivots of earlier rows that occur in it, and their order follows from the stored order and the stored list. Replay therefore repeats neither the search for the pivots nor the sorting of the rows at a new sample point. The plan never reads a coefficient, so replay reproduces the direct elimination.

To replay the plan, the worker writes the coefficient values of the current point into the store and follows the stored order. For each row it takes the pivot fixed by the order, merges into the row the stored solutions of its solved unknowns, normalizes the row, and stores the result. After the forward pass, the solutions are substituted back in reverse order, which is the usual back substitution. Two details reduce the arithmetic. The inverse of the pivot coefficient is computed once and applied to every coefficient of the row, so the division in the normalization becomes a multiplication,

where $a_\pi$ is the pivot coefficient and $\operatorname{inv}_p$ is its inverse modulo $p$ (Shoup precomputation [75]). A merge of two sorted rows is linear in their lengths, and its result is again sorted.

The plan remains valid while the number of rows is unchanged and no coefficient that the plan uses vanishes at the current point. The check is conservative: a zero anywhere in the store, not only at a pivot, invalidates the order and triggers a fresh computation from the nonzero entries. For random coefficients the expected number of zeros is $N_t/p$, where $N_t$ is the number of stored terms. At $p \simeq2^{31}$ this is 0.18 for the F1 system, so the plan is reused for most primes and rebuilt for the rest. A rebuild costs one sort of the row list and does not change the result.

The same plan and the same store serve every solve of the global stage: the solves of the reconstruction at its primes and sample points, the smaller re-solve that checks the selected equations against the targets, and the solve that the symbolic backend performs before it reconstructs.

The plan was measured on F234 with symbolic $d$ in six runs that combine three plan settings with the two states of the seed truncation. The settings are the full plan, the plan without the recorded merges, which keeps the row order and finds the pivots at each point, and no plan at all, which orders the rows and finds the pivots at each point. The six runs return the same 309 masters and byte-identical master and reduction files, so the plan changes the cost and not the result. With the truncation off, the global stage takes 1m14s with the full plan, 1m44s with the row order alone, and 2m1s with no plan. With the truncation on, the times are 1m55s, 2m26s and 2m35s. Most of the gain over the per-point elimination comes from recording the merges, not only the order of the rows.

Head score ordering

The number of equations selected for reconstruction directly affects memory, but it is not simply positively correlated with the reconstruction time. In the finite-field reconstruction of one family, fewer selected equations do not always mean a faster reconstruction in practice, and a preliminary conclusion is that the relevant quantity is the total number of terms moved during the elimination at each sample point. We introduce a head score into the row ordering to reduce that number.

The head fixes the first unknown of a row but not the order of rows that share it. Several rows can share the same head, and they form a head group. A common choice orders the rows of a group by length, shortest first. The head score below takes precedence over the length inside the group, which remains as a tie-breaker. A row contains only unknowns whose number does not exceed that of its head. Consider an unknown of the row that is the head of another row. If it is smaller than the head of the current row, its row has a smaller head and is processed earlier, so the unknown is solved by the time the current row is reached and its solution is merged into the row. The head of the current row is solved before the current row is processed only when another row of the same group precedes it, and it is the pivot of the first row of the group. The unknowns that enter a row through merges are the fill-in of the elimination. The more head unknowns a row contains, the more merges it receives before it reaches its pivot, and the more terms the elimination moves.

Let $H = \{h(r) : r\ \text{nonempty}\}$ be the set of heads and let $S_r$ be the set of unknowns in row $r$. The head score of the row is

When the row is processed, every head unknown smaller than its own head has already been solved and contributes one merge, and the head itself contributes one unless the row is the first of its group. Unknowns outside $H$ may also be solved before the row is reached, but they are not counted. A merge can introduce further head unknowns through the stored solution, so $\sigma(r)$ is an approximation rather than an exact count. It measures the part of the row growth that can be predicted from the lists of unknowns before the elimination starts, and a small score identifies a row that reaches its pivot after few additions.

The rows are sorted by the tuple

compared lexicographically, where empty rows are placed last, $h(r)$ brings the rows of a head group together, $\sigma(r)$ orders them inside the group in increasing order, and the remaining three entries break ties. Here $\mathcal{L}(r)$ is the list of the unknowns of row $r$ in the stored order. The increasing score places the row with the fewest head unknowns first in its group. That row becomes the stored solution of the shared head, so a small first row keeps the solution that the other rows merge small. The score is computed from the lists of unknowns and the heads alone, so it does not depend on the coefficient values: it is used when the plan is built and recomputed by the same rule when the plan is invalidated. It changes only the order in which the rows of a head group are processed, not which integrals are reducible.

On F234 with the truncation on, dropping $\sigma(r)$ from the ordering key raises the selected set from 261,145 to 270,633 equations, the terms moved per sample point from $5.35 \times10^{8}$ to $6.68 \times10^{8}$, and the reconstruction time from 1m28s to 1m53s. With the truncation off the score is neutral: it selects 555,833 equations against 551,148 and moves $2.24 \times10^{8}$ terms against $2.20 \times10^{8}$, and the two runs differ by about 1 s. The ordering is greedy: each row is placed by its own score, without a search over the possible orders, so it is not optimal. It is chosen because it is computable from the positions of the nonzero coefficients alone and therefore adds no cost at replay time. Because the elimination efficiency directly determines the reconstruction efficiency, this direction merits further study.

Installation and usage

Building the program

SecRed is written in C++ and builds with CMake. It requires GMP and MPFR, FLINT for arbitrary-precision rational arithmetic, GiNaC for symbolic manipulation, yaml-cpp for input parsing, jemalloc as the allocator, and two external libraries that are linked at build time: FiniteFlow for modular linear algebra and rational reconstruction, and FireFly for reconstruction of symbolic systems. On Debian or Ubuntu the packaged dependencies are installed with

sudo apt install build-essential cmake ninja-build libgmp-dev \
libmpfr-dev libflint-dev libyaml-cpp-dev libginac-dev \
libcln-dev libjemalloc-dev

FiniteFlow and FireFly are built from source and installed with their own configure scripts. The SecRed source is cloned from the repository listed in the Program Summary with

git clone https://github.com/chiyutuci/secred.git
cd secred

and the program is built with

cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build -j

which produces build/secred and the test binary build/secred_tests. A dependency installed under a non-standard prefix is passed with -DCMAKE_PREFIX_PATH. The build links FiniteFlow and FireFly into the binary, so neither has to be on PATH at run time. The test suite is run with ./build/secred_tests. It covers the scalar-product construction, the IBP and LI generators, the seed box and truncation, the subsector filter, the finite-field solver, and the symmetry maps.

Input file and options

A run is started as secred input.yaml. The input has three blocks. The family block gives the topology: its name, the space–time dimension, the loop and external momenta, the kinematic invariants with their mass dimensions, the scalar-product rules, the propagators as [momentum, mass], and optional cut flags. The targets block lists the integrals to reduce, one index vector per line. The optional preferred block names integrals to be kept as masters when the mathematics allows it. It expresses a preference rather than a guarantee. A minimal input for F234, with four of its 860 benchmark targets, is

family:
  name: F234
  dimension: d
  internals: [l1, l2]
  externals: [p1, p2, p3]
  sps_rules:
    - [p1, p1, 0]
    - [p2, p2, 0]
    - [p3, p3, 0]
    - [p1, p2, 1/3]
    - [p1, p3, 1/7]
    - [p2, p3, 9/2-1/3-1/7]
  propagators:
    - [l1, 3/2]
    - [l1-p2, 3/2]
    - [l1+p1/2-p2/2+p3/2, 0]
    - [l2, 3/2]
    - [l2+p3, 3/2]
    - [l2-p1/2-p2/2+p3/2, 0]
    - [l1+l2-p1/2-p2/2+p3/2, 3/2]
    - [l1+l2+p1/2-p2/2+p3/2, 3/2]
    - [l1+p1, 0]
  settings:
    sthreads: 64        # sector threads
    gthreads: 64        # reconstruction threads

parallel_mode: fine  # fine | layer
ordering_mode: balance  # balance | normal
subsector_filter_n: 1 # filter offset
truncate: true        # seed truncation (probe_mode: targets by default)
symmetry: true        # symmetry relations
max_memory: 512       # GB
memory_reserve: 30    # GB
cache_dir: cache_F234 # output directory
targets:
- [1, 1, 1, 1, 1, 0, 1, 1, -4]
- [1, -2, 1, 1, 1, 1, 1, 1, -2]
- [1, 3, 1, 1, 0, -2, 0, 1, 0]
- [1, 1, 2, 0, 0, 1, 1, 1, -3]

The third target carries a double dot, the fourth a dot, and the first two carry rank 4. The settings block is optional. The options below are grouped by what they control, and each entry gives its default.

Parallelism and memory.

sthreads (default 64). Number of worker threads of the sector stage.

gthreads (default sthreads). Number of concurrent reconstructions in the global stage, i.e. the number of FireFly sample points or FiniteFlow rational reconstructions.

parallel_mode (default fine). In fine mode a sector is scheduled as soon as its direct parents finish. In layer mode the sectors of one popcount run as a batch: their upstream targets are merged and deduplicated before propagation. Fine mode keeps more threads occupied when sector run times differ widely.

max_memory (default 512 GB). Memory budget of the process.

memory_reserve (default 30 GB). Headroom kept below max_memory for allocations outside the scheduler.

max_primes (default 500). Cap on the number of primes in numeric reconstruction. The CRT stops at convergence, so the cap limits a run rather than determining its length.

cache_dir (default cache). Output directory.

Sector algorithms.

subsector_filter_n (default 1). Offset $n$ of Eq. (34).

truncate (default false). The seed truncation of § 4.3. When it is off, the whole seed box is admitted, no truncation rule is derived and seed_level_offset is void. When it is on, the caps are derived from the targets and the truncation probe defaults to targets. Naming a probe_mode also turns the truncation on, with that scope. An explicit truncate: false overrides this setting: it forces the unscoped all walk.

probe_mode (default targets while the truncation is on). Selects which refused masters may widen the box: targets restricts the check to input targets, all lets every out-of-box master widen the box. It is forced to all under truncate: false.

seed_level_offset (default automatic). Overrides the offset $o$ of Eq. (29). By default it is 0 when all targets are pure denominators and 1 otherwise. It is ignored when the truncation is off or probe_mode: all is in effect.

symmetry (default true). Enables the master-based symmetry relations.

Ordering.

ordering_mode (default balance). In balance mode the seeds, the subsector columns and the global integral numbering all use the balance key: total power first, then depth, rank, the largest denominator power and the largest numerator power. The same comparator is used at the three places, so the boundary, the filter and the elimination agree on which integral is simpler. The normal mode replaces this key by depth, then rank, with the index vector as the final tie-break, and is intended for comparison.

Output

All output is written below cache_dir. reductions contains the result: for each target, one line per term in the form $[a_1,\ldots,a_N]$*coefficient, with the target’s own row first (coefficient $-1$) followed by the masters, and a blank line between targets. masters_sector lists the master candidates left by the sector stage, but this is not yet the final basis. masters_global lists the final basis and is written as soon as the global learn pass has fixed it, before the reconstruction, which can be the longest phase of a run. The directory sectors/ holds the per-sector logs and target lists, secred.log is a copy of the standard output, and firefly.log is written on runs that used FireFly. Lines from the FiniteFlow backend carry the [fflow] tag.

Benchmarks

Setup

All runs were performed on a node with two AMD EPYC 7542 processors (128 logical threads), 512 GB of memory and Ubuntu 24.04, with the program built by GCC 13.3 in release mode. The sector and global stages used 64 threads each in fine mode, $n = 1$ and balance ordering unless stated otherwise. The runs use the default truncate: false unless stated otherwise, and the run with the truncation on is discussed in § 7.2. The reported equation counts are the raw equations generated by all settled sector systems, the equations forwarded to global, and the equations to reconstruction selected from the global system. The wall clock is the full run, and the peak memory is the resident set size reported by the program.

The two families are shown in Figure 3. F234 is the two-loop five-point non-planar family of $J/\psi\to\gamma\gamma\gamma$ with $L = 2$, $N = 9$, three independent external momenta and a massive leg $p_4^2 = 4m^2$. It is a medium-size family, and current NNLO calculations routinely require hundreds or thousands of such families per process. The 860 targets have $\max\tau= 8$ and $\max\rho= 4$, giving $D_{\max} = 9$ and $R_{\max} = 5$ over 219 non-trivial sectors. F1 is the four-loop four-point cut family of the same process with $L = 4$, $N = 14$, $E = 1$, $p^2 = 4m^2$ and the first three propagators cut. Its 15,602 targets have $\max\tau= 13$ and $\max\rho= 3$, giving $D_{\max} = 14$ and $R_{\max} = 4$ over 860 non-trivial sectors. These bounds give $5.11 \times10^8$ seeds without the corner cut. F234 is small enough for the Laporta programs in the comparison to finish, and F1 is large enough that a global Laporta run exhausts the node during equation generation.

Table 1 splits the wall clock into the phases reported by the program. The sector stage is a minority of the total in every run, and the global stage dominates. This is why § 5 is devoted to the global elimination: the sector algorithms decide how large the global system is, but most of the time is spent in the global stage.

phase	F234 ($n = 1$)	F234 (targets)	F1 ($n = 0$)	F1 ($n = 1$)
preprocessing	< 1 s	< 1 s	1m17s	6m52s
sector stage	10.1 s	2.2 s	10m16s	1h24m58s
global stage	1m18s	1m55s	58m20s	1h52m35s
total	1m30s	1m59s	1h14m43s	3h24m30s
peak memory	6.97 GB	8.42 GB	386.92 GB	425.30 GB

The phase timers are independent and do not add exactly to the total.

For F234 the propagators are, with $p_4=p_1+p_2+p_3$,

with $p_1^2=p_2^2=p_3^2=0$, $p_4^2=4m^2$, $p_1\cdot p_2=s$, $p_1\cdot p_3=t$ and $p_2\cdot p_3=2m^2-s-t$. The masses on $D_1,D_2,D_4,D_5,D_7,D_8$ are equal. The numeric cross-check uses $d=1999/500$, $s=1/3$, $t=1/7$ and $m=3/2$.

For F1 the propagators are, with $p^2=4m^2$,

with masses on $D_4,D_5,D_7,D_9,D_{10},D_{11},D_{12},D_{13}$, the first three propagators cut, $m=1$ and $d=1999/500$. The input files for both families are in the repository.

Two-loop five-point non-planar family

Table 2 compares the programs on F234 with symbolic $d$. The first SecRed row uses the adaptive seed box with its corner cut and no truncation, and the second also enables the seed truncation, with probe_mode: targets. SecRed and default Kira return 309 masters each, but the bases differ: 250 integrals are common and each program has 59 of its own, because the two programs use a different integral ordering. The truncation-off run forwards 78% fewer equations to the global stage, uses half the peak memory and finishes in 90 s against 501 s. The number of equations that reach reconstruction is essentially the same, so the gain is not in the last stage: it is in the size of the global system and in the memory needed to hold it.

program	raw	global	recon.	masters	wall	peak
Kira 3.1	19,501,427	5,192,771	565,749	309	501.3 s	13.8 GB
Kira 3.1 (trunc.)	1,122,048	464,034	172,561	346	242.4 s	6.22 GB
SecRed (trunc. off)	14,478,997	1,162,358	555,833	309	90.1 s	6.97 GB
SecRed (trunc. on)	1,113,377	414,834	261,145	309	119.0 s	8.42 GB
FIRE 7.1 (sym.)	3,177,403*	2,525,345*	—	318	487 s	1.0 GB
Blade	—	—	—	—	3600 s†	33 GB
NeatIBP	—	—	—	—	3600 s†	25.6 GB

Kira’s per-sector truncation mode, the counterpart of the seed truncation of § 4.3, cuts its own raw system to 1,122,048 equations and its global system to 464,034, both below the 1,162,358 forwarded by the truncation-off SecRed run. The smaller system comes at a price. Kira’s basis grows to 346 masters, a superset with 37 extra integrals in proper subsectors. The truncation shortens Kira’s run from 501 s to 242 s and lowers its peak from 13.8 GB to 6.22 GB, but it remains slower than the truncation-off SecRed run. The per-sector truncation parameters follow the same strategy as the Kira 3 interface of AMFlow 2.0 [76].

The targets row of Table 2 restricts the correction to masters on the reduction paths of the input targets. It returns the same 309 masters and a byte-identical reductions file as the truncation-off run, with 1,113,377 raw equations instead of 14,478,997. It is not cheaper than that run: it takes 119 s against 90 s. The reconstruction receives fewer equations and yet takes 1m28s, against 50 s in the truncation-off run, so the smaller system does not translate into a shorter reconstruction here. This is consistent with § 5.2, where the reconstruction time followed the number of terms moved during the elimination more closely than the number of equations. The truncation-on run moves $5.35 \times10^{8}$ terms per sample point, against $2.24 \times10^{8}$ with the truncation off. When the equation counts differ by less than one or two orders of magnitude, fewer equations can correspond to more complicated reduction paths, so the merges fill in more and move more terms. At $n = 0$ the targets mode returns 545 masters: the narrower box accepts 236 extra local masters, the price anticipated in § 4.3, and they survive into the final basis. The level filter at $n = 1$ supplies the missing subsector equations and restores the 309-master basis. The truncation-off runs are the conservative benchmark used for the scans below. The two programs therefore use truncation to different ends: Kira accepts a larger basis for a shorter run, while SecRed keeps the 309-master basis and a smaller system but does not gain time on this family.

To separate the cost of the equation system from the cost of the coefficients, the same truncation-on run was repeated with $s$ and $t$ symbolic and $d$ fixed at 1999/500. The equation system, the 309 masters and the masters_global file are unchanged. The reconstruction is not: the two-parameter run takes 2h23m13s against minutes for the symbolic-$d$ counterpart, so the number of symbolic parameters sets the interpolation cost. Running the same two-parameter input with truncate: false gives a byte-identical reductions file and a global system larger by a factor 2.8, yet the run finishes in 1h38m25s. The same inversion persists: the truncation-on run moves $4.22 \times10^{8}$ terms per sample point in its merges and spends 1.843 s per sample point in forward elimination, the truncation-off run $1.89 \times10^{8}$ terms and 1.251 s.

FIRE 7.1 with its internal Lee symmetries finishes the same target set in 487 s with 318 masters, nine more than SecRed and default Kira, and its process group peaks at 1.0 GB resident. The nine extra masters come from the absence of a global stage: FIRE reduces each sector on its own, so an integral whose elimination needs equations from several sectors stays in the basis.

Blade and NeatIBP were each given a one-hour cap on the same input and did not produce coefficients. Blade completed 10 of its 11 maximal-cut systems before the cap. NeatIBP reported 100 sectors finished and wrote master and IBP files for 153 sectors, but had not started the reduction when it was stopped. These runs measure the cost of their preprocessing at this scale and do not indicate a failure of the methods. Both are designed for families in which their algebraic or numerical preprocessing is advantageous.

Table 3 and Figure 4 scan the filter offset $n$. The master count is 314 at $n = 0$ and 309 for every $n \ge1$. The five extra integrals at $n = 0$ are reducible through cross-subsector relations that the discarded columns would have supplied. Raising $n$ beyond 1 does not change the basis, and the raw equations stay nearly flat, but the global system grows by a factor 3 from $n = 1$ to the unfiltered limit and the wall clock from 1m30s to 12m28s. The extra subsector columns add few independent sector equations. They add equations that only the global stage can use, and the global elimination bears the cost. The default is $n = 1$. The same scan on F1 (§ 7.3) has the opposite outcome, and is discussed there.

The seed ordering is a separate choice from the filter offset. Table 4(a) compares the balance key of § 3.1 with the normal order on F234 at $n = 1$. Both return 309 masters, but the balance order generates 44% fewer global equations and finishes in 1m30s against 3m43s, with 7.0 GB against 12.1 GB. The balance order also lets the truncation probe succeed in 208 of the 219 sectors at its first attempt, against 157 under the normal order, so fewer sectors incur the cost of a discarded round. The rest comes from the elimination itself, which sees integrals in the same order in which the box and the filter ranked them. The current version deliberately does not offer a user-defined ordering of integral complexity, because the balance order is tightly coupled to the sector-level algorithms: the level cap, the subsector filter and the elimination all depend on it.

Table 4: Configuration studies behind the default settings. (a) F234 at $n = 1$ with symbolic $d$ under the two seed orderings. (b) F1 with numerical $d = 1999/500$ and 15,602 targets. The $n = 1$ run set the memory budget to 480 GB with a 60 GB reserve, and the $n = 0$ run used the same reserve. The Kira row uses per-sector truncation set in its input. Its peak is an observed value during equation generation, not a kernel high-water mark. “Raw”, “global” and “recon.” are as in Table 2.

	masters	raw eq.	to global	to recon.	wall	peak
(a) F234, $n = 1$: seed ordering
balance	309	14,478,997	1,162,358	555,833	1m30s	6.97 GB
normal	309	19,053,931	2,083,375	927,520	3m43s	12.10 GB
(b) F1: filter offset and Kira
$n = 0$	566	1,105,719,427	34,557,830	29,590,539	1h14m43s	386.92 GB
$n = 1$	559	1,591,571,476	156,627,968	39,571,032	3h24m30s	425.30 GB
Kira (trunc.)	703	108,139,890	39,118,736	5,865,036	1h40m	37 GB

Four-loop four-point cut family

Table 4(b) reports F1 at $n = 0$ and $n = 1$, together with the Kira per-sector truncation run. A standalone Kira run without truncation exhausts the node’s memory during equation generation, so no comparison without truncation is available. FIRE 7 does not accept cut propagators, so it cannot be compared on this family either. With per-sector truncations set in its input, Kira completes the family in about 1h40m with an observed peak of about 37 GB, which is not a kernel high-water mark. It returns 703 masters, against 566 at $n = 0$ and 559 at $n = 1$, and the basis contains integrals in which a cut propagator carries power 2 or 3. A cut family requires power exactly one on the cut lines, so those integrals cannot be used in the cut calculation. SecRed at $n = 0$ is faster than this run, 1h14m43s against about 1h40m, but it needs 387 GB against 37 GB and returns a basis that satisfies the cut condition. Obtaining a smaller, cut-admissible basis from Kira requires adjusting the truncation sector by sector, which is impractical in automated batch calculations. SecRed computes the caps from the targets and corrects them when the settled basis shows that they are too tight. At $n = 1$ the basis has 559 masters, 7 fewer than at $n = 0$, and the $n = 1$ basis is a strict subset of the $n = 0$ basis. The cost is substantial: the raw equations grow by 44%, the global system by a factor 4.5, and the wall clock by a factor 2.7. For this family the extra completeness does not justify this cost, and the production setting is $n = 0$.

The F1 runs also exercise the memory-aware scheduler. Most of the memory is used in the reconstruction stage. At $n = 1$ the reconstruction pool is sized by Eq. (39) to 48 of 64 workers, with a measured 7.68 GB per worker inside a 420 GB budget. The run peaks at 425.30 GB, 5.3 GB above the 420 GB pool budget, and the 60 GB reserve absorbs the overshoot. Without the pool, the first attempt at this configuration ran out of memory while filling the reconstruction workers. Sizing the pool from the measured store allows the run to fit within the budget. The thread count therefore needs no manual tuning: the user sets only the maximum thread count and the memory budget, and the scheduler sizes the pool to the largest number of workers that fits the budget.

Summary

We present SecRed 1.0, a program that reduces Feynman integrals sector by sector. Each sector is processed independently and only its linearly independent equations are forwarded to a global stage, so that the working set of a sector remains local and the sectors form independent tasks for dynamic scheduling under a memory budget. The sector stage combines an adaptive seed box with a corner cut, a seed truncation with self-correcting caps, a subsector level filter and symmetry relations expanded only from the local masters. The box and the truncation widen or relax themselves when the settled basis shows that they removed too much. The result is a substantially smaller system entering the elimination, with the self-correction preserving the completeness of the basis and the symmetry relations removing most of the expansions. The global stage records the elimination plan once and replays it at every prime and sample point, so that only the arithmetic is repeated, while the head score orders the rows of a head group so that the merges move fewer terms. These design choices and optimizations, combined with the memory control, allow the program to run efficiently within a limited memory budget.

The correctness of the reductions has been cross-checked against independent implementations. SecRed serves as the IBP engine of AMFlow [59] in the NNLO corrections to $J/\psi\to\gamma\gamma\gamma$ [77], in double $J/\psi$ production at $B$ factories [78], and in an ongoing $J/\psi\to ggg$ calculation. The program is expected to support deeper and more demanding calculations in high-energy physics in the years ahead.

Several directions remain for improvement. The seed set remains the main source of avoidable cost. The adaptive box is conservative and the truncation is aggressive. Although both adjust themselves when a needed integral is missed, a scheme that is smaller while remaining safe, in the sense that it retains the information required for the minimal master basis, has yet to be found. Recent proposals point in this direction: a tube of seeds along the path between a target and its masters [52], and formulations that avoid the seed set entirely [56].

The size of the system, however, is not the only factor in the run time. A smaller system is beneficial, but it does not necessarily accelerate the elimination. The run time depends more strongly on the reduction paths and on the fill-in of the merges than on the number of equations, and the mechanism behind this relation is not yet sufficiently understood to guide the choice of the equation set. Machine-learning methods explore the structure of the equation set by searching for heuristics that decide which equations enter the reduction [53, 54]. The elimination algorithm itself admits further optimization, and the finite-field reconstruction stage also offers scope for improvement, for example through GPU acceleration.

SecRed is designed for large-scale, batch automated reduction, such as the numerical evaluation of master integrals with AMFlow, where many families are reduced in a single campaign. The aim is to achieve efficiency without manual tuning and stable behavior across the batch, and the benchmarks and the tests are selected accordingly. For particular mathematical structures or special requirements, for calculations that permit manual parameter tuning, or for reconstructions with many symbolic parameters, further optimization remains possible.

Author contributions

Ruichang Niu designed the algorithms, implemented and optimized the program, and wrote the manuscript. Chao Zeng tested the program and provided the benchmark cases. Bin Gong and Jianxiong Wang contributed to the discussion of the algorithm design and coordinated the work. All authors reviewed the manuscript.

Acknowledgements

The authors thank X. Chang for numerous discussions on the design and development of the program, and X. Huang for discussions on the benchmark cases. This work was supported by the National Natural Science Foundation of China with Grants No. 12135013.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

The SecRed source code and the input files used for the benchmarks are available in the repository listed in the Program Summary. The program is also archived in the CPC Program Library.