No-neighbours recurrence schemes for space-time Green's functions on a 3D simple cubic lattice

On the face of it, continuous space and lattice-based (discrete) space seem worlds apart. However, they meet in physics models from either direction. Recent developments in symbolic discrete mathematics weave through these apparently different branches of physics.

On the one hand, in lattice-based physics research, e.g. in condensed matter physics, crystallography, and statistical physics, periodic arrangements of atoms on lattices dictate a discrete approach involving finite-cell-size interactions. The elementary building block in studying such interactions is the lattice Green's function (LGF), $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} P(\nvK;\xi)$, i.e. the field at a lattice point $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK$ due to a unit excitation at the lattice origin $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK=\nv{0}$, where the parameter ξ may be mapped to temporal frequency. Alternatively, upon expanding P in terms of a power series in ξ, the coefficient of $\xi^{n}$ is the probability that a random walker starting at $\newcommand{\nv}[1]{\vec{#1}} \nv{0}$ on the lattice arrives at $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK$ in n steps. It is impossible to list all contributions to the theory of LGFs. As early as 1928 Courant, in a paper on difference equations on a 3D simple cubic lattice, derived an integral representation for the associated LGF [1]. At the lattice origin, this representation reduces to one of the three Watson integrals [2] for the respective simple, face-centred, and body-centred cubic lattices, further developments and generalisations of which have been described in the overview paper by Zucker [3]. The theory of LGFs has progressed in leaps and bounds over the past decade, following the lucid papers by Guttmann on the LGFs in arbitrary spatial dimensions in an enumerative combinatorics setting [4], and its connection with Picard–Fuchs ordinary differential equations of Calabi–Yau manifolds [5]. This has inspired a fundamental mathematical analysis in terms of the minimal linear differential operators that annihilate diagonals of rational functions in the framework of differential Galois groups [6, 7]. It was recognised that the integrands in the LGF integral representations are δ-finite holonomic functions [8], and hence that creative telescoping (or similar techniques) can be used to construct efficient recurrence relations for consecutive coefficients of a power expansion of $\newcommand{\nv}[1]{\vec{#1}} P(\nv{0};\xi)$ in ξ, but also with respect to the lattice dimension $\newcommand{\Gn}{\nu} \newcommand{\G}{G} \Gn$ (d in the papers cited here) for the polynomial coefficients of the linear differential operators that annihilate the LGFs [8–10].

On the other hand, in continuum physics, finite-difference discretisation is a widely accepted compromise in order to solve practical differential equations, and this was what motivated Courant, Friedrichs, and Lewy in their seminal paper [11], to which the LGF integral representation in [1] was an addition omitted from the former due to lack of space. Although the lattice-based and continuous space physics applications differ, the mathematical underpinning is exactly the same. Nevertheless, in the latter case, in the 1960s the interest shifted away from the elliptic case towards explicit schemes for solving hyperbolic initial-value problems through a marching scheme in discrete space-time. In the electromagnetism community this approach [12] has been dubbed the finite-difference time-domain (FDTD) method [13]. We denote the scalar counterpart of the asssociated discrete space-time-domain Green's function as $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} G_{\nvK}^{n}$, and its Z-transform domain counterpart as $\newcommand{\Gn}{\nu} \newcommand{\G}{G} \newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} G_{\nvK}(Z)=P(-\nvK, \xi)\xi/(2\Gn)$ with $\xi=2/(Z+Z^{-1})$.

Rather than the LGF at the lattice origin (or at other special symmetry lattice points) in arbitrary dimensions, we are primarily interested in LGF discrete time sequences for fixed spatial dimensions at arbitrary lattice points. Our motivation stems from almost two decades of work on discrete Green's function based boundary conditions [14–16], and later, on discrete Green's function diakoptics, in which LGF discrete time sequences are used in a self-consistent combination with FDTD to capture the interactions between possibly disjoint computational FDTD domains and their lattice-based environment [17–20]. Crucial to the success of discrete function diakoptics is the possibility to generate long LGF sequences on the fly. In [21] we constructed $\newcommand{\Gn}{\nu} \newcommand{\G}{G} {\Gn}D$ and $\newcommand{\Gn}{\nu} \newcommand{\G}{G} (\Gn-1)D$ nested combinatorial sums for the discrete time LGF sequences for a $\newcommand{\Gn}{\nu} \newcommand{\G}{G} {\Gn}$D simple cubic lattice. These sums are not practical for number crunching. However, for the 2D case those nested sums were amenable to univariate creative telescoping [22–24], resulting in a fast, no-neighbours, third-order, sixth-degree discrete-time recurrence scheme for LGF sequences at arbitrary points on a square lattice.

Here, we employ multivariate creative telescoping [25] to obtain a fast no-neighbours, seventh-order, eighteenth-degree, discrete-time recurrence scheme for arbitrary lattice points on a 3D simple cubic lattice. This result is far more intricate than we had expected, but also much richer in features. The complexity reduces for lattice points on diagonal symmetry planes, resulting in a fast no-neighbours, fifth-order, twelfth-degree, discrete-time recurrence scheme.

In condensed-matter and other lattice-based physics, one is interested in the Brillouin zone dispersion of time-harmonic waves for wavelengths commensurate to the grid size, whereas in FDTD, one would like to sample sufficiently smooth wavefields fine enough in time so that dispersion due to wave-lattice interactions is suppressed. Our Green's function discrete time sequences apply to either field of research. Here we shall restrict ourselves to presenting the resulting no-neighbours recurrence schemes through compact tables of coefficients, and to discussing their ramifications and the features of an LGF discrete time sequence through an example.

As mentioned in the introduction, our interest in discrete Green's-function diakoptics, [26], sparked the research reported here concerning the scalar 3D discrete space-time lattice Green's function (LGF) on a simple cubic lattice with lattice points $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \newcommand{\Kt}{\K_{2}} \newcommand{\Kd}{\K_{3}} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK=(\Ko, \Kt, \Kd)$. This LGF is also referred to as the scalar FDTD Green's function, and satisfies a unit-excitation discrete wave equation.

In this section, we consider the more general case of a simple hypercubic lattice in $\newcommand{\Gn}{\nu} \newcommand{\G}{G} \Gn$ spatial dimensions, for which the discrete wave equation can be written in terms of a second-order explicit scheme

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\nv}[1]{\vec{#1}} \newcommand{\nm}[1]{\mathsf{#1}} \newcommand{\G}{G} \newcommand{\Gn}{\nu} \newcommand{\Gd}{\delta} \displaystyle \label{nuD-GF-eq} G^{n+1} = -G^{n-1}+\frac{1}{\Gn}\left(\nm{a}G^{n}+\Gd_{\nv{0}}\Gd_{0}\right), \nonumber \end{align} \tag{ 1 }$$

in which n is a discrete time index, and $\newcommand{\nm}[1]{\mathsf{#1}} \nm{a}$ denotes the (unit) adjacency matrix. The source term $\newcommand{\Gd}{\delta} \newcommand{\G}{G} \newcommand{\nv}[1]{\vec{#1}} \Gd_{\nv{0}}\Gd_{0}$ represents the action of a discrete space-time point source. It is 1 at the space-time vertex $\newcommand{\nv}[1]{\vec{#1}} (\nv{0}, 0)$, and 0 otherwise. We have adjusted the time-step $\newcommand{\GD}{\Delta} \newcommand{\G}{G} \GD{t}$ to the unit grid size h  =  1, so that for a constant wavespeed c, the Courant number $\newcommand{\GD}{\Delta} \newcommand{\G}{G} \alpha=\sqrt{\nu}\, c\GD{t}$ is at its maximum value, $\alpha=1$, for the explicit scheme to (just) satisfy the Courant–Friedrichs–Lewy (CFL) condition for stability [27].

Owing to the simple hypercubic lattice symmetry, permutations of, and parity inversions in any of the lattice-point spatial coordinates leave the scalar LGF invariant. Hence, we may restrict our analysis to that part of the lattice for which the lattice points $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK$ are non-negative, i.e. $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK$ has non-negative coordinates⁴. The LGF vanishes prior to the lattice arrival time

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ol}{u} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \newcommand{\G}{G} \newcommand{\Gn}{\nu} \displaystyle \label{lattice-arr-time} n_{O}=1+\KM=1+\sum_{\ol=1}^{\Gn}\ell_{\ol} \nonumber \end{align} \tag{ 2 }$$

in the Manhattan metric. For our choice of the Courant parameter, $\alpha=1$, every other term in the LGF sequence associated with an arbitrary lattice point $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK$ vanishes.

Starting from an integral representation for the LGF in the $\newcommand{\Z}{Z} \Z$-domain, one may readily derive a closed form expression for the $\newcommand{\Gn}{\nu} \newcommand{\G}{G} \Gn$-dimensional space-time LGF in terms of $\newcommand{\Gn}{\nu} \newcommand{\G}{G} \Gn$ nested finite sums [21], i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\nv}[1]{\vec{#1}} \newcommand{\nvK}{\nvk} \newcommand{\nvk}{\nv{\K}} \newcommand{\fracD}[2]{\frac{#1}{#2}} \newcommand{\ol}{u} \newcommand{\Ko}{\K_{1}} \newcommand{\KGn}{\K_{\Gn}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \newcommand{\G}{G} \newcommand{\Gn}{\nu} \displaystyle \label{G_K^n} G_{\nvK}^{1+\K+2n} = \sum\limits_{\sum_{\ol=1}^{\Gn}m_{\ol}=m\leqslant n} \fracD{(-1)^{n-m}\Gn^{-1-\KM-2m}(\KM+n+m)!}{(n-m)!m_{1}!\cdots m_{\Gn}!(\Ko+m_{1})!\cdots(\KGn+m_{\Gn})!}, \nonumber \end{align} \tag{ 3 }$$

where n is assumed non-negative, and $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} G_{\nvK}^{\K+2n}=0$. We have also shown that upon applying the Chu–Vandermonde identity [28] for $\newcommand{\Gn}{\nu} \newcommand{\G}{G} \Gn\geqslant3$, the ν nested finite sums, can be reduced to a representation in terms of $(\nu-1)$ nested finite sums, given by [21, equations (13) or (14)]. Unfortunately, two typos appear in [21, equation (14)], and in the definition of the multinomial in [21, equation (11)]. The correct expressions are

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\GSl}{\lambda} \newcommand{\GSm}{\mu} \newcommand{\nv}[1]{\vec{#1}} \newcommand{\nvK}{\nvk} \newcommand{\nvk}{\nv{\K}} \newcommand{\Ko}{\K_{1}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \newcommand{\G}{G} \newcommand{\Gn}{\nu} \newcommand{\GS}{\Sigma} \displaystyle G_{\nvK}^{1+\KM+2n} =& \, \sum\limits_{m=0}^{n}(-1)^{n-m}\Gn^{-1-\KM-2m}\left(\begin{array}{@{}c@{}}{\KM+n+m} \nonumber \\{n}\end{array}\right)\left(\begin{array}{@{}c@{}}{n} \nonumber \\{m}\end{array}\right)\nonumber \\ & \mbox{}\times \sum\limits_{m_{1}=0}^{m-\GSm_{0}}\cdots\sum\limits_{m_{\Gn-2}=0}^{m-\GSm_{\Gn-3}} \left(\begin{array}{@{}c@{}}{\GSm_{\Gn-2}} \nonumber \\{m_{1},\ldots,m_{\Gn-2}}\end{array}\right) \left(\begin{array}{@{}c@{}}{\GSl_{\Gn-2}+\GSm_{\Gn-2}} \nonumber \\{\Ko+m_{1},\ldots,\K_{\Gn-2}+m_{\Gn-2}}\end{array}\right)\nonumber \\ & \mbox{}\times \left(\begin{array}{@{}c@{}}{\KM+2m-(\GSl_{\Gn-2}+2\GSm_{\Gn-2})} \nonumber \\{\KM+m-(\GSl_{\Gn-1}+\GSm_{\Gn-2})}\end{array}\right) \left(\begin{array}{@{}c@{}}{\KM+m} \nonumber \\{\GSl_{\Gn-2}+\GSm_{\Gn-2}}\end{array}\right)\left(\begin{array}{@{}c@{}}{m} \nonumber \\{\GSm_{\Gn-2}}\end{array}\right),\label{G_K^n_nomials} \nonumber \end{align} \tag{ 4 }$$

in which $\newcommand{\GS}{\Sigma} \newcommand{\G}{G} \newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\KM}{\K} \newcommand{\GSl}{\lambda} \GSl_{i}=\sum_{j=1}^{i}\KM_{j}$ and $\newcommand{\GS}{\Sigma} \newcommand{\G}{G} \newcommand{\GSm}{\mu} \GSm_{i}=\sum_{j=1}^{i}m_{j}$ with $\newcommand{\GS}{\Sigma} \newcommand{\G}{G} \newcommand{\GSl}{\lambda} \GSl_{0}=0$ and $\newcommand{\GS}{\Sigma} \newcommand{\G}{G} \newcommand{\GSm}{\mu} \GSm_{0}=0$, and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\fracD}[2]{\frac{#1}{#2}} \newcommand{\G}{G} \newcommand{\Gn}{\nu} \displaystyle \left(\begin{array}{@{}c@{}}{m} \nonumber \\{m_{1},\ldots,m_{\Gn}}\end{array}\right) = \fracD{m!}{m_{1}!\cdots{m_{\Gn}}!}, \nonumber \end{align} \tag{ 5 }$$

respectively⁵.

Note that (4) is much more efficient than (3), in part due to the O(n²) and O(n³) number of terms in the respective double and triple nested sums, but also due to multinomials being more efficient in memory and number of operations than its factorial counterparts. However, neither of these finite nested sum representations are very practical except for quite modest grid sizes and short LGF sequences computations, because in addition to the growing number of terms in the nested sums, the terms themselves at every new time step grow very fast in magnitude, whereas the corresponding LGF actually decays owing to severe cancellations.

Nevertheless, both (3) and (4) are relevant for the work reported below, in that (3) is instrumental in obtaining recurrence schemes for LGFs, while (4) is most efficient for the determination of the initial conditions that start off an LGF recurrence scheme in case more than one initial condition is required. Below, we briefly review the case for $\newcommand{\Gn}{\nu} \newcommand{\G}{G} \Gn=2$, before proceeding with the recurrence schemes for non-negative lattice points on diagonal symmetry planes in $\mathbb{Z}^{3}$, and finally, for arbitrary points in $\mathbb{Z}^{3}_{\geqslant0}$.

For $\nu=2$, the original 2D nested finite sum reduced to a single sum, to which in [21] we applied the procedure of creative telescoping [22, 23] to arrive at a recurrence operator (or telescoper) that acts on the LGF $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} G(n)=G_{\nvK}^{1+\K+2n}$, of the form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\oL}{\Lambda} \newcommand{\ol}{u} \displaystyle \label{RecOp} L = \sum\limits_{\ol=0}^{\oL}p_{\ol}(n)N^{\ol}, \nonumber \end{align} \tag{ 6 }$$

in which $\newcommand{\oL}{\Lambda} \oL$ denotes the order of the recurrence, while N is the shift operator, i.e. $NG(n)=G(n+1)$. The recurrence operator annihilates the Green's function through the recurrence relation

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\oL}{\Lambda} \newcommand{\fracD}[2]{\frac{#1}{#2}} \newcommand{\ol}{u} \displaystyle \label{RecRel} LG = \sum\limits_{\ol=0}^{\oL}p_{\ol}(n)G(n+\ol)=0 \quad \!\!\Rightarrow\quad \!\! G(n+\oL)=-\sum\limits_{\ol=0}^{\oL-1}\fracD{p_{\ol}(n)}{p_{\oL}(n)}G(n+\ol), \nonumber \end{align} \tag{ 7 }$$

where $\newcommand{\ol}{u} p_{\ol}(n)$ are polynomial coefficients of degree d in n. For $\nu=2$, the Mathematica package Zeil by Petkovsěk [24] applied to the single sum

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\nv}[1]{\vec{#1}} \newcommand{\nvK}{\nvk} \newcommand{\nvk}{\nv{\K}} \newcommand{\Ko}{\K_{1}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle \label{G_K^n_2D} G_{\nvK}^{1+\KM+2n} = \sum\limits_{m=0}^{n}(-1)^{n-m}2^{-1-\KM-2m} \left(\begin{array}{@{}c@{}}{\KM+n+m} \nonumber \\{n}\end{array}\right)\left(\begin{array}{@{}c@{}}{n} \nonumber \\{m}\end{array}\right)\left(\begin{array}{@{}c@{}}{\KM+2m} \nonumber \\{\Ko+m}\end{array}\right) \nonumber \end{align} \tag{ 8 }$$

with generic non-negative values of $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \Ko$, $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kt}{\K_{2}} \Kt$ (and $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\KM}{\K} \newcommand{\Ko}{\K_{1}} \newcommand{\Kt}{\K_{2}} \KM=\Ko+\Kt$) leads to a generic recurrence relation of order $\newcommand{\oL}{\Lambda} \oL=3$ with polynomial coefficients of degree d  =  6 [21].

Interestingly, for $\newcommand{\ol}{u} \ol=1, 2$, the polynomial coefficients $\newcommand{\ol}{u} p_{\ol}(n)$ of degree 6 display a certain structure, in that they may be factorised in term of polynomials that are linear in n, $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \Ko$ and $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kt}{\K_{2}} \Kt$, and polynomials of degree 4 in n, $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \Ko$ and $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kt}{\K_{2}} \Kt$ that are irreducible over the field $\mathbb{Q}$. For $\newcommand{\ol}{u} \ol=0, 3$, the polynomial coefficients $\newcommand{\ol}{u} p_{\ol}(n)$ of degree 6 factorise into six linear polynomials.

For lattice points on the diagonal ($\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \newcommand{\Kt}{\K_{2}} \Kt=\Ko$), Zeil produces a recurrence of order $\newcommand{\oL}{\Lambda} \oL=2$ and degree d  =  2, which is clearly less complex than the generic case. From a range of numerical experiments for lattice points for which $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \newcommand{\Kt}{\K_{2}} \newcommand{\cfc}{s} \vert \Kt-\Ko\vert =\cfc$ with $\newcommand{\cfc}{s} \cfc$ a numeric rather than a symbolic value, we found that the order remained $\newcommand{\oL}{\Lambda} \oL=2$, while the degree grows as $\newcommand{\cfc}{s} 2\lfloor\cfc/2\rfloor+2+2\min(\cfc, 1)$, with $\lfloor\ldots\rfloor$ denoting the greatest lower bound integer. Hence, although the order may be reduced by considering lattice points at a fixed numeric distance to the main diagonal, the degree grows rapidly away from that diagonal and soon the computational complexity exceeds that of the generic third-order case. Moreover, the lattice-point dependent degree precludes any straightforward construction of the generic recurrence operator from the specific off-diagonal cases. Assigning both $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \Ko$ and $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kt}{\K_{2}} \Kt$ (and hence $\newcommand{\cfc}{s} \cfc$) a numeric rather than a symbolic value does not reduce the order or degree either.

Now, let us address the elephant in the room, viz., we must investigate whether it is possible that $\newcommand{\oL}{\Lambda} p_{\oL}(n)=0$ for some $n\geqslant0$, which would lead to breakdown of the recurrence scheme in (7). For the generic recurrence relation ($\newcommand{\oL}{\Lambda} \oL=3$) for the 2D case with $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\KM}{\K} \newcommand{\Ko}{\K_{1}} \newcommand{\Kt}{\K_{2}} \KM=\Ko+\Kt$, we have (see, [21])

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Kt}{\K_{2}} \newcommand{\Ko}{\K_{1}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle \label{2DpoL3} p_{3}(n) = -4(n+3)(\Ko+n+3)(\Kt+n+3)(\KM+n+3)(\KM+2n+2)(\KM+2n+3), \nonumber \end{align} \tag{ 9 }$$

which is negative for non-negative n (given that $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \Ko$ and $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kt}{\K_{2}} \Kt$ should actually read $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \vert \Ko\vert$ and $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kt}{\K_{2}} \vert \Kt\vert$, respectively). Hence, p₃ never vanishes, and the recurrence scheme never breaks down.

So, we have established that the generic third-order sixth-degree recurrence scheme in terms of the symbolic 2D spatial coordinates $\newcommand{\e}{{\rm e}} \ell_{1}$, $\newcommand{\e}{{\rm e}} \ell_{2}$ is computationally less complex than that of some of its lower order, higher degree counterparts, which is an example of trading order for degree [25, 29, 30]. Moreover, our generic 2D square lattice scheme is guaranteed to never break down. In [21], we have also shown that the recurrence may be evaluated in double precision arithmetic, provided that we augment the dynamic range on the fly. Further, we demonstrated that grid dispersion can be kept at bay for large grids by doubling the number of samples per wavelength for every factor of ten increase in grid distance.

Now, the question is whether we can generalise the recurrence scheme to three spatial dimensions. Below, we demonstrate that we can. However, we do not retain all of the nice properties germane to the 2D case. Nevertheless, every problem that we shall encounter can be overcome [31], and unexpected lattice artefacts can be explained [32] by drawing an analogy with singularity theory developed in the latter half of the twentieth century for continuous wave physics [33, 34].

Although Zeil only works for single sums, multivariate creative telescoping algorithms have also been developed [25, 35, 36]. For three spatial dimensions ($\nu=3$), it might seem prudent to first reduce the initial triple nested sums in (3) to the double nested sums in (4). For non-negative lattice points on the diagonal symmetry plane $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \newcommand{\Kt}{\K_{2}} \newcommand{\Kd}{\K_{3}} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK=(\Ko, \Kt=\Ko, \Kd)$, we have used the Maple package MultiZeilberger⁶ by Apagodu and Zeilberger [25] on the double nested sum in (4) (actually, to its factorial equivalent [21, equation 13]), which, after predicting a recurrence relation of order 9 produced a behemoth of a result. Fortunately, when we applied MultiZeilberger to the more symmetric triple nested sums in (3), we obtained a recurrence relation of the form (7) of order $\newcommand{\oL}{\Lambda} \oL=5$, and degree d  =  12, which amounts to a considerable reduction in complexity.

The pertaining polynomial coefficients $\newcommand{\ol}{u} p_{\ol}(n)$, with $\newcommand{\ol}{u} \ol=0, \ldots, 5$, may be factorised into linear polynomials in n, $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \Ko$, and $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kd}{\K_{3}} \Kd$, and six irreducible polynomials $\newcommand{\ol}{u} q_{\ol}(n)$ in n, $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \Ko$, and $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kd}{\K_{3}} \Kd$. They are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Kd}{\K_{3}} \newcommand{\Ko}{\K_{1}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{0}(n) \!=\!&\,\, {81\,q_{0}(n)} (2n+\KM+10)(2n+\KM+9)(2n+\KM+8)(2n+\KM+7)(n+1)\nonumber \\ & \times (n+\Ko+1)(n+\Kd+1)(n+2\Ko+1)(n+\Ko+\Kd+1)\nonumber \\ & \times (n+\KM+1) \nonumber \end{align} \tag{ 10a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{1}(n) \!=\! {9\,q_{1}(n)}(2n+\KM+10)(2n+\KM+9)(2n+\KM+8)(2n+\KM+3) \nonumber \end{align} \tag{ 10b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{2}(n) \!=\! {+2\,q_{2}(n)}(2n+\KM+10)(2n+\KM+9)(2n+\KM+5)(2n+\KM+2) \nonumber \end{align} \tag{ 10c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{3}(n) \!=\! {-2\,q_{3}(n)}(2n+\KM+10)(2n+\KM+7)(2n+\KM+3)(2n+\KM+2) \nonumber \end{align} \tag{ 10d }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{4}(n) \!=\! {-9\,q_{4}(n)}(2n+\KM+9)(2n+\KM+4)(2n+\KM+3)(2n+\KM+2) \nonumber \end{align} \tag{ 10e }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Kd}{\K_{3}} \newcommand{\Ko}{\K_{1}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{5}(n) \!=\! &\,\, {-81\,q_{5}(n)} (2n+\KM+5)(2n+\KM+4)(2n+\KM+3)(2n+\KM+2)(n+5)\nonumber \\ & {}\times (n+\Ko+5)(n+\Kd+5)(n+2\Ko+5)(n+\Ko+\Kd+5)\nonumber \\ & {}\times (n+\KM+5). \nonumber \end{align} \tag{ 10f }$$

Here, $\newcommand{\e}{{\rm e}} \ell=2\ell_{1}+\ell_{3}$, and the irreducible polynomials take the following form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ol}{u} \displaystyle \label{qolnSym} q_{\ol}(n) = \left\{\begin{array}{@{}ll@{}} \sum\limits_{i=0}^{2} q_{\ol,i}n^{i} &~{\rm for}~\ol\in\{0,5\}, \nonumber \\ \sum\limits_{i=0}^{8}q_{\ol,i}n^{i} & ~{\rm for}~ \ol\in\{1,\ldots,4\}, \end{array}\right. \nonumber \end{align} \tag{ 11 }$$

where the coefficients $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\ol}{u} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} q_{\ol, i}=q_{\ol, i}(\nvK)$ are polynomials of the form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ol}{u} \newcommand{\K}{\ell} \displaystyle \label{qoliSym} q_{\ol,i} = \sum_{j_{1},j_{3}} q_{\ol,i;j_{1},j_{3}}\K_{1}^{\,j_{1}}\K_{3}^{\,j_{3}}. \nonumber \end{align} \tag{ 12 }$$

In particular, $\newcommand{\ol}{u} q_{\ol, i}$ are polynomials in $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \Ko$ and $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kd}{\K_{3}} \Kd$ of degree 2  −  i, and 8  −  i for $\newcommand{\ol}{u} \ol\in\{0, 5\}$, and $\newcommand{\ol}{u} \ol\in\{1, \ldots\, 4\}$, respectively.

The construction of the second degree polynomials q₀(n) and q₅(n) is straightforward. For example, the polynomial coefficient q_0,0 in (11) can be assembled from (12) and the column below i  =  0 in the left-hand table in table 1, and reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\kd}{\Kd} \newcommand{\ko}{\Ko} \newcommand{\Kd}{\K_{3}} \newcommand{\Ko}{\K_{1}} \newcommand{\K}{\ell} \displaystyle q_{0,0}=2\kd^2+2\ko\kd-4\ko^{2}+7\kd+14\ko+20. \nonumber \end{align} \tag{ 13 }$$

Table 1. The polynomial coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \newcommand{\qCoefS}[1]{q_{#1,i;j_{1},j_{3}}} \qCoefS{0}$, and $\newcommand{\qCoef}[1]{q_{#1,i;j}} \newcommand{\qCoefS}[1]{q_{#1,i;j_{1},j_{3}}} \qCoefS{5}$ of the fifth-order recurrence scheme for lattice points on the symmetry planes in a simple cubic 3D lattice.

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \newcommand{\qCoefS}[1]{q_{#1,i;j_{1},j_{3}}} \qCoefS{0}$
$\newcommand{\jod}{j_{1},j_{3}} \jod$	i  =  0	i  =  1	i  =  2
0, 2	2
1, 1	2
2, 0	−4
0, 1	7	2
1, 0	14	4
0, 0	20	14	2

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \newcommand{\qCoefS}[1]{q_{#1,i;j_{1},j_{3}}} \qCoefS{5}$
$\newcommand{\jod}{j_{1},j_{3}} \jod$	i  =  0	i  =  1	i  =  2
0, 2	2
1, 1	2
2, 0	−4
0, 1	5	2
1, 0	10	4
0, 0	8	10	2

The polynomial coefficients $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\ol}{u} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} q_{\ol, i}(\nvK)$ for $\newcommand{\ol}{u} \ol\in\{1, \ldots, 4\}$ take up too much room if written out in full. Instead, the integer coefficients $\newcommand{\ol}{u} q_{\ol, i;j_{1}, j_{3}}$ are given in appendix A in tabular form.

For generic 3D lattice points, we again invoke the simple cubic lattice symmetry. In particular, we restrict (without loss of generality) our description to non-negative lattice points $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK$ that are ordered according to $\newcommand{\Gn}{\nu} \newcommand{\G}{G} \newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} 0\leqslant\Ko\leqslant\ldots\leqslant\K_{\Gn}$. Again, we tried to feed the double nested sums in (4) to MultiZeilberger. However, after predicting a recurrence relation of order 11, the package never produced a recurrence scheme. In an attempt to tackle the triple nested sums in (3) for generic symbolic lattice points, a test run of MultiZeilberger predicted a recurrence of order $\newcommand{\oL}{\Lambda} \oL=7$, but after consuming all of the available memory $32\, {\rm GB}$ within a couple of days, the process seemed to stall, and had to be terminated after a week.

Fortunately, it turned out that for an arbitrary selection of numeric off-diagonal-plane lattice points. MultiZeilberger did produce a recurrence relation of the form (7) of order $\newcommand{\oL}{\Lambda} \oL=7$, and fixed degree d  =  18.

Moreover, in contrast to the second-order, variable degree results returned by the single sum package Zeil in the 2D case, the numeric eighteenth-degree polynomial coefficients $\newcommand{\ol}{u} p_{\ol}(n)$ for randomly picked numeric lattice points as returned by MultiZeilberger all had the same structure, in that they presented as products of polynomials with integer coefficients that are linear in n and polynomials

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ol}{u} \displaystyle \label{qoLn} q_{\ol}(n) = \left\{\begin{array}{@{}ll@{}} \sum\limits_{i=0}^{4} q_{\ol,i}n^{i} & {\rm{for~}} \ol\in\{0,7\}, \nonumber \\ \sum\limits_{i=0}^{12}q_{\ol,i}n^{i} & {\rm{for~}} \ol\in\{1,\ldots,6\} \end{array}\right. \nonumber \end{align} \tag{ 14 }$$

that appear to be irreducible over $n\in\mathbb{Q}$, and have (numeric) integer coefficients $\newcommand{\ol}{u} q_{\ol, i}$.

Since the structure of the recurrences turned out to be the same for all the numeric lattice points we tested, we formulated an ansatz that the factorisation of the ${18}^{\mathrm th}$ degree polynomials for generic symbolic lattice points involves linear polynomials polynomials in n, $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \Ko$, $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kt}{\K_{2}} \Kt$, and $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kd}{\K_{3}} \Kd$, and polynomials $\newcommand{\ol}{u} q_{\ol}(n)$ with integer coefficients $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\ol}{u} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} q_{\ol, i}=q_{\ol, i}(\nvK)$ that are polynomials in $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \Ko$, $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kt}{\K_{2}} \Kt$ and $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kd}{\K_{3}} \Kd$ of degree $\newcommand{\oL}{\Lambda} 4-d_{\oL}$ with $\newcommand{\oL}{\Lambda} d_{\oL}=4$ for $\newcommand{\ol}{u} \ol\in\{0, 7\}$, and $\newcommand{\oL}{\Lambda} d_{\oL}=12$ for $\newcommand{\ol}{u} \ol\in\{1, \ldots\, 6\}$, respectively. Further, we invoked the invariance of the scalar LGF upon interchanging any of the spatial coordinates $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \Ko$, $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kt}{\K_{2}} \Kt$ and $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Kd}{\K_{3}} \Kd$ to restrict ourselves to a symmetric polynomial expansion

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ol}{u} \newcommand{\qCoef}[1]{q_{#1,i;j}} \displaystyle \label{qCoef} q_{\ol,i} = \sum_{j} \qCoef{\ol}r_{j}, \nonumber \end{align} \tag{ 15 }$$

in which j denotes a multi-index, i.e. $j=\{\,j_{1}, j_{2}, j_{3}\}$, where $0\leqslant j_{1}\leqslant j_{2}\leqslant j_{3}$ and $\newcommand{\oL}{\Lambda} j_{1}+j_{2}+j_{3}\leqslant d_{\oL}$, while the elementary symmetric polynomials are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\aKd}{\Kd} \newcommand{\aKt}{\Kt} \newcommand{\aKo}{\Ko} \newcommand{\Kd}{\K_{3}} \newcommand{\Kt}{\K_{2}} \newcommand{\Ko}{\K_{1}} \newcommand{\K}{\ell} \displaystyle r_{j} = \aKo^{\,j_{1}} \aKt^{\,j_{1}} \aKd^{\,j_{1}} \qquad \qquad \qquad \qquad \qquad {\rm{for~}} j_{1}=j_{2}=j_{3}, \nonumber \end{align} \tag{ 16a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\aKd}{\Kd} \newcommand{\aKt}{\Kt} \newcommand{\aKo}{\Ko} \newcommand{\Kd}{\K_{3}} \newcommand{\Kt}{\K_{2}} \newcommand{\Ko}{\K_{1}} \newcommand{\K}{\ell} \displaystyle r_{j} = \aKo^{\,j_{1}} \aKt^{\,j_{2}} \aKd^{\,j_{2}} +\aKo^{\,j_{2}} \aKt^{\,j_{1}} \aKd^{\,j_{2}} +\aKo^{\,j_{2}} \aKt^{\,j_{2}} \aKd^{\,j_{1}} \quad {\rm{for~}} j_{1}\neq j_{2}=j_{3}, \nonumber \end{align} \tag{ 16b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\aKd}{\Kd} \newcommand{\aKt}{\Kt} \newcommand{\aKo}{\Ko} \newcommand{\Kd}{\K_{3}} \newcommand{\Kt}{\K_{2}} \newcommand{\Ko}{\K_{1}} \newcommand{\K}{\ell} \displaystyle r_{j} = \aKo^{\,j_{1}} \aKt^{\,j_{1}} \aKd^{\,j_{3}} +\aKo^{\,j_{1}} \aKt^{\,j_{3}} \aKd^{\,j_{1}} +\aKo^{\,j_{3}} \aKt^{\,j_{1}} \aKd^{\,j_{1}} \qquad\quad\,\, \qquad \qquad \qquad \qquad {\rm{for~}} j_{1}=j_{2}\neq j_{3}, \nonumber \end{align} \tag{ 16c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\aKd}{\Kd} \newcommand{\aKt}{\Kt} \newcommand{\aKo}{\Ko} \newcommand{\Kd}{\K_{3}} \newcommand{\Kt}{\K_{2}} \newcommand{\Ko}{\K_{1}} \newcommand{\K}{\ell} \displaystyle r_{j} = & \aKo^{\,j_{1}} \aKt^{\,j_{2}} \aKd^{\,j_{3}} +\aKo^{\,j_{2}} \aKt^{\,j_{3}} \aKd^{\,j_{1}} +\aKo^{\,j_{3}} \aKt^{\,j_{1}} \aKd^{\,j_{2}} \quad & \nonumber \\ & \mbox{}+\aKo^{\,j_{2}} \aKt^{\,j_{1}} \aKd^{\,j_{3}} +\aKo^{\,j_{3}} \aKt^{\,j_{2}} \aKd^{\,j_{1}} +\aKo^{\,j_{1}} \aKt^{\,j_{3}} \aKd^{\,j_{2}} \quad & {\rm{for~}} j_{1}<j_{2}<j_{3}, \nonumber \end{align} \tag{ 16d }$$

where we would like to recall that the spatial coordinates $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \K_{i}$ are assumed non-negative.

The ansatz that the structure of the recurrence operator remains the same regardless of the specific (off-diagonal-plane) lattice point prompted us to undertake a systematic search for the generic form of the coefficients. There are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sum\limits_{{0\leqslant j_{1}\leqslant j_{2}\leqslant j_{3}}\atop{j_{1}+j_{2}+j_{3}\leqslant12}} 1 = 102 \nonumber \end{align} \tag{ 17 }$$

elementary symmetric polynomials. In order to extract the symbolic polynomial coefficients, we generated the recurrence relations pertaining to an off-diagonal cube of 125 numeric lattice points, $(s_{1}, 1600+s_{2}, 9900+s_{3})$, $s_{1}, s_{2}, s_{3}=0, \ldots, 4$. The resulting overdetermined linear system of equations was indeed consistent and resulted in the unique solution

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Kd}{\K_{3}} \newcommand{\Kt}{\K_{2}} \newcommand{\Ko}{\K_{1}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{0}(n) \!=\! &\,\,{729\,q_{0}(n)}(2n+\KM+14)(2n+\KM+13)(2n+\KM+12)(2n+\KM+11)\nonumber \\ & {}\times (2n+\KM+10)(2n+\KM+9) (n+1)\nonumber \\ & {}\times (n+\Ko+1)(n+\Kt+1)(n+\Kd+1)\nonumber \\ & {}\times (n+\Kt+\Kd+1)(n+\Kd+\Ko+1)(n+\Ko+\Kt+1)\nonumber \\ & {}\times (n+\KM+1), \nonumber \end{align} \tag{ 18a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{1}(n) \!=\!&\,\, {81\,q_{1}(n)}(2n+\KM+14)(2n+\KM+13)(2n+\KM+12)(2n+\KM+11)\nonumber \\ & {}\times (2n+\KM+10) (2n+\KM+3), \nonumber \end{align} \tag{ 18b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{2}(n) \!=\!&\,\, {9\,q_{2}(n)}(2n+\KM+14)(2n+\KM+13)(2n+\KM+12)(2n+\KM+11)\nonumber \\ & {}\times (2n+\KM+5) (2n+\KM+2), \nonumber \end{align} \tag{ 18c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{3}(n) \!=\!&\,\, {\, q_{3}(n)}(2n+\KM+14)(2n+\KM+13)(2n+\KM+12) \nonumber \\ & {}\times (2n+\KM+7) (2n+\KM+3)(2n+\KM+2), \nonumber \end{align} \tag{ 18d }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{4}(n) \!=\!&\,\, {- q_{4}(n)}(2n+\KM+14)(2n+\KM+13) (2n+\KM+9)\nonumber \\ & {}\times (2n+\KM+4)(2n+\KM+3)(2n+\KM+2), \nonumber \end{align} \tag{ 18e }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{5}(n) \!=\!&\,\, {-9\,q_{5}(n)}(2n+\KM+14) (2n+\KM+11) (2n+\KM+5)\nonumber \\ & {}\times (2n+\KM+4)(2n+\KM+3)(2n+\KM+2), \nonumber \end{align} \tag{ 18f }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{6}(n) \!=\!&\,\, {-81\,q_{6}(n)}(2n+\KM+13) (2n+\KM+6)(2n+\KM+5)\nonumber \\ & {}\times (2n+\KM+4)(2n+\KM+3)(2n+\KM+2), \nonumber \end{align} \tag{ 18g }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Kd}{\K_{3}} \newcommand{\Kt}{\K_{2}} \newcommand{\Ko}{\K_{1}} \newcommand{\KM}{\K} \newcommand{\K}{\ell} \displaystyle p_{7}(n) \!=\!&\,\, {-729\,q_{7}(n)}\, (2n+\KM+7)(2n+\KM+6)(2n+\KM+5)\nonumber \\ & {}\times (2n+\KM+4)(2n+\KM+3)(2n+\KM+2)(n+7)\nonumber \\ & {}\times (n+\Ko+7)(n+\Kt+7)(n+\Kd+7)\nonumber \\ & {}\times (n+\Kt+\Kd+7)(n+\Kd+\Ko+7)(n+\Ko+\Kt+7)\nonumber \\ & {}\times (n+\KM+7), \nonumber \end{align} \tag{ 18h }$$

in which $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\KM}{\K} \newcommand{\Ko}{\K_{1}} \newcommand{\Kt}{\K_{2}} \newcommand{\Kd}{\K_{3}} \KM=\Ko+\Kt+\Kd$. For $\newcommand{\ol}{u} \ol\in\{1, \ldots, 6\}$, the polynomial coefficients $\newcommand{\ol}{u} q_{\ol}(n)$ defined through (14)–(16) take up way too much room if written out in full. Instead, the integer coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \newcommand{\ol}{u} \qCoef{\ol}$ of the elementary symmetric polynomials r_j are given in appendix B in tabular form. For $\newcommand{\ol}{u} \ol\in\{0, 7\}$, the polynomial coefficients do not take up quite as much room. For example, $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} q_{0, 1}=q_{0, 1}(\nvK)$ can be constructed from (16) and the column below i  =  1 in the left-hand table in table 2, and reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\kd}{\Kd} \newcommand{\kt}{\Kt} \newcommand{\ko}{\Ko} \newcommand{\Kd}{\K_{3}} \newcommand{\Kt}{\K_{2}} \newcommand{\Ko}{\K_{1}} \newcommand{\K}{\ell} \displaystyle q_{0,1}\!=\!-2(\ko^3+\kt^3+\kd^3)+6\kt\kd\ko+54(\kt\kd+\ko\kd+\ko\kt) +244(\ko+\kt+\kd)+738. \nonumber \end{align} \tag{ 19 }$$

Table 2. The polynomial coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{0}$, and $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{7}$ of the seventh order recurrence scheme for generic points on a simple cubic 3D lattice. Here, j is the multi-index of r_j in (16).

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{0}$
j	i  =  0	i  =  1	i  =  2	i  =  3	i  =  4
0, 0, 4	−4
0, 1, 3	−1
0, 2, 2	6
1, 1, 2
0, 0, 3	−9	−2
0, 1, 2
1, 1, 1	27	6
0, 0, 2	19
0, 1, 1	122	54	6
0, 0, 1	369	244	54	4
0, 0, 0	756	738	244	36	2

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{7}$
j	i  =  0	i  =  1	i  =  2	i  =  3	i  =  4
0, 0, 4	−4
0, 1, 3	−1
0, 2, 2	6
1, 1, 2
0, 0, 3	−7	−2
0, 1, 2
1, 1, 1	21	6
0, 0, 2	19
0, 1, 1	74	42	6
0, 0, 1	175	148	42	4
0, 0, 0	228	350	148	28	2

Below, we comment on the LGF recurrence scheme through an elucidating example.

Let us illustrate our LGF recurrence scheme through a (not quite) numerical example. Not quite numerical in that we have used integer arithmetic in C++ using the GNU Multiple Precision Arithmetic Library GMP [37] to generate an LGF sequence of rational numbers that are exact within the context of the simple cubic lattice FDTD model with unit Courant number ($\newcommand{\Ga}{\alpha} \newcommand{\G}{G} \Ga=1$). Traditionally, the absolute discrete time index in FDTD is denoted as n. Following this convention, we choose to relabel the sequence index in (7) as m, and rewrite $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} G(m)=G_{\nvK}^{n}$ with n  =  n_O  +  2m, and n_O defined in (2) is the travel time in the Manhattan metric associated with the first exponentially small wave motion.

The various features of the lattice Green's function are more easily distinguished at relatively large lattice distances than at shorter ones. In our example, we have chosen the lattice point $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK=(11\, 267, 40\, 622, 599\, 89)$, and processed the recurrence scheme to generate the LGF sequence for three million time steps. In figure 1, we have depicted the most interesting part of the LGF sequence with three segments magnified by a factor 100. The LGF sequence starts at $n=n_{O}=111\, 879$, with $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} G^{n}_{\nvK}\approx1.129\times10^{-8042}$, followed by rapid growth until the 'arrival' of the dominant wave motion at about $\newcommand{\Gn}{\nu} \newcommand{\G}{G} \newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} n=n_{S}=\sqrt{\Gn}\vert \nvK\vert =\sqrt{3(l_{1}^{2}+l_{2}^{2}+l_{3}^{2})}\approx126\, 993$. One might expect everything to fizzle out after that, but in addition to a very long tail, a distinct anomalous event seems to occur at about $n=n_{A}\approx179\, 284$. Finally, we have zoomed in on a segment about $n=n_{B}=184\, 972$, which hardly seems of any interest, but still merits discussion in the context of recurrence breakdown. Below, we briefly comment on our observations.

Zoom In Zoom Out Reset image size

• (i)  
  If we had wanted to compute exactly the same LGF sequence for three million time steps using the standard (albeit with integer arithmetic) FDTD scheme given by (1) on a simple cubic grid with imperfect boundary conditions, then we would have had to use a computational domain with more than $28\times10^{18}$ grid points, so as to prevent reflections at the boundaries from contaminating the exact results. Of course, FDTD computations on such huge grids are (currently) inconceivable, and irrelevant if sufficiently accurate absorbing boundary conditions are available, or better, if one applies discrete Green's function diakoptics on disjoint computational domains proposed in [14, 17, 18, 20, 38].Through a mixed multi-precision and floating point implementation Stefanski [19] optimised the triple combinatorial nested sums for 3D electromagnetic LGFs derived by Jeng [39], so as to demonstrate the feasibility of discrete Green's function diakoptics for moderately sized domains and short sequences. The computational overhead we envisaged in following such a brute-force approach for computing LGF sequences was what motivated us to embark on a slow route towards our current considerably faster alternative [17, 18] for the scalar LGF, and for its 3D electromagnetic counterpart [40].
• (ii)  
  For lattice points near to the origin, where standard FDTD is feasible, our recurrence scheme and FDTD yield identical results in integer arithmetic as long as reflections from the boundaries have not returned to the lattice point yet. Here, we should point out that exact integer arithmetic computations are not essential. However, integer computations are useful for generating absolute benchmarks for the much faster finite precision arithmetic computations to be discussed in [31].
• (iii)  
  Figure 1 clearly shows that the discrete counterpart of the continuous space-time 3D Green's function for the scalar wave equation $\newcommand{\Gd}{\delta} \newcommand{\G}{G} \Gd(t-R/c)/(4\pi R)$ is not a localised distribution. A meaningful comparison and a quantitative assessment of the numerical grid dispersion at large distances, similar to that performed in [21] for 2D LGFs can be made in a weak sense through convolution with sufficiently smooth pulses.
• (iv)  
  In contrast to the 2D case, a 3D LGF sequence appears to have an additional feature. We are alluding to the apparent distinct anomalous event at about $n=n_{A}\approx179\, 284$. We have argued in [41, 42] that this is a high-frequency feature of 3D FDTD, rather than some artefact induced by the recurrence scheme, and that it can be fully suppressed provided the Green's function is convolved with a sufficiently smooth pulse. In [32], we shall take inspiration from singularity theory for continuous-wave physics [33, 34], and the discrete Gabor transform [43] to explore asymptotically, and explain comprehensively, the various features of the lattice Green's function for large lattice distances. In particular, we shall identify the late-time apparently anomalous event as a high-frequency interplay of several resonances.
• (v)  
  The chosen lattice point $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK=\nvK_{B}$ is special in the sense that the generic seventh-order recurrence breaks down at $m=m_{B}=36\, 540$, because $q_{7}(m_{B})=0$, and hence $p_{7}(m_{B})=0$. At m  =  m_B the recurrence relation is said to be singular [44–46]. The sequence index m  =  m_B corresponds to the time index $n=n_{B}=184\, 972$. In [31], we shall show that the chosen singular lattice point $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK_{B}$ is atypical in that it does not belong to a larger class of singular lattice points for which $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\Ko}{\K_{1}} \newcommand{\Kt}{\K_{2}} \newcommand{\Kd}{\K_{3}} \Kd=\Ko+\Kt$.We have circumvented the recurrence breakdown by running our no-neighbours recurrence scheme for the six neighbours of $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK_{B}$ until n  =  n_B, and by using the FDTD scheme (1) to determine the LGF at $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK=\nvK_{B}$. Beyond n  =  n_B, we proceeded using the recurrence scheme for $\newcommand{\e}{{\rm e}} \newcommand{\K}{\ell} \newcommand{\nvk}{\nv{\K}} \newcommand{\nvK}{\nvk} \newcommand{\nv}[1]{\vec{#1}} \nvK_{B}$. The FDTD scheme is a nonsingular second-order recurrence scheme in its own right, albeit a non-local one (involving neighbouring lattice points).As can be seen from the zoomed-in segment of the LGF sequence about n  =  n_B, the recurrence breakdown seems to have no impact on the behaviour of the LGF. This observation is supported by the fact that for lattice points on the symmetry planes, we may choose to use either the generic seventh-order scheme, or the fifth-order one for symmetry planes, which can be shown to never break down simultaneously [31]. This implies that at the singular lattice points for the fifth-order symmetry-plane scheme, the generic seventh-order recurrence operator is a nonsingular left multiple of the singular fifth-order scheme, and hence that these integer singularities must be apparent, or removable [47]. Our circumvention of recurrence breakdown through the FDTD scheme suggests that all integer singularities of the generic seventh-order recurrence scheme are removable, provided that the singular lattice points are isolated.The use of the FDTD scheme up to n  =  n_B to circumvent recurrence breakdown does add to the computational expense. However, in [31], we shall show how sparsely the singular lattice points for the generic recurrence scheme are distributed in $\mathbb{Z}^{3}$, and we shall provide evidence for what we call the strong nearest-neighbour conjecture, which, if proven, would imply that the singular lattice points are always isolated and hence removable. Further, in [31] we shall discuss the details concerning the implementation of the recurrence scheme, e.g. the required initial conditions, and the number of digits needed in fixed-precision computations for retaining double precision accuracy.

In this paper we have discussed no-neighbours recurrence schemes for on-the-fly generation of discrete space-time lattice Green's function (LGF) sequences for 3D scalar wavefields on a simple cubic lattice at the maximum time step allowed for stability ($\newcommand{\Ga}{\alpha} \newcommand{\G}{G} \Ga=1$). Thus, applications such as discrete Green's function diakoptics for distant disjoint computational domains have become feasible.

In addition to opening up a field of applications, the LGF sequences are rich in features, in particular a late-time high-frequency effect that will be fully explained in [32]. At certain lattice points the recurrence scheme breaks down at some discrete instant in time n  =  n_B. However, we have described how to circumvent this through application of the FDTD stencil to the neighbours of the singular lattice point up to the instant n_B. In [31], we shall discuss recurrence breakdown and finite-precision implementations aspects in detail. Finally, in [40], we shall elaborate on the evaluation of vectorial lattice Green's function sequences for electromagnetic wavefields.

Other extensions and applications, especially, the consideration of smaller time steps ($\newcommand{\Ga}{\alpha} \newcommand{\G}{G} \Ga<1$), and a comprehensive implementation of discrete Green's function diakoptics remain to be investigated.

The research presented above has been made possible through the generous support of Maplesoft. In particular, we would like to thank Gosé Fischer of CANdiensten for making a 30-day evaluation of Maple 18 available in July 2014. The research of the second author (JMA) was partly supported by Engineering and Physical Sciences Research Council Science and Innovation Research Grant: EP/D501288/1 'Electronics Design Centre for Heterogeneous Systems'.

In tables A1–A4, we have tabulated the coefficients for the construction, through (11) and (12), of the eighth degree polynomials $\newcommand{\ol}{u} q_{\ol}(n)$ with $\newcommand{\ol}{u} \ol\in\{1, \ldots, 4\}$. See the supplementary information (stacks.iop.org/JPhysA/51/085201/mmedia) for further details.

Table A1. The polynomial coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \newcommand{\qCoefS}[1]{q_{#1,i;j_{1},j_{3}}} \qCoefS{1}$ of the fifth-order recurrence scheme for lattice points on the symmetry planes in a simple cubic 3D lattice.

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \newcommand{\qCoefS}[1]{q_{#1,i;j_{1},j_{3}}} \qCoefS{1}$
$\newcommand{\jod}{j_{1},j_{3}} \jod$	i  =  0	i  =  1	i  =  2	i  =  3	i  =  4	i  =  5	i  =  6	i  =  7	i  =  8
0, 8
1, 7	−8
2, 6	−32
3, 5
4, 4	112
5, 3	40
6, 2	−144
7, 1	−32
8, 0	64
0, 7	−12	−8
1, 6	−76	−32
2, 5	−32	44
3, 4	580	384
4, 3	1020	492
5, 2	−184	−224
6, 1	−1072	−592
7, 0	−224	−64
0, 6	−74	−76	−16
1, 5	32	108	66
2, 4	2392	2726	768
3, 3	6802	7038	1742
4, 2	5020	4440	888
5, 1	−1984	−2696	−912
6, 0	−2144	−2144	−592
0, 5	14	32	54	22
1, 4	3649	5898	3219	576
2, 3	18 303	28 418	14 459	2354
3, 2	29 962	43 828	20 748	3128
4, 1	14 888	18 672	7200	792
5, 0	−1584	−3968	−2696	−608
0, 4	1735	3649	2949	1073	144
1, 3	19 936	41 800	32 424	10 940	1330
2, 2	60 862	124 742	93 270	29 936	3442
3, 1	64 710	127 336	90 672	27 420	2950
4, 0	16 980	29 776	18 672	4800	396
0, 3	7509	19 936	20 900	10 808	2735	266
1, 2	51 767	137 460	142 170	71 356	17 265	1596
2, 1	99 007	259 886	264 408	129 704	30 445	2718
3, 0	51 018	129 420	127 336	60 448	13 710	1180
0, 2	15 827	51 767	68 730	47 390	17 839	3453	266
1, 1	67 158	220 494	292 764	200 740	74 610	14 154	1064
2, 0	61 028	198 014	259 886	176 272	64 852	12 178	906
0, 1	17 029	67 158	110 247	97 588	50 185	14 922	2359	152
1, 0	34 058	134 316	220 494	195 176	100 370	29 844	4718	304
0, 0	7340	34 058	67 158	73 498	48 794	20 074	4974	674	38

Table A2. The polynomial coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \newcommand{\qCoefS}[1]{q_{#1,i;j_{1},j_{3}}} \qCoefS{2}$ of the fifth-order recurrence scheme for lattice points on the symmetry planes in a simple cubic 3D lattice.

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \newcommand{\qCoefS}[1]{q_{#1,i;j_{1},j_{3}}} \qCoefS{2}$
$\newcommand{\jod}{j_{1},j_{3}} \jod$	i  =  0	i  =  1	i  =  2	i  =  3	i  =  4	i  =  5	i  =  6	i  =  7	i  =  8
0, 8	−32
1, 7	−116
2, 6	16
3, 5	352
4, 4	88
5, 3	−284
6, 2	−104
7, 1	−80
8, 0	160
0, 7	−370	−116
1, 6	−1322	−400
2, 5	402	228
3, 4	5606	2008
4, 3	3712	1424
5, 2	−4668	−1416
6, 1	−4048	−1568
7, 0	688	−160
0, 6	−1813	−1322	−200
1, 5	−2751	−1728	−186
2, 4	16 677	12 792	2364
3, 3	41 369	30 730	5530
4, 2	12 462	11 052	2436
5, 1	$-27\, 324$	$-17\, 184$	−2544
6, 0	$-12\, 376$	−8096	−1568
0, 5	−1924	−2751	−864	−62
1, 4	27 968	29 076	10 779	1360
2, 3	143 771	158 913	57 585	6790
3, 2	176 281	200 922	74 880	9080
4, 1	3496	15 492	10 884	2024
5, 0	$-57\, 020$	$-54\, 648$	$-17\, 184$	−1696
0, 4	21 170	27 968	14 538	3593	340
1, 3	219 673	322 834	176 316	42 220	3710
2, 2	552 144	831 297	459 945	110 460	9690
3, 1	406 594	623 770	354 552	88 300	8090
4, 0	−2644	6992	15 492	7256	1012
0, 3	117 752	219 673	161 417	58 772	10 555	742
1, 2	680 934	1310 742	982 524	359 484	64 125	4452
2, 1	1099 518	2149 386	1631 265	602 064	108 105	7554
3, 0	417 316	813 188	623 770	236 368	44 150	3236
0, 2	285 343	680 934	655 371	327 508	89 871	12 825	742
1, 1	1097 902	2672 580	2614 188	1319 380	362 910	51 618	2968
2, 0	905 392	2199 036	2149 386	1087 510	301 032	43 242	2518
0, 1	371 082	1097 902	1336 290	871 396	329 845	72 582	8603	424
1, 0	742 164	2195 804	2672 580	1742 792	659 690	145 164	17 206	848
0, 0	209 520	742 164	1097 902	890 860	435 698	131 938	24 194	2458	106

Table A3. The polynomial coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \newcommand{\qCoefS}[1]{q_{#1,i;j_{1},j_{3}}} \qCoefS{3}$ of the fifth-order recurrence scheme for lattice points on the symmetry planes in a simple cubic 3D lattice.

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \newcommand{\qCoefS}[1]{q_{#1,i;j_{1},j_{3}}} \qCoefS{3}$
$\newcommand{\jod}{j_{1},j_{3}} \jod$	i  =  0	i  =  1	i  =  2	i  =  3	i  =  4	i  =  5	i  =  6	i  =  7	i  =  8
0, 8	−32
1, 7	−116
2, 6	16
3, 5	352
4, 4	88
5, 3	−284
6, 2	−104
7, 1	−80
8, 0	160
0, 7	−326	−116
1, 6	−1078	−400
2, 5	966	228
3, 4	6442	2008
4, 3	4832	1424
5, 2	−3828	−1416
6, 1	−5360	−1568
7, 0	−1648	−160
0, 6	−1081	−1078	−200
1, 5	921	−504	−186
2, 4	25 029	15 576	2364
3, 3	56 069	35 630	5530
4, 2	33 846	18 180	2436
5, 1	$-15\, 804$	$-13\, 344$	−2544
6, 0	$-20\, 248$	$-10\, 720$	−1568
0, 5	3130	921	−252	−62
1, 4	52 204	46 608	13 701	1360
2, 3	203 287	201 213	64 635	6790
3, 2	294 851	283 002	88 560	9080
4, 1	134 816	103 476	25 548	2024
5, 0	$-18\, 580$	$-31\, 608$	$-13\, 344$	−1696
0, 4	41 282	52 204	23 304	4567	340
1, 3	318 685	438 638	217 716	46 820	3710
2, 2	821 262	1130 523	564 705	122 100	9690
3, 1	839 686	1084 214	512 592	105 860	8090
4, 0	257 372	269 632	103 476	17 032	1012
0, 3	174 538	318 685	219 319	72 572	11 705	742
1, 2	973 950	1789 686	1253 508	423 204	69 435	4452
2, 1	1752 906	3144 318	2171 847	726 984	118 515	7554
3, 0	1006 316	1679 372	1084 214	341 728	52 930	3236
0, 2	415 735	973 950	894 843	417 836	105 801	13 887	742
1, 1	1611 910	3749 244	3456 948	1629 740	417 090	55 230	2968
2, 0	1554 400	3505 812	3144 318	1447 898	363 492	47 406	2518
0, 1	561 234	1611 910	1874 622	1152 316	407 435	83 418	9205	424
1, 0	1122 468	3223 820	3749 244	2304 632	814 870	166 836	18 410	848
0, 0	321 840	1122 468	1611 910	1249 748	576 158	162 974	27 806	2630	106

Table A4. The polynomial coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \newcommand{\qCoefS}[1]{q_{#1,i;j_{1},j_{3}}} \qCoefS{4}$ of the fifth-order recurrence scheme for lattice points on the symmetry planes in a simple cubic 3D lattice.

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \newcommand{\qCoefS}[1]{q_{#1,i;j_{1},j_{3}}} \qCoefS{4}$
$\newcommand{\jod}{j_{1},j_{3}} \jod$	i  =  0	i  =  1	i  =  2	i  =  3	i  =  4	i  =  5	i  =  6	i  =  7	i  =  8
0, 8
1, 7	−8
2, 6	−32
3, 5
4, 4	112
5, 3	40
6, 2	−144
7, 1	−32
8, 0	64
0, 7	−36	−8
1, 6	−116	−32
2, 5	296	44
3, 4	1724	384
4, 3	1932	492
5, 2	−1160	−224
6, 1	−2480	−592
7, 0	−160	−64
0, 6	−194	−116	−16
1, 5	1760	684	66
2, 4	13 684	6490	768
3, 3	27 286	13 866	1742
4, 2	10 348	6216	888
5, 1	$-18\, 640$	−8248	−912
6, 0	$-10\, 592$	−4960	−592
0, 5	2986	1760	342	22
1, 4	40 271	29 478	7149	576
2, 3	140 145	109 142	27 913	2354
3, 2	161 726	132 676	35 556	3128
4, 1	9016	17 808	7056	792
5, 0	$-56\, 496$	$-37\, 280$	−8248	−608
0, 4	40 861	40 271	14 739	2383	144
1, 3	297 040	314 888	122 784	20 980	1330
2, 2	664 786	735 298	297 894	52 672	3442
3, 1	465 366	548 168	234 312	43 380	2950
4, 0	$-13\, 068$	18 032	17 808	4704	396
0, 3	218 091	297 040	157 444	40 928	5245	266
1, 2	1102 825	1562 988	860 358	231 556	30 615	1596
2, 1	1636 133	2403 182	1365 024	377 504	51 095	2718
3, 0	605 694	930 732	548 168	156 208	21 690	1180
0, 2	622 577	1102 825	781 494	286 786	57 889	6123	266
1, 1	2198 898	3999 858	2906 724	1090 940	224 550	24 150	1064
2, 0	1776 488	3272 266	2403 182	910 016	188 752	20 438	906
0, 1	978 515	2198 898	1999 929	968 908	272 735	44 910	4025	152
1, 0	1957 030	4397 796	3999 858	1937 816	545 470	89 820	8050	304
0, 0	702 200	1957 030	2198 898	1333 286	484 454	109 094	14 970	1150	38

In tables B1–B6, we have tabulated the coefficients for the construction, through (14)–(16), of the twelfth degree polynomials $\newcommand{\ol}{u} q_{\ol}(n)$ with $\newcommand{\ol}{u} \ol\in\{1, \ldots, 6\}$. See the supplementary information for further details.

Table B1. The polynomial coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{1}$ of the seventh-order recurrence scheme for generic points on a simple cubic 3D lattice. Here, j is the multi-index of r_j in (16).

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{1}$
j	i  =  0	i  =  1	i  =  2	i  =  3	i  =  4	i  =  5
0, 2, 10	4
1, 1, 10	8
0, 3, 9	17
1, 2, 9	19
0, 4, 8	22
1, 3, 8	4
2, 2, 8	−36
0, 5, 7	−1
1, 4, 7	−6
2, 3, 7	−141
0, 6, 6	−20
1, 5, 6	23
2, 4, 6	−36
3, 3, 6	−158
2, 5, 5	108
3, 4, 5	106
4, 4, 4	228
0, 1, 10	20	16
0, 2, 9	58	38
1, 1, 9	20	12
0, 3, 8	49	8
1, 2, 8	−333	−204
0, 4, 7	33	−12
1, 3, 7	−711	−403
2, 2, 7	−1872	−1038
0, 5, 6	80	46
1, 4, 6	−27	16
2, 3, 6	−1797	−918
1, 5, 5	702	438
2, 4, 5	1116	666
3, 3, 5	700	356
3, 4, 4	3474	1726
0, 0, 10	24	40	16
0, 1, 9	−1	40	12
0, 2, 8	−681	−666	−204
1, 1, 8	−2262	−2292	−636
0, 3, 7	−1492	−1422	−403
1, 2, 7	−8833	−8857	−2268
0, 4, 6	−455	−54	16
1, 3, 6	−8429	−7760	−1806
2, 2, 6	$-18\, 116$	$-18\, 044$	−4380
0, 5, 5	762	1404	438
1, 4, 5	3053	4601	1362
2, 3, 5	106	−275	6
2, 4, 4	20 140	19 156	4644
3, 3, 4	32 916	28 522	6488
0, 0, 9	−66	−2	40	8
0, 1, 8	−3062	−4524	−2292	−424
0, 2, 7	$-11\, 815$	$-17\, 666$	−8857	−1512
1, 1, 7	$-28\, 797$	$-44\, 046$	$-22\, 305$	−3730
0, 3, 6	$-12\, 613$	$-16\, 858$	−7760	−1204
1, 2, 6	$-56\, 229$	$-86\, 505$	$-43\, 350$	−6924
0, 4, 5	396	6106	4601	908
1, 3, 5	−4799	−9150	−5025	−700
2, 2, 5	$-14\, 874$	$-40\, 998$	$-26\, 310$	−4476
1, 4, 4	51 390	75 666	36 624	5836
2, 3, 4	159 965	192 307	77 387	10 812
3, 3, 3	306 087	370 992	149 805	20 790
0, 0, 8	−3084	−6124	−4524	−1528	−212
0, 1, 7	$-27\, 738$	$-57\, 594$	$-44\, 046$	$-14\, 870$	−1865
0, 2, 6	$-52\, 916$	$-112\, 458$	$-86\, 505$	$-28\, 900$	−3462
1, 1, 6	$-117\, 043$	$-261\, 696$	$-208\, 941$	$-71\, 180$	−8580
0, 3, 5	−1373	−9598	−9150	−3350	−350
1, 2, 5	7556	$-90\, 485$	$-123\, 630$	$-52\, 070$	−6720
0, 4, 4	53 778	102 780	75 666	24 416	2918
1, 3, 4	428 438	660 039	379 425	97 730	9730
2, 2, 4	766 349	1130 068	603 357	140 308	12 546
2, 3, 3	1418 750	2228 484	1287 402	328 000	31 710
0, 0, 7	$-20\, 767$	$-55\, 476$	$-57\, 594$	$-29\, 364$	−7435	−746
0, 1, 6	$-77\, 545$	$-234\, 086$	$-261\, 696$	$-139\, 294$	$-35\, 590$	−3432
0, 2, 5	62 523	15 112	$-90\, 485$	$-82\, 420$	$-26\, 035$	−2688
1, 1, 5	170 955	53 574	$-242\, 661$	$-228\, 960$	$-73\, 080$	−7644
0, 3, 4	441 436	856 876	660 039	252 950	48 865	3892
1, 2, 4	1817 208	3404 982	2430 414	822 100	133 075	8568
1, 3, 3	3148 927	6264 086	4848 930	1832 820	342 195	25 578
2, 2, 3	5531 535	11 067 562	8578 159	3228 658	595 690	43 680
0, 0, 6	$-37\, 279$	$-155\, 090$	$-234\, 086$	$-174\, 464$	$-69\, 647$	$-14\, 236$
0, 1, 5	232 450	341 910	53 574	$-161\, 774$	$-114\, 480$	$-29\, 232$
0, 2, 4	1548 076	3634 416	3404 982	1620 276	411 050	53 230
1, 1, 4	3522 117	8254 632	7644 318	3540 750	853 725	101 052
0, 3, 3	2563 498	6297 854	6264 086	3232 620	916 410	136 878
1, 2, 3	9948 301	24 710 184	24 690 186	12 699 350	3555 585	518 028
2, 2, 2	17 064 909	42 936 588	43 498 584	22 714 476	6464 199	956 112
0, 0, 5	208 763	464 900	341 910	35 716	$-80\, 887$	$-45\, 792$
0, 1, 4	2461 558	7044 234	8254 632	5096 212	1770 375	341 490
0, 2, 3	6636 482	19 896 602	24 710 184	16 460 124	6349 675	1422 234
1, 1, 3	14 699 257	44 308 818	55 236 990	36 852 200	14 200 095	3167 010
1, 2, 2	24 771 021	75 617 568	95 620 554	64 862 044	25 493 555	5819 844
0, 0, 4	1445 256	4923 116	7044 234	5503 088	2548 106	708 150
0, 1, 3	8250 003	29 398 514	44 308 818	36 824 660	18 426 100	5680 038
0, 2, 2	13 744 818	49 542 042	75 617 568	63 747 036	32 431 022	10 197 422
1, 1, 2	29 987 586	108 807 234	167 162 898	141 820 740	72 602 880	22 972 656
0, 0, 3	3985 461	16 500 006	29 398 514	29 539 212	18 412 330	7370 440
0, 1, 2	14 241 024	59 975 172	108 807 234	111 441 932	70 910 370	29 041 152
1, 1, 1	30 766 689	130 425 498	238 226 160	245 708 160	157 494 120	65 012 136
0, 0, 2	5907 591	28 482 048	59 975 172	72 538 156	55 720 966	28 364 148
0, 1, 1	12 682 314	61 533 378	130 425 498	158 817 440	122 854 080	62 997 648
0, 0, 1	4619 313	25 364 628	61 533 378	86 950 332	79 408 720	49 141 632
0, 0, 0	1500 660	9238 626	25 364 628	41 022 252	43 475 166	31 763 488

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{1}$
j	i  =  6	i  =  7	i  =  8	i  =  9	i  =  10	i  =  11	i  =  12
0, 0, 6	−1144
0, 1, 5	−2548
0, 2, 4	2856
1, 1, 4	4676
0, 3, 3	8526
1, 2, 3	30 996
2, 2, 2	57 792
0, 0, 5	−9744	−728
0, 1, 4	33 684	1336
0, 2, 3	172 676	8856
1, 1, 3	381 150	19 260
1, 2, 2	717 836	37 056
0, 0, 4	113 830	9624	334
0, 1, 3	1055 670	108 900	4815
0, 2, 2	1939 948	205 096	9264
1, 1, 2	4397 610	467 280	21 150
0, 0, 3	1893 346	301 620	27 225	1070
0, 1, 2	7657 552	1256 460	116 820	4700
1, 1, 1	17 292 618	2864 520	268 785	10 890
0, 0, 2	9680 384	2187 872	314 115	25 960	940
0, 1, 1	21 670 712	4940 748	716 130	59 730	2178
0, 0, 1	20 999 216	6191 632	1235 187	159 140	11 946	396
0, 0, 0	16 380 544	5999 776	1547 908	274 486	31 828	2172	66

Table B2. The polynomial coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{2}$ of the seventh-order recurrence scheme for generic points on a simple cubic 3D lattice. Here, j is the multi-index of r_j in (16).

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{2}$
j	i  =  0	i  =  1	i  =  2	i  =  3	i  =  4	i  =  5
0, 1, 11	128
0, 2, 10	364
1, 1, 10	472
0, 3, 9	67
1, 2, 9	297
0, 4, 8	−478
1, 3, 8	−692
2, 2, 8	−1452
0, 5, 7	109
1, 4, 7	−498
2, 3, 7	−2167
0, 6, 6	836
1, 5, 6	1205
2, 4, 6	52
3, 3, 6	−1402
2, 5, 5	1668
3, 4, 5	1182
4, 4, 4	2220
0, 0,11	448	256
0, 1, 10	2436	944
0, 2, 9	2866	594
1, 1, 9	4356	1380
0, 3, 8	−2187	−1384
1, 2, 8	$-10\, 761$	−4460
0, 4, 7	−2619	−996
1, 3, 7	$-25\, 845$	−9419
2, 2, 7	$-51\, 252$	$-17\, 102$
0, 5, 6	6144	2410
1, 4, 6	2265	944
2, 3, 6	$-41\, 001$	$-12\, 738$
1, 5, 5	30 858	10 934
2, 4, 5	32 136	10 730
3, 3, 5	18 012	7292
3, 4, 4	70 710	23 990
0, 0, 10	5256	4872	944
0, 1, 9	8845	8712	1380
0, 2, 8	$-31\, 043$	$-21\, 522$	−4460
1, 1, 8	$-79\, 578$	$-51\, 444$	−9228
0, 3, 7	$-71\, 354$	$-51\, 690$	−9419
1, 2, 7	$-326\, 409$	$-224\, 523$	$-38\, 304$
0, 4, 6	4091	4530	944
1, 3, 6	$-264\, 675$	$-179\, 100$	$-29\, 802$
2, 2, 6	$-589\, 124$	$-401\, 412$	$-65\, 012$
0, 5, 5	81 626	61 716	10 934
1, 4, 5	201 369	146 799	25 098
2, 3, 5	21 272	13 251	3662
2, 4, 4	630 060	423 324	69 660
3, 3, 4	867 568	570 030	94 504
0, 0, 9	7238	17 690	8712	920
0, 1, 8	$-175\, 166$	$-159\, 156$	$-51\, 444$	−6152
0, 2, 7	$-647\, 081$	$-652\, 818$	$-224\, 523$	$-25\, 536$
1, 1, 7	$-1499\, 799$	$-1530\, 330$	$-518\, 499$	$-57\, 770$
0, 3, 6	$-533\, 757$	$-529\, 350$	$-179\, 100$	$-19\, 868$
1, 2, 6	$-2705\, 027$	$-2822\, 469$	$-958\, 230$	$-104\, 548$
0, 4, 5	342 830	402 738	146 799	16 732
1, 3, 5	$-149\, 413$	$-184\, 542$	$-68\, 319$	−7260
2, 2, 5	$-1267\, 686$	$-1523\, 970$	$-548\, 802$	$-59\, 292$
1, 4, 4	2461 962	2533 506	845 712	91 340
2, 3, 4	5093 177	4833 423	1528 965	161 756
3, 3, 3	9375 357	8936 940	2832 975	298 590
0, 0, 8	$-281\, 388$	$-350\, 332$	$-159\, 156$	$-34\, 296$	−3076
0, 1, 7	$-2159\, 672$	$-2999\, 598$	$-1530\, 330$	$-345\, 666$	$-28\, 885$
0, 2, 6	$-3778\, 498$	$-5410\, 054$	$-2822\, 469$	$-638\, 820$	$-52\, 274$
1, 1, 6	$-8900\, 551$	$-13\, 027\, 620$	$-6879\, 357$	$-1558\, 260$	$-127\, 020$
0, 3, 5	$-135\, 055$	$-298\, 826$	$-184\, 542$	$-45\, 546$	−3630
1, 2, 5	$-4393\, 310$	$-7942\, 023$	$-4699\, 542$	$-1124\, 106$	$-92\, 600$
0, 4, 4	3530 318	4923 924	2533 506	563 808	45 670
1, 3, 4	16 104 342	19 972 701	9294 861	1917 870	148 130
2, 2, 4	24 563 089	28 940 436	12 769 305	2525 340	190 618
2, 3, 3	50 272 162	62 746 144	29 155 470	5991 840	459 710
0, 0, 7	$-2394\, 781$	$-4319\, 344$	$-2999\, 598$	$-1020\, 220$	$-172\, 833$	$-11\, 554$
0, 1, 6	$-9247\, 065$	$-17\, 801\, 102$	$-13\, 027\, 620$	$-4586\, 238$	$-779\, 130$	$-50\, 808$
0, 2, 5	$-2838\, 505$	$-8786\, 620$	$-7942\, 023$	$-3133\, 028$	$-562\, 053$	$-37\, 040$
1, 1, 5	$-8787\, 919$	$-25\, 549\, 530$	$-22\, 929\, 591$	$-9030\, 000$	$-1617\, 840$	$-106\, 764$
0, 3, 4	20 948 072	32 208 684	19 972 701	6196 574	958 935	59 252
1, 2, 4	69 096 062	99 234 474	56 262 246	15 871 764	2259 045	131 688
1, 3, 3	137 265 561	213 764 246	131 986 290	40 437 060	6154 785	372 498
2, 2, 3	234 255 941	365 678 846	225 361 089	68 713 090	10 396 950	625 360
0, 0, 6	$-7492\, 671$	$-18\, 494\, 130$	$-17\, 801\, 102$	$-8685\, 080$	$-2293\, 119$	$-311\, 652$
0, 1, 5	$-1486\, 702$	$-17\, 575\, 838$	$-25\, 549\, 530$	$-15\, 286\, 394$	$-4515\, 000$	$-647\, 136$
0, 2, 4	78 611 240	138 192 124	99 234 474	37 508 164	7935 882	903 618
1, 1, 4	167 093 949	288 887 940	201 077 370	72 579 090	14 512 785	1560 132
0, 3, 3	145 249 098	274 531 122	213 764 246	87 990 860	20 218 530	2461 914
1, 2, 3	536 224 445	1021 468 612	795 403 566	325 229 890	73 914 165	8881 236
2, 2, 2	925 883 649	1792 103 436	1416 380 508	586 686 276	134 672 895	16 279 056
0, 0, 5	5045 749	$-2973\, 404$	$-17\, 575\, 838$	$-17\, 033\, 020$	$-7643\, 197$	$-1806\, 000$
0, 1, 4	159 446 888	334 187 898	288 887 940	134 051 580	36 289 545	5805 114
0, 2, 3	471 598 432	1072 448 890	1021 468 612	530 269 044	162 614 945	29 565 666
1, 1, 3	1025 989 187	2345 852 082	2240 812 470	1163 278 640	355 858 545	64 434 762
1, 2, 2	1751 571 111	4078 346 544	3970 774 662	2101 403 324	654 668 325	120 379 620
0, 0, 4	127 800 660	318 893 776	334 187 898	192 591 960	67 025 790	14 515 818
0, 1, 3	763 219 155	2051 978 374	2345 852 082	1493 874 980	581 639 320	142 343 418
0, 2, 2	1281 194 946	3503 142 222	4078 346 544	2647 183 108	1050 701 662	261 867 330
1, 1, 2	2783 564 118	7679 836 074	9017 692 962	5898 612 100	2356 701 420	590 464 056
0, 0, 3	489 341 403	1526 438 310	2051 978 374	1563 901 388	746 937 490	232 655 728
0, 1, 2	1741 663 944	5567 128 236	7679 836 074	6011 795 308	2949 306 050	942 680 568
1, 1, 1	3756 967 623	12 122 069 778	16 883 573 100	13 343 681 760	6607 105 680	2130 074 520
0, 0, 2	950 909 571	3483 327 888	5567 128 236	5119 890 716	3005 897 654	1179 722 420
0, 1, 1	2032 807 482	7513 935 246	12 122 069 778	11 255 715 400	6671 840 880	2642 842 272
0, 0, 1	968 405 463	4065 614 964	7513 935 246	8081 379 852	5627 857 700	2668 736 352
0, 0, 0	408 831 948	1936 810 926	4065 614 964	5009 290 164	4040 689 926	2251 143 080

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{2}$
j	i  =  6	i  =  7	i  =  8	i  =  9	i  =  10	i  =  11	i  =  12
0, 0, 6	$-16\, 936$
0, 1, 5	$-35\, 588$
0, 2, 4	43 896
1, 1, 4	72 436
0, 3, 3	124 166
1, 2, 3	441 476
2, 2, 2	810 192
0, 0, 5	$-215\, 712$	$-10\, 168$
0, 1, 4	520 044	20 696
0, 2, 3	2960 412	126 136
1, 1, 3	6419 322	271 980
1, 2, 2	12 120 276	515 776
0, 0, 4	1935 038	148 584	5174
0, 1, 3	21 478 254	1834 092	67 995
0, 2, 2	40 126 540	3462 936	128 944
1, 1, 2	90 813 954	7851 312	292 230
0, 0, 3	47 447 806	6136 644	458 523	15 110
0, 1, 2	196 821 352	25 946 844	1962 828	64 940
1, 1, 1	448 120 638	59 430 600	4512 915	149 490
0, 0, 2	314 226 856	56 234 672	6486 711	436 184	12 988
0, 1, 1	710 024 840	128 034 468	14 857 650	1002 870	29 898
0, 0, 1	880 947 424	202 864 240	32 008 617	3301 700	200 574	5436
0, 0, 0	889 578 784	251 699 264	50 716 060	7113 026	660 340	36 468	906

Table B3. The polynomial coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{3}$ of the seventh-order recurrence scheme for generic points on a simple cubic 3D lattice. Here, j is the multi-index of r_j in (16).

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{3}$
j	i  =  0	i  =  1	i  =  2	i  =  3
0, 0, 12	1024
0, 1, 11	1920
0, 2, 10	−1464
1, 1, 10	3984
0, 3, 9	−2294
1, 2, 9	−1122
0, 4, 8	3996
1, 3, 8	−5976
2, 2, 8	−1512
0, 5, 7	2070
1, 4, 7	−252
2, 3, 7	−8802
0, 6, 6	−3720
1, 5, 6	6534
2, 4, 6	−2952
3, 3, 6	−5676
2, 5, 5	7992
3, 4, 5	5796
4, 4, 4	9000
0, 0, 11	18 496	3840
0, 1, 10	41 864	7968
0, 2, 9	−908	−2244
1, 1, 9	53 672	7032
0, 3, 8	$-50\, 910$	$-11\, 952$
1, 2, 8	$-111\, 726$	$-29\, 988$
0, 4, 7	−3198	−504
1, 3, 7	$-191\, 595$	$-43\, 425$
2, 2, 7	$-348\, 570$	$-81\, 234$
0, 5, 6	65 280	13 068
1, 4, 6	51 306	9504
2, 3, 6	$-254\, 202$	$-55\, 116$
1, 5, 5	238 734	47 754
2, 4, 5	206 334	51 102
3, 3, 5	102 048	32 832
3, 4, 4	431 130	112 950
0, 0, 10	161 392	83 728	7968
0, 1, 9	241 654	107 344	7032
0, 2, 8	$-379\, 206$	$-223\, 452$	$-29\, 988$
1, 1, 8	$-685\, 152$	$-364\, 896$	$-54\, 108$
0, 3, 7	$-760\, 695$	$-383\, 190$	$-43\, 425$
1, 2, 7	$-3341\, 493$	$-1634\, 841$	$-190\, 890$
0, 4, 6	201 738	102 612	9504
1, 3, 6	$-2432\, 364$	$-1170\, 582$	$-131\, 292$
2, 2, 6	$-5827\, 212$	$-2709\, 828$	$-299\, 880$
0, 5, 5	1005 930	477 468	47 754
1, 4, 5	2439 423	1142 091	118 530
2, 3, 5	59 661	92 475	24 102
2, 4, 4	5849 172	2781 828	326 160
3, 3, 4	7438 434	3666 786	444 960
0, 0, 9	440 892	483 308	107 344	4688
0, 1, 8	$-2140\, 692$	$-1370\, 304$	$-364\, 896$	$-36\, 072$
0, 2, 7	$-8884\, 005$	$-6682\, 986$	$-1634\, 841$	$-127\, 260$
1, 1, 7	$-21\, 112\, 605$	$-15\, 484\, 788$	$-3754\, 953$	$-294\, 930$
0, 3, 6	$-6429\, 528$	$-4864\, 728$	$-1170\, 582$	$-87\, 528$
1, 2, 6	$-37\, 583\, 655$	$-27\, 920\, 673$	$-6604\, 272$	$-489\, 528$
0, 4, 5	6178 569	4878 846	1142 091	79 020
1, 3, 5	$-2461\, 971$	$-1632\, 150$	$-334\, 863$	$-17\, 460$
2, 2, 5	$-23\, 729\, 166$	$-17\, 479\, 908$	$-3900\, 366$	$-254\, 772$
1, 4, 4	32 299 542	24 348 078	5758 704	426 420
2, 3, 4	52 467 522	38 973 339	9546 381	766 260
3, 3, 3	97 608 483	72 326 514	17 640 105	1405 830
0, 0, 8	$-5443\, 524$	$-4281\, 384$	$-1370\, 304$	$-243\, 264$
0, 1, 7	$-42\, 915\, 585$	$-42\, 225\, 210$	$-15\, 484\, 788$	$-2503\, 302$
0, 2, 6	$-73\, 421\, 307$	$-75\, 167\, 310$	$-27\, 920\, 673$	$-4402\, 848$
1, 1, 6	$-182\, 021\, 079$	$-186\, 924\, 762$	$-69\, 167\, 835$	$-10\, 867\, 380$
0, 3, 5	$-5391\, 378$	$-4923\, 942$	$-1632\, 150$	$-223\, 242$
1, 2, 5	$-134\, 887\, 941$	$-142\, 928\, 073$	$-52\, 797\, 690$	$-7985\, 682$
0, 4, 4	61 115 736	64 599 084	24 348 078	3839 136
1, 3, 4	192 861 660	194 160 483	72 289 413	11 759 190
2, 2, 4	238 468 221	232 046 280	86 206 005	14 493 924
2, 3, 3	596 405 664	595 463 646	220 756 644	35 933 520
0, 0, 7	$-67\, 590\, 699$	$-85\, 831\, 170$	$-42\, 225\, 210$	$-10\, 323\, 192$
0, 1, 6	$-270\, 863\, 238$	$-364\, 042\, 158$	$-186\, 924\, 762$	$-46\, 111\, 890$
0, 2, 5	$-187\, 727\, 322$	$-269\, 775\, 882$	$-142\, 928\, 073$	$-35\, 198\, 460$
1, 1, 5	$-537\, 982\, 365$	$-776\, 979\, 378$	$-414\, 012\, 393$	$-102\, 331\, 080$
0, 3, 4	305 304 120	385 723 320	194 160 483	48 192 942
1, 2, 4	696 155 145	801 836 712	381 333 708	94 465 332
1, 3, 3	1913 853 333	2405 206 230	1200 740 040	296 860 260
2, 2, 3	3206 144 283	4020 693 408	2001 521 601	494 486 406
0, 0, 6	$-316\, 850\, 521$	$-541\, 726\, 476$	$-364\, 042\, 158$	$-124\, 616\, 508$
0, 1, 5	$-558\, 893\, 574$	$-1075\, 964\, 730$	$-776\, 979\, 378$	$-276\, 008\, 262$
0, 2, 4	1026 348 966	1392 310 290	801 836 712	254 222 472
1, 1, 4	1855 533 075	2345 106 150	1248 340 176	374 346 450
0, 3, 3	2527 042 978	3827 706 666	2405 206 230	800 493 360
1, 2, 3	8942 209 551	13 555 120 386	8483 646 240	2812 115 910
2, 2, 2	16 013 753 613	24 683 704 668	15 629 922 918	5211 040 284
0, 0, 5	$-425\, 028\, 007$	$-1117\, 787\, 148$	$-1075\, 964\, 730$	$-517\, 986\, 252$
0, 1, 4	2460 861 643	3711 066 150	2345 106 150	832 226 784
0, 2, 3	9977 530 949	17 884 419 102	13 555 120 386	5655 764 160
1, 1, 3	21 399 417 097	38 490 999 438	29 182 241 160	12 156 880 020
1, 2, 2	38 346 916 893	70 722 151 488	54 814 248 306	23 218 459 020
0, 0, 4	2732 264 750	4921 723 286	3711 066 150	1563 404 100
0, 1, 3	20 336 507 375	42 798 834 194	38 490 999 438	19 454 827 440
0, 2, 2	35 515 908 870	76 693 833 786	70 722 151 488	36 542 832 204
1, 1, 2	77 911 117 200	170 216 315 664	158 513 197 158	82 518 255 840
0, 0, 3	16 659 002 103	40 673 014 750	42 798 834 194	25 660 666 292
0, 1, 2	61 471 175 688	155 822 234 400	170 216 315 664	105 675 464 772
1, 1, 1	134 436 520 755	345 447 658 950	382 252 444 410	240 044 395 440
0, 0, 2	42 102 546 939	122 942 351 376	155 822 234 400	113 477 543 776
0, 1, 1	90 818 056 458	268 873 041 510	345 447 658 950	254 834 962 940
0, 0, 1	53 487 879 255	181 636 112 916	268 873 041 510	230 298 439 300
0, 0, 0	27 487 989 900	106 975 758 510	181 636 112 916	179 248 694 340

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{3}$
j	i  =  4	i  =  5	i  =  6	i  =  7
0, 0, 8	$-18\, 036$
0, 1, 7	$-147\, 465$
0, 2, 6	$-244\, 764$
1, 1, 6	$-603\, 000$
0, 3, 5	−8730
1, 2, 5	$-406\, 260$
0, 4, 4	213 210
1, 3, 4	704 430
2, 2, 4	923 490
2, 3, 3	2156 850
0, 0, 7	$-1251\, 651$	$-58\, 986$
0, 1, 6	$-5433\, 690$	$-241\, 200$
0, 2, 5	$-3992\, 841$	$-162\, 504$
1, 1, 5	$-11\, 643\, 660$	$-477\, 036$
0, 3, 4	5879 595	281 772
1, 2, 4	12 194 505	648 432
1, 3, 3	36 260 175	1744 722
2, 2, 3	60 531 030	2925 864
0, 0, 6	$-23\, 055\, 945$	$-2173\, 476$	$-80\, 400$
0, 1, 5	$-51\, 165\, 540$	$-4657\, 464$	$-159\, 012$
0, 2, 4	47 232 666	4877 802	216 144
1, 1, 4	69 408 585	7577 892	366 660
0, 3, 3	148 430 130	14 504 070	581 574
1, 2, 3	520 972 965	51 077 124	2063 124
2, 2, 2	965 029 653	93 991 968	3749 976
0, 0, 5	$-138\, 004\, 131$	$-20\, 466\, 216$	$-1552\, 488$	$-45\, 432$
0, 1, 4	187 173 225	27 763 434	2525 964	104 760
0, 2, 3	1406 057 955	208 389 186	17 025 708	589 464
1, 1, 3	3017 433 825	447 051 402	36 573 054	1269 900
1, 2, 2	5815 038 825	861 343 740	69 779 556	2378 304
0, 0, 4	416 113 392	74 869 290	9254 478	721 704
0, 1, 3	6078 440 010	1206 973 530	149 017 134	10 449 444
0, 2, 2	11 609 229 510	2326 015 530	287 114 580	19 937 016
1, 1, 2	26 344 042 290	5290 777 044	652 954 554	45 217 584
0, 0, 3	9727 413 720	2431 376 004	402 324 510	42 576 324
0, 1, 2	41 259 127 920	10 537 616 916	1763 592 348	186 558 444
1, 1, 1	94 595 142 600	24 318 722 736	4083 803 514	431 946 360
0, 0, 2	52 837 732 386	16 503 651 168	3512 538 972	503 883 528
0, 1, 1	120 022 197 720	37 838 057 040	8106 240 912	1166 801 004
0, 0, 1	127 417 481 470	48 008 879 088	12 612 685 680	2316 068 832
0, 0, 0	115 149 219 650	50 966 992 588	16 002 959 696	3603 624 480

Coefficients $\newcommand{\qCoef}[1]{q_{#1,i;j}} \qCoef{3}$
j	i  =  8	i  =  9	i  =  10	i  =  11	i  =  12
0, 0, 4	26 190
0, 1, 3	317 475
0, 2, 2	594 576
1, 1, 2	1341 630
0, 0, 3	2612 361	70 550
0, 1, 2	11 304 396	298 140
1, 1, 1	26 076 105	682 770
0, 0, 2	46 639 611	2512 088	59 628
0, 1, 1	107 986 590	5794 690	136 554
0, 0, 1	291 700 251	23 997 020	1158 938	24 828
0, 0, 0	579 017 208	64 822 278	4799 404	210 716	4138