# Numerical checks

This file records the numerical comparisons behind the statements of accuracy in the paper. Formulas are written in LaTeX notation, and equation labels refer to the paper source.

## Tests of the results.

 We collect here the numerical tests that the main text refers to. To test $\Lzero$ without numerical differentiation we integrate $\Lzero\Vc$ against polynomial test functions that vanish at the ends of the interval. The relative residuals are below $10^{-18}$, which is the level set by the accuracy of the values. With one coefficient of $\Lzero$ changed by one part in $10^{6}$ they rise to $10^{-11}$, and the operator of order 14 gives $10^{-7}$. The operator of order 9 from the second reduction annihilates the values to $10^{-27}$.

On the line~\eqref{eq:line}, eq.~\eqref{eq:closedform} differs from the numerical values by a relative amount with median $5\times10^{-34}$ and maximum $2\times10^{-30}$. It agrees with eight evaluations at higher precision to between 42 and 54 digits. With one rational coefficient of the formula changed, the agreement drops to $10^{-3}$. For general triangles with one leg per site the 14 coefficients were fitted at 45 points. The formula agrees with the integral to better than 34 digits at 21 further generic points, and to between 30 and 34 digits at two nearly degenerate triangles and two soft configurations.

On the second line the representation has 44 integration constants. With these fitted to half of 128 values it reproduces the other half to 38 digits, and a generic space of functions of the same dimension reaches only $10^{-32}$. In general kinematics the count of one symbol for each square root was made modulo two primes. The ten rational numbers were fitted at 140 generic points and agree with the values of eq.~\eqref{eq:closedformgeneral} to $10^{-41}$. With these values the formula reproduces 809 further values with $X_v\ge P_v$ to 43 digits or better. A change of one of the ten numbers by one part in $10^{3}$ changes the agreement to a few per cent. The continuation to $X_v<P_v$ agrees with the integral to 50 digits at 372 further points, 328 of which have at least one $X_v<P_v$. For the de~Sitter coefficient of section~\ref{sec:discussion}, the difference of the finite parts at two kinematic points agrees with the logarithm of eq.~\eqref{eq:dsform} to 29 digits. The loop-integrated correlator equals $-3\pi^2\zeta(3)$ to 40 digits at eight kinematic points.

For the elliptic sector of section~\ref{sec:elliptic}, the independent verification covers the differential equation, the change of basis for the elliptic block, the poles, and a sample of the constants. It also covers the boundary values at $\lambda=\tfrac12$, to more than 30 digits.

## Tests for the box.

 For the soft limit~\eqref{eq:softlimit}, adding $D_0\log 2P_1$ to $V_4$ lowers the tail of its Chebyshev expansion in $P_1$ from $10^{-10}$ to $10^{-25}$. A coefficient of the logarithm that is wrong by one part in $10^{6}$ is detected. Equation~\eqref{eq:D0closed} agrees with a numerical integration of eq.~\eqref{eq:D0} to 33 digits or better at 45 shapes. A second computation, which reduces eq.~\eqref{eq:D0} to the periods of the curve by a different route, reproduces eq.~\eqref{eq:D0closed} at one shape. The exact test of section~\ref{sec:boxmonodromy} was made modulo three primes at five shapes of the first family.

The numerical values of $V_4$ and $C_4$ on the line of section~\ref{sec:boxline} were computed to 32 digits at more than 170 points with $0\le t\le100$. For the correlator we computed 120 further values at larger site energies, 40 with $7\le X\le15.5$ and 80 with $16\le X\le5000$. An unconstrained fit of the 67 or 81 boundary constants on the interval $0.1\le t\le2$ reproduces the values to 33 digits, and it does the same for a different function. With the additional data and the conditions described in section~\ref{sec:boxline}, the fitted solutions predict the remaining values to 33 digits. The same fit to the wrong function fails at the level of $10^{-10}$ for the correlator and $10^{-15}$ for the coefficient. A function that differs from the correct one by the factor $1+10^{-9}\,t$ is rejected, and so is an error of $10^{-6}$ in one entry of the system. The systems obtained by setting $\eps=0$ first, or by keeping simple poles only, fail at the level of $10^{-11}$ and $10^{-15}$. At $t=0$ the solutions agree with a direct evaluation at one leg per site to 30 digits.

Two independent quadratures of $\Omega$ agree to 30 digits. The search for an integer relation used powers of $\pi$, values of the gamma function at arguments $k/8$ and $\tfrac13$, $\log2$, $\log3$, Catalan's constant, $\mathrm{Cl}_2(\pi/3)$, $\arccos\tfrac13$ and $\arctan\sqrt2$, with integer coefficients up to 300 at 28 digits. With the value of $\Omega$ in eq.~\eqref{eq:boxclosed}, the expression $u+\Omega\,v$ reproduces 173 numerical values with $t>0$ to at least 32.3 digits, and the 120 further values to 33 digits or better. For the 80 values with $X\ge16$ the expression was evaluated from the series~\eqref{eq:boxlargeX}, and for the 40 values with $7\le X\le15.5$ from the differential system. The value of $\Omega$ determined from these values alone agrees with the quadrature to $10^{-34}$, and replacing $\Omega$ by $\Omega\,(1+10^{-25})$ lowers the agreement to 23 digits. The solution fitted to the numerical values, which makes no use of the expansion by regions, has $\alpha_4$ equal to $\Omega$ to 30 digits and $\beta_5$, $\alpha_5$ and $\beta_7$ equal to their exact values to 33 digits.

# Method details

These paragraphs describe how the numerical values and the exact operators were obtained. They were part of an earlier version of appendix A.

## Numerical evaluation

The integral over $y_{12}$ in eq.~\eqref{eq:V} can be done in closed form, which leaves a two-dimensional integral over a region bounded by the triangle inequalities. We evaluate it with nested double-exponential quadrature in ball arithmetic, which controls rounding errors rigorously. We estimate the quadrature error from successive refinements of the step. On the line~\eqref{eq:line} we computed $\Vc$ at 288 Chebyshev nodes of the interval $0.345\le\lambda\le0.985$ to 33 to 35 digits. In general kinematics we evaluated $\Vc$ at 1382 points, with errors estimated in the same way. An independent three-dimensional quadrature of the loop integral in double precision agrees with the value at $\lambda=\tfrac7{10}$ to 15 digits, with eq.~\eqref{eq:closedform} to 14 digits at four triangles off the line, and with eq.~\eqref{eq:closedformgeneral} to 14 digits at five points with $X_v>P_v$.

## Exact reduction

The integration-by-parts identities of the family~\eqref{eq:family} follow from the vector fields that preserve $\Bk$ up to a factor. We solve them by linear algebra modulo prime numbers at numerical values of $\lambda$ and $\eps$. The class of the integrand of $\Vc$ and its $\lambda$-derivatives span a finite-dimensional space, and the first linear relation among them is the differential operator. Its coefficients are reconstructed as rational functions of $\lambda$ from many sample points. They are lifted to the rational numbers with the Chinese remainder theorem and rational reconstruction, and checked at primes that were not used in the lift.

## Factorization

Right factors of first order are found from the solutions of the operator and of its adjoint whose logarithmic derivative is a rational function. Such solutions are called hyperexponential, and the corresponding factors are divided out exactly. The test for a factor equivalent to $\Ltwo$ is an exact computation. Two operators are equivalent if a differential operator with rational coefficients maps the solutions of one to those of the other. For a bounded degree of these coefficients the existence of such a map is a question about the rank of a linear system. We repeated it with an enlarged bound.

## Fits

The constants of the closed forms~\eqref{eq:closedform} and~\eqref{eq:closedformgeneral} were obtained by a practice that is standard for multiloop amplitudes \cite{Acres:2021sss}. The function space is fixed from the cuts of the integral, and the integral is evaluated to high precision. An integer-relation algorithm then determines the rational coefficients of an ansatz in that space. The result is accepted only if it reproduces digits that were not used in the fit. We applied the same standard to the constant $-3\pi^2\zeta(3)$ of section~\ref{sec:discussion}. No constants were fitted in the elliptic sector, where the boundary values are numerical and the constants of the iterated integrals have closed forms (section~\ref{sec:ellipticresults}).