#  FiniteFieldSolve: Exactly Solving Large Linear Systems in High-Energy
Theory

James Mangan  0000-0002-9713-7446 Department of Physics and Astronomy,
Northwestern University, Evanston, Illinois 60208, USA

###### Abstract

Large linear systems play an important role in high-energy theory, appearing
in amplitude bootstraps and during integral reduction. This paper introduces
FiniteFieldSolve, a general-purpose toolkit for exactly solving large linear
systems over the rationals. The solver interfaces directly with Mathematica,
is straightforward to install, and seamlessly replaces Mathematica’s native
solvers. In testing, FiniteFieldSolve is approximately two orders of magnitude
faster than Mathematica and uses an order of magnitude less memory. The
package also compares favorably against other public solvers in
FiniteFieldSolve’s intended use cases. As the name of the package suggests,
solutions are obtained via well-known finite field methods. These methods
suffer from introducing an inordinate number of modulo (or integer division)
operations with respect to different primes. By automatically recompiling
itself for each prime, FiniteFieldSolve converts the division operations into
much faster combinations of instructions, dramatically improving performance.
The technique of compiling the prime can be applied to any finite field
solver, where the time savings will be solver dependent. The operation of the
package is illustrated through a detailed example of an amplitude bootstrap.

###### keywords:

Large linear systems; exact solutions; finite fields; computer algebra
software

††journal: Computer Physics Communications

PROGRAM SUMMARY/NEW VERSION PROGRAM SUMMARY

Program Title: FiniteFieldSolve
CPC Library link to program files: (to be added by Technical Editor)
Developer’s repository link: https://github.com/jfmangan/FiniteFieldSolve
Code Ocean capsule: (to be added by Technical Editor)
Licensing provisions: GPLv3
Programming language: Mathematica, C++
Nature of problem: Exactly solving large linear systems over the rationals
occurs in various settings in high-energy theory, for example when performing
integral reduction or bootstrapping an amplitude.
Solution method: The linear system is solved by repeatedly row reducing over
different finite fields (see Ref [1] and references therein). Finite fields
avoid the intermediary expression swell inherent to arbitrary precision
rationals and bypass roundoff errors from floating point numbers. A downside
to using modular arithmetic is that it introduces a tremendous number of
integer divisions, but this can be mitigated by compiling the divisions down
to simpler instructions. The solver is designed to handle arbitrarily dense
systems such as those that appear in certain amplitudes bootstraps.

## References

  * [1] M. Kauers, “Fast solvers for dense linear systems,” Nucl. Phys. B Proc. Suppl. 183, 245-250 (2008) doi:10.1016/j.nuclphysbps.2008.09.111

##  1 Introduction

Exactly solving large linear systems is a key tool in the theorist’s arsenal.
The most time consuming and computationally expensive steps in many quantum
field theory calculations are integral reduction and integration, where the
former can be reduced to a linear algebra problem [1]. Large linear systems
also appear prominently in amplitude bootstraps, including the double copy [2,
3], the soft bootstrap [4], as well as others. The double copy in particular
has been crucial in understanding the ultraviolet behavior of supergravity [5,
6, 7, 8, 9, 10, 11, 12] and in making precision gravitational wave predictions
[13, 14, 15, 16]. There are a number of publicly distributed linear solvers
[17, 18, 19, 20, 21, 22], including some integrated into integral reduction
software [23, 24, 25]. However, interfacing with these solvers can be a non-
trivial task and some of them are intended purely for very sparse systems.
This paper presents FiniteFieldSolve, a general-purpose exact large linear
system solver for arbitrary density equations over the rationals that can be
called directly from Mathematica. The solver presented in this paper is open
source, a drop in replacement for Mathematica’s Solve and Reduce, easy to
install and use, and roughly two orders of magnitude faster than Mathematica
and uses an order of magnitude less memory in applications tested (see Tables
1(c) and 2(c)).111The author is only aware of two other solvers that link
directly into Mathematica and are intended for exact large linear systems over
the rationals. The first is LinearSystemSolver [22] based on [26] and the
second is spaSMLink [21] which is based on SpaSM [20]. Some of the frontend of
FiniteFieldSolve was adapted from spaSMLink but FiniteFieldSolve is unrelated
to SpaSM. Although FiniteFieldSolve is designed for arbitrarily dense systems,
some of the methods presented here are applicable to other solvers, including
the technique of compiling the prime as described below. FiniteFieldSolve also
has the tremendous advantage that its memory footprint is very stable and
predictable and the time to row reduce a matrix over a given finite field can
be estimated very accurately. This makes FiniteFieldSolve highly predictable
while providing users access to a wide range of physics problems that would
otherwise be intractable, all with the convenience of Mathematica’s high-level
syntax.

As a general-purpose solver, FiniteFieldSolve is equipped to handle fully
dense systems. While sparse solvers will outperform FiniteFieldSolve for very
sparse systems, the ability to manipulate high density equations allows the
user to offload certain challenges to the computer. Specifically, the user is
no longer responsible for carefully generating sparse systems. For example,
imposing factorization on the ansatz for some scattering amplitude naturally
involves rational functions of kinematics. Combining fractions in order to
obtain sparse equations can become prohibitively complicated and is difficult
to parallelize. Instead, the equations can be extracted by repeatedly sampling
random configurations of the kinematics, _i.e._ , by plugging in random
integers for the generalized Mandelstams. While this will generally produce
dense equations, it is trivial to parallelize, and all of the complications of
solving the system can be foisted onto FiniteFieldSolve. High density
equations can also show up in bootstraps of scalar theories since there is no
helicity structure that breaks the equations into sparse sub-sectors.

While FiniteFieldSolve can natively solve dense systems, some strategies
employed by the solver are not restricted to such systems. Understanding these
strategies requires a brief foray into finite field methods. The standard
approach to exactly solving linear systems is to repeatedly row reduce the
coefficient matrix over different primes and then reconstruct the full
rational solution at the end. Modular arithmetic avoids the round off errors
inherent to floats and is impervious to the intermediary expression swell that
plagues arbitrary precision rationals.222Expression swell should eventually
subside since the final answer is expected to be simple on physical grounds.
However, there are two main drawbacks to finite fields.

First, modular arithmetic introduces a tremendous number of integer division
(or modulo) operations, which is by far the slowest of the four basic
arithmetic operations. This issue can be dramatically improved by “compiling”
the prime. FiniteFieldSolve uses the compiler’s optimizer to convert modulo
operations with respect to a constant prime into much faster combinations of
simpler instructions. This means that the row reduction algorithm must be
_recompiled_ for each prime, but this up front cost is negligible for large
systems. Compiling the prime results in overall performance increasing by up
to a factor of 4 in testing. All solvers that rely on modular arithmetic can
benefit from compiling the primes, but the performance gain will depend on how
much time the solver spends on actual arithmetic as opposed to memory
allocation or pivot selection.

The second issue with finite field methods is that the matrix must be solved
multiple times over different primes in order to reconstruct the full rational
solution. This can lead to extraneous work if large portions of the answer
require only a few primes to reconstruct but certain pieces of the answer
require many primes to capture. FiniteFieldSolve attempts to minimize wasted
effort by gleaning as much information as possible from the first row
reduction. For example, with high probability, it only takes one or a small
number of primes to determine: the linearly independent equations in a system,
the rank of the system, which variables are (in)dependent, which variables are
set to zero, if the system is inconsistent, etc. This information can be
reused over subsequent primes to reduce the size of the matrix that needs to
be row reduced, which can dramatically improve performance as the time to row
reduce an $n\times n$ matrix is naively $\mathcal{O}(n^{3})$. Wisely reusing
information from the first few primes can be applied to all finite field
solvers and is already in use in at least some [25, 27].

Section 2 gives a brief review of finite field methods in the context of
exactly solving linear systems. Technical details of FiniteFieldSolve are
presented in Section 3 along with comments about which techniques in
FiniteFieldSolve may be applicable to other solvers. Besides compiling the
primes and reusing information from early primes, FiniteFieldSolve improves
performance by: using two forms of parallelization, statically allocating
memory for faster access, and using 16-bit primes to improve speed and memory
consumption. Section 4 explains installing FiniteFieldSolve, which amounts to
running an installer script, installing a compiler if necessary, and
downloading an example file if desired. Section 4 also provides a guide to the
most important, high-level functions in FiniteFieldSolve. Example calculations
and benchmarks against other solvers are provided in Section 5. The section
includes a detailed example of bootstrapping the 8pt Dirac-Born-Infeld (DBI)
tree amplitude [4] as well as an overview of constructing the color-dual 4pt
two-loop non-linear sigma model (NLSM) integrand [28]. The examples
demonstrate a workflow for FiniteFieldSolve where constraints are implemented
one at a time and the linear system is periodically pruned of unnecessary
information. Section 6 contains some concluding remarks.

##  2 Finite field methods

Finite fields are a standard tool in exactly solving linear systems so only a
brief review is presented here. For background and related material see Refs
[26, 29, 30, 31, 32, 27, 33, 34, 35, 36, 37] and the references therein.
Finite fields are the method of choice because they are unaffected by
intermediary expression swell and roundoff errors. The tradeoff is that the
solutions produced are only probabilistically correct. At a very high level,
these solvers work by: first row reducing the linear system over different
finite fields $\mathbb{Z}_{p}$, then stitching together the results with the
Chinese remainder theorem (CRT), and then finally reconstructing the rational
solution using the extended Euclidean algorithm (EEA). For convenience, the
relevant number theoretic algorithms are summarized in A.

In more detail, FiniteFieldSolve begins by converting an inhomogeneous system
of equations $A\mathbf{x}=\mathbf{b}$ into a homogeneous one by appending
$\mathbf{-b}$ to $A$ and solving

$\displaystyle(A,\mathbf{-b})\begin{pmatrix}\mathbf{x}\\\
x_{0}\end{pmatrix}=0,$ (1)

where $x_{0}$ is taken to be 1 at the end of the calculation. The system is
inconsistent if solving the system sets $x_{0}$ to a particular value. In this
way any problem can be converted to a homogeneous one of the form
$A\mathbf{x}=0$. The approach to solving this homogeneous system is summarized
in Algorithm 1. In order to solve the system, the matrix $A$ is projected onto
a finite field $\mathbb{Z}_{p}$ where $p$ is prime. Any denominators in $A$
can be projected onto $\mathbb{Z}_{p}$ using Fermat’s little theorem. Of
course, if $A$ contains any factors of $1/p$ then a different prime must be
selected, but this is relatively rare in practice. The main algorithm loops
over primes $p_{1},p_{2}...$ starting at $p_{1}$. Using any desired row
reduction algorithm, the system $A\mathbf{x}_{p_{1}}\equiv 0\mod p_{1}$ is
solved for $\mathbf{x}_{p_{1}}$. (In Algorithm 1 and what follows,
$\mathbf{x}_{a}$ denotes the solution to $A\mathbf{x}_{a}\equiv 0\mod a$ where
$a$ could be any integer.) The algorithm then attempts to reconstruct the
rational solution $\mathbf{x}$ to $A\mathbf{x}=0$ using the extended Euclidean
algorithm (EEA). A fraction $a/b$ will be successfully reconstructed when
$|a|$ and $b$ are smaller than about $\sqrt{p_{1}}$ [38, 39, 40]. Since larger
primes allow more fractions to be reconstructed, conventional wisdom is to use
the largest primes that fit in the word size of the machine.
Counterintuitively, FiniteFieldSolve reaps considerable performance benefits
from choosing slightly smaller primes. In any case, if reconstruction fails,
then the algorithm proceeds to solve $A\mathbf{x}_{p_{2}}\equiv 0\mod p_{2}$
over a second prime. Using the Chinese remainder theorem (CRT), the
information from $\mathbf{x}_{p_{1}}$ and $\mathbf{x}_{p_{2}}$ is combined to
obtain $\mathbf{x}_{p_{1}\cdot p_{2}}$, which solves $A\mathbf{x}_{p_{1}\cdot
p_{2}}\equiv 0\mod(p_{1}p_{2})$ by construction. Using $\mathbf{x}_{p_{1}\cdot
p_{2}}$ is enough to reconstruct rational solutions with numerators and
denominators less than about $\sqrt{p_{1}p_{2}}$. The algorithm row reduces
$A$ over more and more primes until $\mathbf{x}$ solves the original
(rational) problem $A\mathbf{x}=0$. If primes $p_{1},p_{2}...p_{i}$ have been
used, it is often sufficient to check that $\mathbf{x}$ is a solution over the
next prime that would be used, _i.e._ , $A\mathbf{x}\equiv 0\mod p_{i+1}$.
With enough primes it is possible to reconstruct arbitrarily complicated
rational solutions. In practice, typically only a few primes are needed
because the rational numbers in physics problems tend to be simpler than those
that come from solving a random matrix.

Algorithm 1 Solve $A\mathbf{x}=0$

Matrix $A$ and primes $p_{1},p_{2}...$

General solution $\mathbf{x}$ satisfying $A\mathbf{x}=0$

$i\leftarrow 1$

repeat

Solve $A\mathbf{x}_{p_{i}}\equiv 0\mod p_{i}$ for $\mathbf{x}_{p_{i}}$

Build $\mathbf{x}_{p_{1}\cdot p_{2}...p_{i}}$ from CRT on
$\mathbf{x}_{p_{1}\cdot p_{2}...p_{i-1}}$ and $\mathbf{x}_{p_{i}}$

Reconstruct rational $\mathbf{x}$ from EEA on $\mathbf{x}_{p_{1}\cdot
p_{2}...p_{i}}$

$i\leftarrow i+1$

until $A\mathbf{x}=0$

##  3 Technical details

FiniteFieldSolve implements the solver described in the previous section
through two distinct components: a high-level frontend and a low-level
backend. The backend is responsible for performing the actual row reduction
over a prime. As this is the most computationally intensive part of the
algorithm, the backend is written in C++ to be as performant as possible. The
frontend of FiniteFieldSolve is a Mathematica wrapper that interfaces with the
row reducer and is responsible for reconstructing the full rational solution.
The core structure of the front end is based on spaSMLink [21] but it has been
dramatically altered to include the features mentioned in this section.

The performance of the backend comes from leveraging modern CPU architecture
and optimizing compilers. The C++ backend employs a modified version of
ordinary Gauss-Jordan elimination with unsigned 16-bit integers and is
_recompiled_ each time it is used. These unusual design features are centered
around the following five performance considerations.

  1. 1.

Gauss-Jordan elimination is a highly parallelizable algorithm since many row
operations can be performed independently. Multi-core parallelization comes
with a small overhead but the testing in Section 5 demonstrates that the
benefits outweigh the costs, sometimes by a generous margin. FiniteFieldSolve
uses OpenMP for parallelization where more details are given in Section 4.
Unlike more specialized algorithms, Gauss-Jordan elimination is highly
predictable, so it is easy to reliably bound the time it will take to solve a
system, at least over a single prime.

  2. 2.

Gauss-Jordan elimination also benefits from vectorization where multiple
elements in a row are updated simultaneously using the SIMD instruction set
extension of the processor. For x86 processors, AVX2 is the relevant
instruction set extension since it provides SIMD instructions for _integer_
registers.

  3. 3.

Representing the matrix over a finite field incurs a tremendous number of
modulo (or division) operations during row reduction. Since division is both
slow and common in this algorithm, it is a prime target for optimization. The
frontend of FiniteFieldSolve automatically recompiles the backend each time
the row reduction routine is called, effectively making the prime a compile-
time constant. The optimizer can then turn divisions with respect to the prime
into much faster combinations of instructions, significantly improving overall
performance [41].333If the divisor is known at compile time, gcc uses the
methods of Ref [41] to optimize the division. Ref [42] proposed a faster
method if just the remainder of the division is required, but gcc does not
appear to make use of this approach. If the divisor is not known at compile
time, there are still methods for optimizing the divisions [43]. Mersenne and
Mersenne-like primes compile down to simpler instructions, but this does not
appear to result in improved real-world performance, likely due to cache
misses and memory bandwidth limitations. As will be shown in Section 5, the
single decision to compile the prime leads to overall performance increasing
by up to a factor of 4.

  4. 4.

By virtue of recompiling the backend for each row reduction, the matrix
dimensions are also known at compile time, allowing the compiler to statically
allocate memory for the matrix. Statically allocated memory is typically much
faster to manipulate than dynamically allocated memory. Forcing the backend to
dynamically allocate the matrix results in overall performance dropping by a
factor of 2 to 3 in testing.

  5. 5.

Unsigned 16-bit integers are used for representing elements of the finite
field. 16-bit integers have the advantage that the product of two such
integers fits in a machine word, avoiding complicated overflow logic, and 16
bits is enough to reconstruct many of the simple fractions like 1/2 that
pervade physics. Choosing small integers improves the memory footprint of
FiniteFieldSolve – since multiple integers can be packed into a machine word –
and thereby improves performance by alleviating memory bandwidth bottlenecks.
16-bit integers also happen to be faster than their 32-bit counterparts.
Because 32-bit integers are twice as long, it typically requires half as many
primes to reconstruct the same rational solution. To make up for this, 32-bit
integers would only need to be half as fast as 16-bit integers, but this is
not the case in testing. Likewise, 8-bit integers would need to be twice as
fast as 16-bit integers, but again this is not borne out in testing. 16-bit
integers thus achieve the best balance of speed and size. FiniteFieldSolve
uses 16-bit integers by default but there are options to use 32-bit integers
if desired.

The five considerations above contribute substantially to the performance of
the backend but there are several features of the Mathematica frontend that
also increase performance. Understanding these features requires a brief
detour into fill-in and pivoting. The backend uses the first available pivot,
working its way from left to right and top to bottom in the matrix. This has
the advantage that the algorithm cannot stall in the process of searching for
optimal pivots as can occur for more complicated pivoting approaches. However,
this naive pivoting algorithm may incur serious fill-in. As a simple but
effective measure to mitigate fill-in, the frontend sorts the rows of the
matrix by their density so that the sparsest rows are at the top of the
matrix. Although fill-in will not negatively impact memory usage, as it would
in a sparse algorithm, minimizing fill-in improves performance since entire
rows can be skipped if they have a zero in the pivot column. Said another way,
the frontend attempts to avoid situations where a very dense row at the top of
the matrix would contaminate subsequent sparse rows. The frontend also
improves performance by using information from the first few primes to shrink
the matrix it must solve over subsequent primes. Specifically, the frontend
uses the first prime to find the linearly independent equations and only
solves that subset of equations over the remaining primes. Last, the row
reduction over the first prime provides a good indication of which variables
the equations set to zero, allowing the frontend to limit the number of
columns used with future primes.

So far, this section has mostly focused on the techniques used to improve the
speed of FiniteFieldSolve. The discussion would not be complete without an
overview of the memory layout of the package as memory can become the limiting
resource for large matrices. Care was taken in FiniteFieldSolve to reduce the
memory footprint and produce stable, predictable memory usage, so that the
package is predictable in both time and memory. The majority of the memory is
used for representing three matrices. First, Mathematica stores the original
input matrix as a List or SparseArray of arbitrary precision rationals. The
total memory used to represent this matrix varies based on its dimensions and
density but is constant in time. The second matrix is a copy of the first one
after reduction modulo a prime $p$. The C++ backend operates on this second
matrix. To conserve memory, the matrix is passed from the frontend to the
backend one row at a time. When FiniteFieldSolve is used in its default 16-bit
mode, this $m\times n$ matrix requires $2mn$ bytes of memory. Since the solver
tries to eliminate unnecessary rows and columns when operating on subsequent
primes, the memory needed to store the backend matrix tends to decrease over
time. For anything except the densest input systems, the majority of the
memory will be used by the 16-bit backend matrix. The third and final matrix
is used to represent the rational solution to the input system.444At various
stages in the calculation the solution is actually made of multiple matrices
that are combined using the Chinese remainder theorem. However, this does not
alter the overall picture of the memory usage. This third matrix is stored in
Mathematica as a SparseArray of arbitrary precision rationals. Since a generic
solution looks schematically like a diagonal of ones with some number of non-
zero columns, this matrix tends to be very sparse and consumes the least
memory of any of the matrices. However, the memory used to store this matrix
will increase slowly over time as the rational numbers inside of it become
more complex. Aside from a small overhead, the total memory consumption is
very close to the sum of the memory used to represent each of these three
matrices.

As a final note, FiniteFieldSolve is not intended for very small systems since
the backend needs to be recompiled for each prime. The overhead from calling
the compiler is independent of the matrix size so this only negatively impacts
small matrices, which is perfectly acceptable since these matrices are already
trivial to solve. Once the time to solve the system overshadows the compile
time then FiniteFieldSolve becomes more performant than Mathematica’s built in
Solve or Reduce. Depending on the computer, the tipping point occurs in the
high hundreds to low thousands of parameters when the time to solve the matrix
takes of order one second.

##  4 User guide

This section describes how to install FiniteFieldSolve along with the major
functions included in the package. The package is available at

https://github.com/jfmangan/FiniteFieldSolve.

The repository contains four relevant files: InstallScript.m, Examples.m,
FiniteFieldSolve.m, and RowReduceLink.cpp. Rather than cloning the repository
with git, it is recommended to download files individually as described below.
Download Examples.m and InstallScript.m and place them in any desired
directory. As the name suggests, Examples.m provides several examples that can
be run once FiniteFieldSolve has been properly installed. One of the included
examples is discussed in Section 5. To install FiniteFieldSolve run
InstallScript.m from any directory. This script will create a folder called
FiniteFieldSolve in $UserBaseDirectory/Applications and then it will download
FiniteFieldSolve.m and RowReduceLink.cpp and place them inside the new
directory.555On startup, FiniteFieldSolve.m will automatically determine its
location in the filesystem, so it is possible for users to install the package
outside of $UserBaseDirectory. However, it is still recommended to use
InstallScript.m as this will ensure that Mathematica can load the package.
Running FiniteFieldSolve requires a C++ compiler. Even though FiniteFieldSolve
automatically detects the operating system, setting up the compiler varies
between machines.

  * 1.

_Windows_ is not natively supported. In order to run FiniteFieldSolve on
Windows it is suggested to use Windows Subsystem for Linux (WSL) [44], Cygwin
[45], or MinGW [46], at which point the Windows installation becomes the same
as for Linux. FiniteFieldSolve has been tested with WSL2.

  * 2.

_Linux_ ships with gcc so when FiniteFieldSolve detects Linux as the operating
system, this is the compiler it will use.

  * 3.

_MacOS_ does not ship with a compiler so one must be installed. The simplest
option is to use clang. If it is not already installed, it can be installed by
running

xcode-select --install

in a terminal. By default, FiniteFieldSolve will use clang on Mac. Using a
different compiler, like gcc, requires modifying the ‘gccString’ in
FiniteFieldSolve.m. FiniteFieldSolve has been tested on Macs with both Intel
and Apple silicon processors.

Although OpenMP is recommended, it is disabled by default because its proper
configuration is beyond the scope of this installation guide. If OpenMP is
installed and the compiler and environment variables are configured properly,
OpenMP can be enabled by changing the ‘IsOpenMPInstalled’ flag in
FiniteFieldSolve.m from False to True. The speed improvements from OpenMP are
discussed in Section 5. Once FiniteFieldSolve has been installed, it can be
loaded in Mathematica with

<<FiniteFieldSolve`

just like any other package. With installation complete, it is time to discuss
the three most important high-level functions in the package.

FiniteFieldSolve[equations, (optional: optionsList)] takes a list of equations
and returns the solution as a list of replacement rules. For example,

FiniteFieldSolve[{a==b, a==1}]

returns {a$\rightarrow$1, b$\rightarrow$1}. FiniteFieldSolve will
automatically detect the variables in the equations. If an equation does not
contain ‘==’, then it is assumed that the equation is equal to zero.
optionsList is an optional parameter of the form {string1, string2...} where
the strings are treated as case independent. FiniteFieldSolve will print more
output if ‘verbose’ is included in optionsList. By default the equations are
sorted by their density before they are solved. To disable this feature
include ‘NoRowSorting’ in optionsList. By default FiniteFieldSolve will use
the first prime it solves over to determine the linearly independent rows and
only use those rows when solving over subsequent primes. This behavior can be
disabled by including ‘KeepLinearDepRows’ in optionsList. By default
FiniteFieldSolve will use the first prime to detect if the equations set any
variables to zero. These variables will not be solved for in the future and
the ensuing linearly dependent rows will be removed by row reducing over the
second prime. To disable this behavior include ‘KeepZeroVariables’ in
optionsList. Row reduction will default to 16-bit primes which are generally
faster and use less memory. To use 32-bit primes include ‘32bit’ in
optionsList. For dense equations – such as those generated by random numerical
sampling over the rationals – it can be useful to clear denominators to avoid
$1/p$ factors. To clear all denominators add ‘ClearDenominators’ to
optionsList.

FindLinearlyIndependentEquations[equations, (optional: prime), (optional:
optionsList)] takes a list of equations and returns the linearly independent
ones. The advantage of this function is that it is typically faster than
solving the whole system. When building up a complicated system of equations,
it can be helpful to periodically prune unnecessary equations with this
function. Failing to do so can result in a final system that is vastly
overcomplete, increasing the time and (peak) memory needed to solve it. This
function also gives an accurate estimate of the rank of a system. If equations
are being generated by numerical sampling, then the number of variables minus
the rank gives an upper bound on the number of equations that must be
generated. By default the function uses the prime 65,521, the largest 16-bit
prime, but a different prime can be specified. This function also defaults to
sorting the equations by density before row reduction so that the function
returns the sparest linearly independent equations. This behavior can be
disabled by including ‘NoRowSorting’ in optionsList.

ConsistentEquationsQ[equations, (optional: prime)] takes an inhomogeneous
system of equations and returns False if the system is inconsistent and True
otherwise. This function operates by row reducing the system over a given
prime, which defaults to 65,521, and then checking if the row reduction is
inconsistent. If the equations are inconsistent over the prime, then they are
guaranteed to be inconsistent over the rationals. To detect an inconsistency,
it is faster to run this function than to reconstruct the full rational
solution since the former only requires a single prime but the latter
typically requires several. For this reason, this function can be a useful
diagnostic tool when an inconsistent system is expected.

The three functions described above share certain features in common. They all
internally convert the equations to a matrix (in the form of a SparseArray) by
calling CoefficientArrays. The matrix may take up less memory than the
equations themselves, so for advanced usage users may wish to work directly
with the matrix. In this case, the three relevant functions are
FiniteFieldSolveMatrix, FindLinearlyIndependentRows, and ConsistentMatrixQ.
The precise documentation for these functions can be found in
FiniteFieldSolve.m.

The workhorse behind all of the functions mentioned above is
RowReduceOverPrime which takes a matrix and row reduces it over a given prime.
Critically, the rows are never rearranged during row reduction, rows of zeros
are never deleted, and the first available pivot is always selected when
scanning the matrix from left to right and top to bottom. This means that
RowReduceOverPrime does not return a matrix in strict row echelon form.
However, this algorithm has the advantage that its output can be used to
trivially identify the linearly independent rows in the input
matrix.666Specifically, if row $i$ in the row reduced matrix is a row of
zeros, then row $i$ in the original matrix represents a linearly dependent
equation. Furthermore, these rows will come from the top of the matrix where
FiniteFieldSolve stores the sparsest equations. So when solving a full system
over multiple primes FiniteFieldSolve can use RowReduceOverPrime to find the
linearly independent, sparsest equations over the first prime to speed up row
reduction over later primes. Details for RowReduceOverPrime can be found
inside FiniteFieldSolve.m.

##  5 Worked example and benchmarks

This section describes two applications of FiniteFieldSolve to bootstrapping
relativistic scattering amplitudes. The first example shows how to determine
the 8pt tree amplitude for Dirac-Born-Infeld (DBI) theory using its soft
behavior [4]. This calculation is worked through in full detail in the example
file included with FiniteFieldSolve. The second example is an overview of
determining the 4pt two-loop color-dual integrand for the non-linear sigma
model (NLSM) [28]. FiniteFieldSolve is benchmarked against other solvers in
these examples.

###  5.1 Soft bootstrap of 8pt DBI

DBI is a massless scalar field theory that is uniquely determined by its
infrared behavior [4]. With a mostly plus metric signature and in units where
the coupling constant is unity, the DBI Lagrangian is

$\displaystyle\mathcal{L}=-\sqrt{1+(\partial\phi)^{2}}\,.$ (2)

In order to obtain a non-trivial linear system, we will bootstrap the 8pt tree
amplitude $A_{8}$ assuming that the 4pt and 6pt amplitudes have already been
computed using the bootstrap. The amplitude is built from three separate types
of graphs, two of which have factorization channels and the third represents a
contact term

The two graphs on the left can be obtained by re-using the amplitudes obtained
from earlier stages in the bootstrap, leaving only the contact diagram to be
determined. From the DBI Lagrangian, the contact term should be proportional
to

$\displaystyle(p_{1}p_{2})(p_{3}p_{4})(p_{5}p_{6})(p_{7}p_{8})+\text{perms}$
(3)

where $p_{i}$ is the momentum of particle $i$ and the sum extends over all
permutations of external labels. However, we will confirm by explicit
calculation that the contact term is uniquely fixed by the soft behavior of
DBI without any reference to the Lagrangian.

The ansatz for the contact term is constructed as follows. There are 20
generalized Mandelstam invariants $p_{i}\cdot p_{j}$ that are independent on-
shell, that is, after imposing

$\displaystyle\left\\{p_{i}^{2}=0~{}|~{}1\leq i\leq n\right\\}\text{ and
}\sum\limits_{i=1}^{n}p_{i}=0\,.$ (4)

The exact basis chosen does not matter. The basis in the example file includes
$(p_{2}p_{4})$ and $(p_{2}p_{5})$ as well as 18 more dot products. There are
8,855 independent terms with mass dimension $p^{8}$ so the contact term looks
schematically like

$\displaystyle
A_{8,\,\text{contact}}(p_{1},p_{2},...p_{8})=c_{1}(p_{2}p_{4})^{4}+c_{2}(p_{2}p_{4})^{3}(p_{2}p_{5})+...$
(5)

The first constraint on this ansatz is that the contact term should be Bose
symmetric, that is, the contact term should be invariant under any relabeling
of the legs. Full permutation invariance can be guaranteed by imposing two
independent symmetry relations. First, the contact term should be invariant
under cyclically rotating all of the labels by one position,

$\displaystyle
A_{8,\,\text{contact}}(p_{1},p_{2},p_{3},p_{4},p_{5},p_{6},p_{7},p_{8})=A_{8,\,\text{contact}}(p_{2},p_{3},p_{4},p_{5},p_{6},p_{7},p_{8},p_{1}),$
(6)

and, second, the contact term should be invariant under swapping particles 1
and 2,

$\displaystyle
A_{8,\,\text{contact}}(p_{1},p_{2},p_{3},p_{4},p_{5},p_{6},p_{7},p_{8})=A_{8,\,\text{contact}}(p_{2},p_{1},p_{3},p_{4},p_{5},p_{6},p_{7},p_{8}).$
(7)

Naively each of the conditions in (6) and (7) produces 8,855 equations – as
many equations as there are terms in the ansatz – but some of the equations
become trivial in the right Mandelstam basis. Table 1(c) documents how long it
takes FiniteFieldSolve to solve the linear system along with the dimensions
and density of the system. The table also contains timings against
Mathematica’s Solve, LinearSystemSolver [22], spaSMLink [21], and Kira 2 [25].
The testing was done on two different machines, one against spaSMLink and one
against Kira. When comparing to spaSMLink and Kira it is worth remembering
that these packages are intended for very sparse systems. Kira does not
directly interface with Mathematica but the time to convert between the two
programs was not included in the timings in the table. In the benchmarks
presented in Table 1(c), OpenMP was enabled for FiniteFieldSolve on both
machines. Disabling OpenMP has virtually no effect on overall performance in
this instance, but this will not be the case for the next example. By
declaring the prime in the C++ backend as ‘volatile’ it is possible to
estimate the performance that would be lost if FiniteFieldSolve did not
compile its primes. Failing to compile the primes results in the overall solve
time increasing by roughly a factor of 2. The time savings from compiling the
prime will be more dramatic in the next example. As the table demonstrates,
FiniteFieldSolve is about two orders of magnitude faster than Mathematica’s
Solve and uses about an order of magnitude less memory. In this setting,
FiniteFieldSolve is an order of magnitude faster than the two sparse solvers
and uses the same order of magnitude of memory. These results are not
surprising given the intended use case of the sparse solvers.

Table 1: FiniteFieldSolve was benchmarked against other solvers in
bootstrapping the 8pt DBI contact term. Note that spaSMLink and Kira are
intended for very sparse systems. The ratios in the table are with respect to
the time and memory it takes FiniteFieldSolve to solve the system. (a)
Benchmarks on MacOS with i5-7267U and 16GB memory

| FiniteFieldSolve | Mathematica Solve | LinearSystemSolver | spaSMLink
---|---|---|---|---
Time | 11 s | 1400 s | 2200 s | 1200 s
Time ratio | 1 | 130 | 200 | 110
Memory | 510 MB | 9.7 GB | 5.6 GB | 980 MB
Memory ratio | 1 | 19 | 11 | 1.9

(b) Benchmarks on Ubuntu under WSL2 with i5-6600 and 32GB memory

| FiniteFieldSolve | Mathematica Solve | LinearSystemSolver |  Kira 2
---|---|---|---|---
Time | 8.5 s | 1500 | 1900 s | 180 s
Time ratio | 1 | 180 | 220 | 21
Memory | 530 MB | 24 GB | 5.7 GB | 210 MB
Memory ratio | 1 | 46 | 11 | 0.40

(c) Matrix properties

Matrix dimensions | 16,705 $\times$ 8,855
---|---
Matrix density | $3.4\times 10^{-3}$
Highest row density | $2.4\times 10^{-2}$

After imposing Bose symmetry there are only 5 free parameters left in the
contact term. All remaining freedom in the ansatz can be fixed by constructing
$A_{8}$ from (5.1) and imposing the soft theorem of DBI.777There are other
principles that fully fix the amplitude including classical conformal
invariance [47]. Amplitudes in this theory vanish as momentum squared in the
soft limit, that is,

$\displaystyle\lim\limits_{z\to 0}A_{n}(zp_{i})=\mathcal{O}\left(z^{2}\right)$
(21)

where particle $i$ is taken soft. At a practical level, when taking the soft
limit of an on-shell amplitude, it is key that $p_{i}$ does not play any
special role in enforcing the on-shell conditions. For example, if momentum
conservation was used to eliminate $p_{n}$ from $A_{n}$, this does not mean
that $A_{n}$ is identically zero in the soft limit of $p_{n}$. After scaling
the momentum of one of the unprivileged particles via $p_{i}\to zp_{i}$, there
are multiple methods to extract the actual equations enforcing the soft limit.
One option would be to call Together on the entire amplitude. A more
appropriate method in this case is to sample different kinematical
configurations numerically by repeatedly mapping the Mandelstams
$(p_{j}p_{k})$ to random integers. Each choice of random Mandelstams will
generally generate a different fully dense equation enforcing the soft limit.
One of the convenient features of FiniteFieldSolve is that it can easily
handle these dense equations, sidestepping an expensive call to Together.
However, since there are only 5 remaining parameters, Mathematica’s Solve is
more than sufficient to complete the bootstrap. After imposing the soft limit,
the amplitude is fully fixed and the contact term corresponds to the 8pt
interaction in (2).

###  5.2 Color-dual bootstrap of 4pt two-loop NLSM

The next example is the construction of color-dual 4pt two-loop numerators for
the non-linear sigma model (NLSM). This example is more complex than the
previous one, but a brief overview is provided below. For the full calculation
see Ref [28]. The benchmarks in Table 2(c) can be understood without any
knowledge of how the system was constructed. The following discussion assumes
some familiarity with color-kinematics duality and the double copy [2, 3].
Recent reviews of these topics include Refs [48, 49, 50].

Color-kinematics duality applies to certain special theories, including Yang-
Mills and NLSM, that carry color (or flavor) indices. The $n$-pt $L$-loop
scattering amplitude in any colored theory can be written as

$\displaystyle\mathcal{A}_{n}^{L}=\sum\limits_{\Gamma}\frac{1}{S_{\Gamma}}\int\left(\prod\limits_{i=1}^{L}\frac{d^{D}\ell_{i}}{(2\pi)^{D}}\right)\frac{C_{\Gamma}N_{\Gamma}}{d_{\Gamma}}$
(22)

where $D$ is the spacetime dimension and the sum runs over _cubic_ graphs
$\Gamma$. The amplitude can always be written as a sum over purely cubic
graphs by multiplying and dividing by propagators. Here $S_{\Gamma}$ is the
symmetry factor associated with the graph $\Gamma$, $d_{\Gamma}$ is the
product of propagators, $C_{\Gamma}$ is the color factor made of a product of
structure constants $f^{abc}$, and $N_{\Gamma}$ is the kinematic numerator
made of dot products of momentum and polarization vectors. The color factors
in (22) obey a set of Jacobi identities. Color-kinematics duality states that
there exists a way to write the kinematic numerators so that they obey the
same Jacobi relations as the color factors. Only a few Lagrangians are known
to manifest color-kinematics duality [28, 51, 52, 53, 54, 55, 56, 57, 58], so
color-dual numerators are typically determined using an ansatz. Once color-
dual numerators are found, the amplitude for the “double-copy” theory is
obtained by replacing $C_{\Gamma}$ with another copy of $N_{\Gamma}$ in (22).

As a brief example of color-kinematics duality, consider the NLSM 4pt tree
amplitude.888Tree-level amplitudes in NLSM are fully fixed by color-kinematics
duality [59]. The only graph topology contributing to this process is

$\displaystyle{\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\feynman \vertex(v0) at (0,0);
\vertex(v1) at (1,0); \vertex(p1) at (-0.707,-0.707){1}; \vertex(p2) at
(-0.707,0.707){2}; \vertex(p3) at (1+0.707,0.707){3}; \vertex(p4) at
(1+0.707,-0.707){4}; \diagram{ (v0) -- [ultra thick,](v1), (v0) -- [ultra
thick,](p1), (v0) -- [ultra thick,](p2), (v1) -- [ultra thick,](p3), (v1) --
[ultra thick,](p4), }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}\,.$ (23)

The color factor associated with this ($s$-channel) graph is
$C_{(12|34)}=f^{12c}f^{c34}$. The color-dual form of the pion numerator can be
found by making an ansatz with the appropriate power counting,

$\displaystyle N_{(12|34)}=c_{1}s^{2}+c_{2}st+c_{3}t^{2},$ (24)

where $s=(p_{1}+p_{2})^{2}$, $t=(p_{2}+p_{3})^{2}$, and $u=-s-t$. This ansatz
implicitly assumes that the numerator is a local function, which is not
absolutely necessary for color-kinematics duality and can be an overly
restrictive assumption at times [28, 54, 60, 61, 62]. Aside from locality and
power counting, three explicit constraints are imposed on the ansatz. First,
the numerator must respect all of the symmetries of its associated graph. At
4pt tree level the graph symmetries are spanned by

$\displaystyle N_{(12|34)}$ $\displaystyle=-N_{(21|34)}$ (25) $\displaystyle
N_{(12|34)}$ $\displaystyle=N_{(34|12)},$ (26)

where the sign compensates for the signs buried in the structure constants
making up the color factor. Strictly speaking color-kinematics duality does
not require that numerators respect their graph symmetries, but it is a
convenient and natural requirement just like postulating local ansatzë for
numerators. Another benefit of enforcing graph symmetries is that the $t$\-
and $u$-channel numerators can be obtained by relabeling the $s$-channel
numerator $N_{(12|34)}$. The second constraint is that the kinematic
numerators should obey the same Jacobi relations as the color factors. In this
case, the color factors obey

$\displaystyle C_{(12|34)}+C_{(23|14)}+C_{(31|24)}=0$ (27)

so the numerators are required to satisfy the same condition

$\displaystyle N_{(12|34)}+N_{(23|14)}+N_{(31|24)}=0\,.$ (28)

The third and final requirement is that the amplitude factorizes correctly on
its poles, or in general, that the integrand has the correct unitarity cuts
[63, 64, 65]. For 4pt NLSM this is trivial since the on-shell 3pt amplitude
vanishes on general grounds. Together, (25) and (28) fix the ansatz in (24) up
to an overall coupling constant, giving

$\displaystyle N_{(12|34)}\propto s(t-u)\,.$ (29)

The fully color-dressed pion amplitude is

$\displaystyle
A_{4}=\frac{C_{(12|34)}N_{(12|34)}}{(p_{1}+p_{2})^{2}}+\text{cyclic}(1,2,3),$
(30)

where the sum extends over cyclic permutations of labels 1, 2, and 3. With the
color-dual numerators in hand, the double copy produces the 4pt special
Galileon amplitude

$\displaystyle
A_{4}=\frac{N^{2}_{(12|34)}}{(p_{1}+p_{2})^{2}}+\text{cyclic}(1,2,3)\propto
stu\,.$ (31)

For the case of the 4pt two-loop NLSM integrand constructed in Ref [28], power
counting requires that the numerator scales as $p^{12}$. There are 8,008
independent scalar invariants at this mass dimension and since two diagrams
require ansatzë, the entire problem involves 16,016 parameters, at least for
the case of a local numerator.999For technical reasons there are actually two
more parameters in the system, one of which is used to represent inhomogeneous
terms in equations as explained in Section 2. The key constraints are the same
as they were for the 4pt tree example, namely, the numerators must respect
graph symmetries, the numerators must obey the same linear relations as the
color factors, and the integrand must factorize correctly. As a matter of
convenience, the Jacobi relations are imposed first so that the basis of
numerators can be reduced to just two graphs. The Jacobi relations produce
59,544 equations, where the exact number depends on implementation details.
The linearly independent subset of 11,543 equations can be selected using
FindLinearlyIndependentEquations. Solving the equations is unnecessary at this
stage and will actually overcomplicate the problem as this intermediary
solution is more complex than the final solution. Next, imposing graph
symmetries produces 115,494 equations, which are similarly implementation-
dependent. Combining the graph symmetry equations with the Jacobi relations
produces 14,773 linearly independent equations. Yet again, it is unnecessary
to actually solve the equations; finding the linearly independent subset is
enough.

The third key constraint is that the integrand has the correct unitarity cuts.
While the Jacobi relations and the graph symmetries are polynomial
constraints, the unitarity cuts generate constraints that are sums of rational
functions where the denominators come from uncut momenta. Extracting the
actual equations requires combining fractions with something like Together in
Mathematica. While this is straightforward in theory, this can become a severe
bottleneck in practice. Alternatively, the equations can be extracted by
plugging in random integers for the kinematics appearing in the rational
functions. Each choice of random kinematics yields another equation. Since
FindLinearlyIndependentEquations effectively computes the rank of the
equations already generated, it is possible to accurately estimate the number
of equations that need to be produced numerically. In this case, 16,016 $-$
14,773 $=$ 1,243 equations that need to be generated by random sampling. Note
that it is trivial to parallelize numerical sampling but the same is not true
for combining large rational functions. The most important advantage of
numerical sampling is that FiniteFieldSolve allows the user to offload all of
the complications onto the computer; no human effort needs to be put into
finding or maintaining the sparse equations that come from combining rational
functions. The disadvantage of sampling is that, first, it produces dense
equations and, second, it requires a method for computing the rank.
FiniteFieldSolve addresses both issues. The rank of the graph symmetry and
Jacobi equations could be obtained by directly solving them, but this is both
overkill and slow since the intermediary solution requires more primes to
reconstruct than the final solution. Aside from computing the rank,
periodically culling the system of linearly dependent equations also has the
advantage that it reduces the memory required to solve the final system. After
imposing unitarity cuts, the system has 15,251 equations. Since there are
16,016 ansatz parameters, there is still room to enforce additional aesthetic
constraints that simplify the answer. The last constraints come from matching
the NLSM numerators to those of Zakharov-Mikhailov theory in two spacetime
dimensions [55, 66].

Table 2: FiniteFieldSolve was benchmarked against other solvers in
bootstrapping the 4pt two-loop NLSM integrand. Note that spaSMLink and Kira
are intended for very sparse systems. The ratios in the table are with respect
to the time and memory it takes FiniteFieldSolve to solve the system. (a)
Benchmarks on MacOS with i5-7267U and 16GB memory

| FiniteFieldSolve | spaSMLink
---|---|---
Time | 380 s | 4.4 hr
Time ratio | 1 | 42
Memory | 830 MB | 4.4 GB
Memory ratio | 1 | 5.4

(b) Benchmarks on Ubuntu under WSL2 with i5-6600 and 32GB memory

| FiniteFieldSolve |  Kira 2
---|---|---
Time | 420 s | 67 h
Time ratio | 1 | 570
Memory | 1.0 GB | 2.6 GB
Memory ratio | 1 | 2.6

(c) Matrix properties

Matrix dimensions | 15,959 $\times$ 16,018
---|---
Matrix density | $6.4\times 10^{-2}$
Highest row density | 1

Once all of the equations have been assembled, it is finally time to solve
them. Since the system was assembled with FiniteFieldSolve, all of the
linearly dependent equations have been removed, leaving only the essential
equations left to solve. FiniteFieldSolve was timed against spaSMLink and Kira
in solving the system and the results are presented in Table 2(c). Again,
recall that spaSMLink and Kira are intended for very sparse systems, which is
not the case for the fully dense rows generated by numerically sampling the
unitarity constraints. Even though the linear system has been timed on
different solvers, FiniteFieldSolve was instrumental in constructing it.
Similar to the DBI example, FiniteFieldSolve uses a comparable amount of
memory to spaSMLink and Kira and is one to two orders of magnitude faster than
the sparse solvers. For the DBI calculation, Kira was faster than spaSMLink,
but the opposite is true here. This may be due to the different pivoting
algorithms they use. Since the NLSM system contains roughly twice as many
parameters as the DBI system, it might be reasonable to expect that it would
take very roughly $2^{3}$ as much time to solve. The NLSM system takes longer
to solve than expected because the solution requires more primes to
reconstruct and the system is denser so FiniteFieldSolve cannot skip as many
rows during row reduction. OpenMP was enabled for FiniteFieldSolve when
constructing the results in Table 2(c). If OpenMP is disabled, the solve times
increase by roughly a factor of 1.4 on both machines, which is reasonable
since both computers have dual-core CPUs. More importantly, if the prime is
not compiled, then overall performance drops by roughly a factor of 4 on both
computers.

##  6 Discussion

FiniteFieldSolve is a general-purpose solver for exact linear systems over the
rationals. Systems of equations of this type appear prominently in high energy
theory in integral reduction routines and amplitude bootstraps. A worked
example of a bootstrap problem was provided in the text and accompanying
repository. FiniteFieldSolve uses well-known finite field methods to avoid
intermediary expression swell and roundoff errors that might otherwise plague
a traditional solver. In testing, the new package is roughly two orders of
magnitude faster and uses an order of magnitude less memory than Mathematica,
opening the door to new physics problems.

A variety of techniques are used to improve performance in FiniteFieldSolve.
First, the C++ backend is recompiled every time it is called. Since the matrix
dimensions are known at compile time, the memory can be statically allocated,
improving memory access. This technique is confined to dense solvers where
memory for new nonzero values never needs to be allocated. Compiling the
backend also converts the ubiquitous and expensive modulo operations into a
much faster combination of instructions. Compiling the prime leads to overall
speed increasing by up to a factor of 4 in testing. Generally, any piece of
software that performs a lot of modular arithmetic with respect to a fixed
prime may benefit from re-compiling for each prime. The actual time savings
will depend on the proportion of time that the software spends on arithmetic.
The backend also improves memory usage and speed by using 16-bit primes, where
32-bit and 8-bit alternatives are slower in testing. Due to its simple row
reduction algorithm, the package benefits from two layers of parallelization.
Moving to the frontend, it improves performance by extracting as much
information as possible over the first few primes, including finding the
linearly independent equations of the system. This can sometimes dramatically
simplify the system that must be solved over subsequent primes. The ability to
easily identify linearly independent rows in conjunction with the solver’s
ability to efficiently manipulate dense equations frees the user to generate
equations by any and all means, including numerical sampling.

Although some of the techniques in FiniteFieldSolve may be applicable to other
solvers, FiniteFieldSolve is most directly suited to dense systems including
those that appear in certain amplitude bootstraps. One of the most promising
applications of FiniteFieldSolve is in bootstrapping higher-dimension color-
dual operators since this can involve large dense systems, especially for
scalar theories [67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81].
Other promising applications may include: bootstrapping higher multiplicity
AdS amplitudes [82, 83, 84, 85, 86, 87]; bootstrapping amplitudes from their
physical properties like soft behavior, gauge invariance, factorization,
conformal invariance, etc., [4, 88, 47, 89, 90, 91, 92, 93, 94, 95]; and
enumeration of higher-dimension or higher-multiplicity effective field theory
operators [96].

Last, and perhaps most importantly, FiniteFieldSolve is simple to use and
straightforward to install. FiniteFieldSolve and Mathematica’s Reduce and
Solve share similar syntax so it is trivial to swap one out for the other.
Installing FiniteFieldSolve on Linux amounts to running the installer script
included in the repository. The procedure is identical on MacOS except that
the user may need to install a compiler by running a single command in the
terminal. FiniteFieldSolve has been tested under Windows, but the installation
process is more involved.

#### Acknowledgments

The author would like to thank JJ Carrasco, Sasank Chava, Alex Edison, Kezhu
Guo, and Nic Pavao for feedback on the manuscript and software, insightful
conversations, and related collaboration. This work was supported by the DOE
under contract DE-SC0015910 and by the Alfred P. Sloan Foundation. J.M.
additionally acknowledges the Northwestern University Amplitudes and Insight
group, the Department of Physics and Astronomy, and Weinberg College for their
generous support. Feynman diagrams in this paper were typeset using TikZ-
Feynman [97].

##  Appendix A Number theoretic algorithms

This appendix provides a brief summary of the number theoretic algorithms
mentioned in the text including the extended Euclidean algorithm, Fermat’s
little theorem, and the Chinese remainder theorem.

Given integers $a$ and $b$, the Euclidean algorithm is a method for computing
integers $x$ and $y$ such that $ax+by=(a,b)$, where $(a,b)$ denotes the
greatest common divisor of $a$ and $b$. The extended variant of the Euclidean
algorithm starts by setting

$\displaystyle r_{0}=a\quad\quad\quad$ $\displaystyle r_{1}=b$ $\displaystyle
s_{0}=1\quad\quad\quad$ $\displaystyle s_{1}=0$ (45) $\displaystyle
t_{0}=0\quad\quad\quad$ $\displaystyle t_{1}=1.$

The algorithm proceeds by repeated integer division

$\displaystyle q_{i}=\left\lfloor\frac{r_{i-1}}{r_{i}}\right\rfloor$ (46)
$\displaystyle r_{i+1}=r_{i-1}-q_{i}r_{i}$ (47) $\displaystyle
s_{i+1}=s_{i-1}-q_{i}s_{i}$ (48) $\displaystyle t_{i+1}=t_{i-1}-q_{i}t_{i}$
(49)

and terminates when $r_{n+1}=0$. The values of $x$ and $y$ are given by
$s_{n}$ and $t_{n}$, respectively. Since $r_{n}=(a,b)$, the algorithm can be
seen as a method for determining the greatest common divisor. Furthermore, the
extended Euclidean algorithm can be used for rational reconstruction of $a$
over $\mathbb{Z}_{b}$. By induction, $as_{i}+bt_{i}=r_{i}$ so that
$r_{i}/s_{i}$ is a _potential_ rational reconstruction of $a$. To find _the_
rational reconstruction of $a$, the algorithm is stopped when
$r_{i}<\sqrt{b/2}$. If $s_{i}<\sqrt{b/2}$, then $r_{i}/s_{i}$ is the unique
rational reconstruction of $a$ with numerator and denominator less than
$\sqrt{b/2}$ [38, 39]. Here $b$ corresponds to the product of primes
$p_{1}p_{2}...$ in Algorithm 1.

The Euclidean algorithm can also be used to determine multiplicative inverses
over $\mathbb{Z}_{b}$. Assuming $a$ and $b$ are co-prime, the multiplicative
inverse of $a$ is just $x$. Another option for determining multiplicative
inverses is to use Fermat’s little theorem which states that $a^{p}\equiv
a\mod p$ for $p$ prime. In other words $a^{p-2}$ is the multiplicative inverse
of $a$ assuming that $a$ is not a multiple of $p$. The exact method used to
determine inverses is largely irrelevant since row reducing an $n\times n$
matrix takes $\mathcal{O}(n^{3})$ time overall but only requires computing at
most $n$ inverses.

The final number theoretic result used in FiniteFieldSolve is the Chinese
remainder theorem. Assuming $m_{1}$, $m_{2}$… are pairwise co-prime, the
system of congruences

$\displaystyle x\equiv c_{1}\mod{m_{1}}$ $\displaystyle x\equiv
c_{2}\mod{m_{2}}$ (50) $\displaystyle...$

has a unique solution $x\mod m$ where $m=m_{1}m_{2}$… In order to construct
the solution, define $M_{i}$ to be $m/m_{i}$ and observe that $M_{i}$ is
clearly an integer. Since $M_{i}$ and $m_{i}$ are co-prime, each $M_{i}$ has a
multiplicative inverse $n_{i}$ satisfying $M_{i}n_{i}\equiv 1\mod m_{i}$. The
solution to the linear system of congruences is then

$\displaystyle x=c_{1}n_{1}M_{1}+c_{2}n_{2}M_{2}+...$ (51)

As indicated in Algorithm 1, FiniteFieldSolve only ever uses the Chinese
remainder theorem over two module $m_{1}$ and $m_{2}$. In this case,
$M_{1}=m_{2}$ and $M_{2}=m_{1}$. The Euclidean algorithm then provides $n_{1}$
and $n_{2}$ satisfying $m_{2}n_{1}+m_{1}+n_{2}=1$. The solution $x$ is simply
given by

$\displaystyle x=c_{1}n_{1}m_{2}+c_{2}n_{2}m_{1}.$ (52)

## References

  * [1] S. Laporta, _High precision calculation of multiloop Feynman integrals by difference equations_ , _Int. J. Mod. Phys. A_ 15 (2000) 5087 [hep-ph/0102033].
  * [2] Z. Bern, J.J.M. Carrasco and H. Johansson, _Perturbative Quantum Gravity as a Double Copy of Gauge Theory_ , _Phys. Rev. Lett._ 105 (2010) 061602 [1004.0476].
  * [3] Z. Bern, J.J.M. Carrasco and H. Johansson, _New Relations for Gauge-Theory Amplitudes_ , _Phys. Rev. D_ 78 (2008) 085011 [0805.3993].
  * [4] C. Cheung, K. Kampf, J. Novotny and J. Trnka, _Effective Field Theories from Soft Limits of Scattering Amplitudes_ , _Phys. Rev. Lett._ 114 (2015) 221602 [1412.4095].
  * [5] Z. Bern, J.J. Carrasco, W.-M. Chen, A. Edison, H. Johansson, J. Parra-Martinez et al., _Ultraviolet Properties of $\mathcal{N}=8$ Supergravity at Five Loops_, _Phys. Rev. D_ 98 (2018) 086021 [1804.09311].
  * [6] Z. Bern, J.J.M. Carrasco, W.-M. Chen, H. Johansson, R. Roiban and M. Zeng, _Five-loop four-point integrand of $N=8$ supergravity as a generalized double copy_, _Phys. Rev. D_ 96 (2017) 126012 [1708.06807].
  * [7] Z. Bern, S. Davies and T. Dennen, _Enhanced ultraviolet cancellations in $\mathcal{N}=5$ supergravity at four loops_, _Phys. Rev. D_ 90 (2014) 105011 [1409.3089].
  * [8] Z. Bern, S. Davies, T. Dennen, A.V. Smirnov and V.A. Smirnov, _Ultraviolet Properties of N=4 Supergravity at Four Loops_ , _Phys. Rev. Lett._ 111 (2013) 231302 [1309.2498].
  * [9] Z. Bern, S. Davies, T. Dennen and Y.-t. Huang, _Ultraviolet Cancellations in Half-Maximal Supergravity as a Consequence of the Double-Copy Structure_ , _Phys. Rev. D_ 86 (2012) 105014 [1209.2472].
  * [10] Z. Bern, S. Davies, T. Dennen and Y.-t. Huang, _Absence of Three-Loop Four-Point Divergences in N=4 Supergravity_ , _Phys. Rev. Lett._ 108 (2012) 201301 [1202.3423].
  * [11] Z. Bern, J.J.M. Carrasco, L.J. Dixon, H. Johansson and R. Roiban, _Simplifying Multiloop Integrands and Ultraviolet Divergences of Gauge Theory and Gravity Amplitudes_ , _Phys. Rev. D_ 85 (2012) 105014 [1201.5366].
  * [12] Z. Bern, J.J. Carrasco, L.J. Dixon, H. Johansson and R. Roiban, _Amplitudes and Ultraviolet Behavior of N = 8 Supergravity_ , _Fortsch. Phys._ 59 (2011) 561 [1103.1848].
  * [13] C. Cheung, I.Z. Rothstein and M.P. Solon, _From Scattering Amplitudes to Classical Potentials in the Post-Minkowskian Expansion_ , _Phys. Rev. Lett._ 121 (2018) 251101 [1808.02489].
  * [14] Z. Bern, C. Cheung, R. Roiban, C.-H. Shen, M.P. Solon and M. Zeng, _Scattering Amplitudes and the Conservative Hamiltonian for Binary Systems at Third Post-Minkowskian Order_ , _Phys. Rev. Lett._ 122 (2019) 201603 [1901.04424].
  * [15] Z. Bern, C. Cheung, R. Roiban, C.-H. Shen, M.P. Solon and M. Zeng, _Black Hole Binary Dynamics from the Double Copy and Effective Theory_ , _JHEP_ 10 (2019) 206 [1908.01493].
  * [16] Z. Bern, J. Parra-Martinez, R. Roiban, M.S. Ruf, C.-H. Shen, M.P. Solon et al., _Scattering Amplitudes and Conservative Binary Dynamics at ${\cal O}(G^{4})$_, _Phys. Rev. Lett._ 126 (2021) 171601 [2101.07254].
  * [17] https://gitlab.com/libeigen/eigen.
  * [18] https://github.com/flintlib/flint.
  * [19] T.L. group, _LinBox_ , v1.6.3 ed., 2019.
  * [20] C. Bouillaguet, _SpaSM: a Sparse direct Solver Modulo $p$_, v1.3 ed., 2023\.
  * [21] https://gitlab.com/kaelingre/spasmlink.
  * [22] https://www3.risc.jku.at/education/courses/ss2012/ca2/LinearSystemSolver.m.
  * [23] A.V. Smirnov and F.S. Chuharev, _FIRE6: Feynman Integral REduction with Modular Arithmetic_ , _Comput. Phys. Commun._ 247 (2020) 106877 [1901.07808].
  * [24] A. von Manteuffel and C. Studerus, _Reduze 2 - Distributed Feynman Integral Reduction_ , 1201.4330.
  * [25] J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch, _Integral reduction with Kira 2.0 and finite field methods_ , _Comput. Phys. Commun._ 266 (2021) 108024 [2008.06494].
  * [26] M. Kauers, _Fast solvers for dense linear systems_ , _Nucl. Phys. B Proc. Suppl._ 183 (2008) 245.
  * [27] T. Peraro, _FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs_ , _JHEP_ 07 (2019) 031 [1905.08019].
  * [28] A. Edison, J. Mangan and N.H. Pavao, _Revealing the Landscape of Globally Color-Dual Multi-loop Integrands_ , 2309.16558.
  * [29] A. von Manteuffel and R.M. Schabinger, _A novel approach to integration by parts reduction_ , _Phys. Lett. B_ 744 (2015) 101 [1406.4513].
  * [30] H. Cohen, _A course in computational algebraic number theory_ , Springer (2000), 10.1007/978-3-662-02945-9.
  * [31] G.H. Hardy and E.M. Wright, _An introduction to the theory of numbers_ , Oxford University Press (2008).
  * [32] T. Peraro, _Scattering amplitudes over finite fields and multivariate functional reconstruction_ , _JHEP_ 12 (2016) 030 [1608.01902].
  * [33] V. Magerya, _Rational Tracer: a Tool for Faster Rational Function Reconstruction_ , 2211.03572.
  * [34] J. Klappert and F. Lange, _Reconstructing rational functions with FireFly_ , _Comput. Phys. Commun._ 247 (2020) 106951 [1904.00009].
  * [35] R. Zippel, _Probabilistic algorithms for sparse polynomials_ , in _Symbolic and Algebraic Computation_ , E.W. Ng, ed., (Berlin, Heidelberg), pp. 216–226, Springer Berlin Heidelberg, 1979, DOI.
  * [36] R.A. Demillo and R.J. Lipton, _A probabilistic remark on algebraic program testing_ , _Information Processing Letters_ 7 (1978) 193.
  * [37] J.T. Schwartz, _Fast probabilistic algorithms for verification of polynomial identities_ , _J. ACM_ 27 (1980) 701.
  * [38] P.S. Wang, _A p-adic algorithm for univariate partial fractions_ , in _Proceedings of the Fourth ACM Symposium on Symbolic and Algebraic Computation_ , SYMSAC ’81, (New York, NY, USA), pp. 212–217, Association for Computing Machinery, 1981, DOI.
  * [39] P.S. Wang, M.J.T. Guy and J.H. Davenport, _P-adic reconstruction of rational numbers_ , _SIGSAM Bull._ 16 (1982) 2.
  * [40] M. Monagan, _Maximal quotient rational reconstruction: An almost optimal algorithm for rational reconstruction_ , in _Proceedings of the 2004 International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’04, (New York, NY, USA), pp. 243–249, Association for Computing Machinery, 2004, DOI.
  * [41] T. Granlund and P.L. Montgomery, _Division by invariant integers using multiplication_ , _SIGPLAN Not._ 29 (1994) 61.
  * [42] D. Lemire, O. Kaser and N. Kurz, _Faster Remainder by Direct Computation: Applications to Compilers and Software Libraries_ , _arXiv e-prints_ (2019) arXiv:1902.01961 [1902.01961].
  * [43] https://libdivide.com/.
  * [44] https://learn.microsoft.com/windows/wsl.
  * [45] https://www.cygwin.com.
  * [46] https://www.mingw-w64.org.
  * [47] C. Cheung, J. Mangan and C.-H. Shen, _Hidden Conformal Invariance of Scalar Effective Field Theories_ , _Phys. Rev. D_ 102 (2020) 125009 [2005.13027].
  * [48] Z. Bern, J.J. Carrasco, M. Chiodaroli, H. Johansson and R. Roiban, _The Duality Between Color and Kinematics and its Applications_ , 1909.01358.
  * [49] Z. Bern, J.J. Carrasco, M. Chiodaroli, H. Johansson and R. Roiban, _The SAGEX review on scattering amplitudes Chapter 2: An invitation to color-kinematics duality and the double copy_ , _J. Phys. A_ 55 (2022) 443003 [2203.13013].
  * [50] T. Adamo, J.J.M. Carrasco, M. Carrillo-González, M. Chiodaroli, H. Elvang, H. Johansson et al., _Snowmass White Paper: the Double Copy and its Applications_ , in _Snowmass 2021_ , 4, 2022 [2204.06547].
  * [51] R. Monteiro and D. O’Connell, _The Kinematic Algebra From the Self-Dual Sector_ , _JHEP_ 07 (2011) 007 [1105.2565].
  * [52] M. Ben-Shahar and H. Johansson, _Off-shell color-kinematics duality for Chern-Simons_ , _JHEP_ 08 (2022) 035 [2112.11452].
  * [53] C. Cheung and J. Mangan, _Scattering Amplitudes and the Navier-Stokes Equation_ , 2010.15970.
  * [54] C. Cheung and J. Mangan, _Covariant color-kinematics duality_ , _JHEP_ 11 (2021) 069 [2108.02276].
  * [55] C. Cheung, J. Mangan, J. Parra-Martinez and N. Shah, _Non-perturbative Double Copy in Flatland_ , _Phys. Rev. Lett._ 129 (2022) 221602 [2204.07130].
  * [56] C. Cheung and C.-H. Shen, _Symmetry for Flavor-Kinematics Duality from an Action_ , _Phys. Rev. Lett._ 118 (2017) 121601 [1612.00868].
  * [57] M. Ben-Shahar and M. Guillen, _10D super-Yang-Mills scattering amplitudes from its pure spinor action_ , _JHEP_ 12 (2021) 014 [2108.11708].
  * [58] M. Ben-Shahar, L. Garozzo and H. Johansson, _Lagrangians manifesting color-kinematics duality in the NMHV sector of Yang-Mills_ , _JHEP_ 08 (2023) 222 [2301.00233].
  * [59] J.J.M. Carrasco and L. Rodina, _UV considerations on scattering amplitudes in a web of theories_ , _Phys. Rev. D_ 100 (2019) 125007 [1908.08033].
  * [60] A. Brandhuber, G. Chen, H. Johansson, G. Travaglini and C. Wen, _Kinematic Hopf Algebra for Bern-Carrasco-Johansson Numerators in Heavy-Mass Effective Field Theory and Yang-Mills Theory_ , _Phys. Rev. Lett._ 128 (2022) 121601 [2111.15649].
  * [61] Z. Bern, S. Davies and J. Nohle, _Double-Copy Constructions and Unitarity Cuts_ , _Phys. Rev. D_ 93 (2016) 105015 [1510.03448].
  * [62] G. Mogull and D. O’Connell, _Overcoming Obstacles to Colour-Kinematics Duality at Two Loops_ , _JHEP_ 12 (2015) 135 [1511.06652].
  * [63] Z. Bern, L.J. Dixon, D.C. Dunbar and D.A. Kosower, _Fusing gauge theory tree amplitudes into loop amplitudes_ , _Nucl. Phys. B_ 435 (1995) 59 [hep-ph/9409265].
  * [64] R. Britto, F. Cachazo and B. Feng, _Generalized unitarity and one-loop amplitudes in N=4 super-Yang-Mills_ , _Nucl. Phys. B_ 725 (2005) 275 [hep-th/0412103].
  * [65] Z. Bern and Y.-t. Huang, _Basics of Generalized Unitarity_ , _J. Phys. A_ 44 (2011) 454003 [1103.1869].
  * [66] V.E. Zakharov and A.V. Mikhailov, _Relativistically Invariant Two-Dimensional Models in Field Theory Integrable by the Inverse Problem Technique. (In Russian)_ , _Sov. Phys. JETP_ 47 (1978) 1017\.
  * [67] J.J.M. Carrasco, M. Lewandowski and N.H. Pavao, _Color-Dual Fates of F3, R3, and N=4 Supergravity_ , _Phys. Rev. Lett._ 131 (2023) 051601 [2203.03592].
  * [68] J.J.M. Carrasco, M. Lewandowski and N.H. Pavao, _Double-copy towards supergravity inflation with $\alpha$-attractor models_, _JHEP_ 02 (2023) 015 [2211.04441].
  * [69] J.J.M. Carrasco and N.H. Pavao, _UV Massive Resonance from IR Double Copy Consistency_ , 2310.06316.
  * [70] A.S.-K. Chen, H. Elvang and A. Herderschee, _Bootstrapping the String Kawai-Lewellen-Tye Kernel_ , _Phys. Rev. Lett._ 131 (2023) 031602 [2302.04895].
  * [71] H.-H. Chi, H. Elvang, A. Herderschee, C.R.T. Jones and S. Paranjape, _Generalizations of the double-copy: the KLT bootstrap_ , _JHEP_ 03 (2022) 077 [2106.12600].
  * [72] T.V. Brown, K. Kampf, U. Oktem, S. Paranjape and J. Trnka, _Scalar BCJ Bootstrap_ , 2305.05688.
  * [73] Y. Li, D. Roest and T. ter Veldhuis, _Hybrid Goldstone Modes from the Double Copy Bootstrap_ , 2307.13418.
  * [74] J.J.M. Carrasco, L. Rodina and S. Zekioglu, _Composing effective prediction at five points_ , _JHEP_ 06 (2021) 169 [2104.08370].
  * [75] J.J.M. Carrasco, L. Rodina, Z. Yin and S. Zekioglu, _Simple encoding of higher derivative gauge and gravity counterterms_ , _Phys. Rev. Lett._ 125 (2020) 251602 [1910.12850].
  * [76] J. Broedel and L.J. Dixon, _Color-kinematics duality and double-copy construction for amplitudes from higher-dimension operators_ , _JHEP_ 10 (2012) 091 [1208.0876].
  * [77] Q. Bonnefoy, G. Durieux, C. Grojean, C.S. Machado and J. Roosmale Nepveu, _The seeds of EFT double copy_ , _JHEP_ 05 (2022) 042 [2112.11453].
  * [78] J.J.M. Carrasco, C.R. Mafra and O. Schlotterer, _Semi-abelian Z-theory: NLSM $+\phi^{3}$ from the open string_, _JHEP_ 08 (2017) 135 [1612.06446].
  * [79] J.J.M. Carrasco, C.R. Mafra and O. Schlotterer, _Abelian Z-theory: NLSM amplitudes and $\alpha$’-corrections from the open string_, _JHEP_ 06 (2017) 093 [1608.02569].
  * [80] I. Low, L. Rodina and Z. Yin, _Double Copy in Higher Derivative Operators of Nambu-Goldstone Bosons_ , _Phys. Rev. D_ 103 (2021) 025004 [2009.00008].
  * [81] I. Low and Z. Yin, _New flavor-kinematics dualities and extensions of nonlinear sigma models_ , _Phys. Lett. B_ 807 (2020) 135544 [1911.08490].
  * [82] L.F. Alday, V. Gonçalves, M. Nocchi and X. Zhou, _Six-Point AdS Gluon Amplitudes from Flat Space and Factorization_ , 2307.06884.
  * [83] L.F. Alday and T. Hansen, _The AdS Virasoro-Shapiro amplitude_ , _JHEP_ 10 (2023) 023 [2306.12786].
  * [84] L.F. Alday, T. Hansen and J.A. Silva, _Emergent Worldsheet for the AdS Virasoro-Shapiro Amplitude_ , _Phys. Rev. Lett._ 131 (2023) 161603 [2305.03593].
  * [85] L.F. Alday, T. Hansen and J.A. Silva, _AdS Virasoro-Shapiro from single-valued periods_ , _JHEP_ 12 (2022) 010 [2209.06223].
  * [86] L.F. Alday, T. Hansen and J.A. Silva, _AdS Virasoro-Shapiro from dispersive sum rules_ , _JHEP_ 10 (2022) 036 [2204.07542].
  * [87] L.F. Alday, C. Behan, P. Ferrero and X. Zhou, _Gluon Scattering in AdS from CFT_ , _JHEP_ 06 (2021) 020 [2103.15830].
  * [88] C. Cheung, K. Kampf, J. Novotny, C.-H. Shen and J. Trnka, _A Periodic Table of Effective Field Theories_ , _JHEP_ 02 (2017) 020 [1611.03137].
  * [89] A. Nützi and M. Reiterer, _Scattering amplitude annihilators_ , _JHEP_ 02 (2020) 020 [1905.02224].
  * [90] F. Loebbert, M. Mojaza and J. Plefka, _Hidden Conformal Symmetry in Tree-Level Graviton Scattering_ , _JHEP_ 05 (2018) 208 [1802.05999].
  * [91] C. Cheung, K. Kampf, J. Novotny, C.-H. Shen, J. Trnka and C. Wen, _Vector Effective Field Theories from Soft Limits_ , _Phys. Rev. Lett._ 120 (2018) 261602 [1801.01496].
  * [92] L. Ren, M. Spradlin, A. Yelleshpur Srikant and A. Volovich, _On effective field theories with celestial duals_ , _JHEP_ 08 (2022) 251 [2206.08322].
  * [93] L. Rodina, _Scattering Amplitudes from Soft Theorems and Infrared Behavior_ , _Phys. Rev. Lett._ 122 (2019) 071601 [1807.09738].
  * [94] N. Arkani-Hamed, L. Rodina and J. Trnka, _Locality and Unitarity of Scattering Amplitudes from Singularities and Gauge Invariance_ , _Phys. Rev. Lett._ 120 (2018) 231602 [1612.02797].
  * [95] I. Low and Z. Yin, _Soft Bootstrap and Effective Field Theories_ , _JHEP_ 11 (2019) 078 [1904.12859].
  * [96] I. Low, J. Shu, M.-L. Xiao and Y.-H. Zheng, _Amplitude/operator basis in chiral perturbation theory_ , _JHEP_ 01 (2023) 031 [2209.00198].
  * [97] J. Ellis, _TikZ-Feynman: Feynman diagrams with TikZ_ , _Comput. Phys. Commun._ 210 (2017) 103 [1601.05437].