The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
28,491 characters
Positive weight Hermite and Legendre quadrature rules
\title{Positive weight Hermite and Legendre quadrature rules}
\author{Joris Pinkse\thanks{\href{[email removed]}{[email removed]}}\\
Department of Economics\\
Penn State}
\date{September 2026}
\maketitle
\thispagestyle{empty}
\begin{abstract}
This paper contains new 80-digit positive-weight Gauss-Hermite and positive-weight interior-node Gauss-Legendre quadrature rules for up to five dimensions and varying polynomial degree accuracy (depending on quadrature type and dimension). Some of these rules improve on the best available rules in the literature and some offer rules where none (other than the tensor product) existed. The results were produced by combining methodology developed by previous researchers with new approaches. The full precision rules themselves are at \url{https://doi.org/10.5281/zenodo.22881864}. Software in Julia, Python, and R providing the rules is available via a package registry and/or github: \href{https://github.com/NittanyLion/Quadriceps.jl}{Quadriceps.jl} (both 64-bit and 128-bit), \href{https://github.com/NittanyLion/quadriceps-py}{quadriceps-py} (64-bit only), and \href{https://github.com/NittanyLion/quadriceps-r}{quadriceps-r} (64-bit only). The software used to create these rules is available via \href{https://github.com/NittanyLion/PositiveWeightQuadratureSolvers.jl}{PositiveWeightQuadratureSolvers.jl}. The replication package is at \href{https://github.com/NittanyLion/QuadricepsReplicationPackage.jl}{QuadricepsReplicationPackage.jl}. A snapshot of all five packages is archived at \url{https://doi.org/10.5281/zenodo.22883240}.
\end{abstract}
\clearpage
\onehalfspacing
\section{Introduction}
This paper contains new 80-digit positive-weight Gauss-Hermite and positive-weight interior-node Gauss-Legendre quadrature rules for up to five dimensions and varying polynomial degree accuracy (depending on quadrature type and dimension). Some of these rules improve on the best available rules in the literature and some offer rules where none (other than the tensor product) existed. The results were produced by combining methodology developed by previous researchers with new approaches. The full precision rules themselves are at \url{https://doi.org/10.5281/zenodo.22881864}. Software in Julia, Python, and R providing the rules is available via a package registry and/or github: \href{https://github.com/NittanyLion/Quadriceps.jl}{Quadriceps.jl} (both 64-bit and 128-bit), \href{https://github.com/NittanyLion/quadriceps-py}{quadriceps-py} (64-bit only), and \href{https://github.com/NittanyLion/quadriceps-r}{quadriceps-r} (64-bit only). The software used to create these rules is available via \href{https://github.com/NittanyLion/PositiveWeightQuadratureSolvers.jl}{PositiveWeightQuadratureSolvers.jl}. The replication package is at \href{https://github.com/NittanyLion/QuadricepsReplicationPackage.jl}{QuadricepsReplicationPackage.jl}. A snapshot of all five packages is archived at \url{https://doi.org/10.5281/zenodo.22883240}.\footnote{The Python and R packages are 64-bit only: neither language has a standard floating point type beyond double precision. Julia has none either, but its generic typing lets Quadriceps.jl return a rule in whatever number type the caller asks for, \texttt{Float128} through Quadmath.jl or \texttt{BigFloat} from Base.}
The objective, for given $f,d$ and odd $p$, is to construct a rule of $N$ (node, weight) pairs $(n,w)$ for which
\begin{equation} \label{eq:integral equality}
∫ ∏_{j=1}ᵈ xⱼ^{aⱼ} f(x) \mathrm{d}x = ∑_{k=1}^N wₖ ∏_{j=1}ᵈ n_{kj}^{aⱼ}
\end{equation}
for every vector $a$ of integers with $aⱼ ≥ 0$ and $|a| = ∑_{j=1}ᵈ aⱼ ≤ p$. These `rules' can then be used to approximate $∫ g(x) f(x) \mathrm{d}x$ by $∑_{k=1}^N wₖ g(nₖ)$. The results in this paper apply to the case in which $f$ is one of the standard normal (Gauss-Hermite) and standard uniform (Gauss-Legendre) density functions.\footnote{The Gauss-Legendre domain is $[0,1]ᵈ$ in this paper.}
I focus on the case of positive weights such that the integral approximation can be used as an argument to a logarithmic function. In particular, this paper is motivated by random coefficients logit models as they are used in economics where the integrals are choice probabilities that enter a loglikelihood function.
Positive weights matter beyond ensuring that the approximation lies between the minimum and maximum function values of the function being integrated. Positive weights ensure that there is no loss of digits due to cancellation of positive and negative weight contributions. It has been shown \citep[see e.g.][section 4.1]{glaubitz21} that in a cube, the error of a quadrature rule is at most $1+∑ₖ|wₖ|$ times the error of the best uniform approximation by a polynomial of the same degree: only if the weights are nonnegative is $∑ₖ |wₖ| = ∑ₖ wₖ = 1$. For the sparse grids that \citet{heiss08} provide for the cube, $∑ₖ |wₖ|$ is as much as 34 ($d=5$, $p=21$).
The easiest way of finding $(n,w)$ satisfying \cref{eq:integral equality} is to solve the problem in one dimension and then take the tensor product of the solution. Doing so ensures positive weights, but it also requires $N = q^d$ nodes where $q=(p+1)/2$, i.e.\ the number of nodes then grows exponentially in $d$.
Except in very low $p$ cases, there likely exist solutions to \cref{eq:integral equality} that require fewer nodes since the tensor product imposes more than just the restrictions in \cref{eq:integral equality}. Indeed, if for instance $d=2,p=7$ then for the tensor product rule, \cref{eq:integral equality} also holds exactly when $a₁=a₂ = 6$, although $a₁+a₂ =12>7 = p$. The number of nodes can be reduced significantly by letting go of cross moments of order greater than $p$. The hardest part of this exercise is to achieve a low node count while imposing the positive-weight constraint.
None of the above is to say that the extra constraints provided by the tensor product have no value. Indeed, one would expect $∫ g(x) f(x) \mathrm{d}x$ to be approximated better by a tensor product-based approximation than by a reduced node quadrature rule of the same order. But there is no reason to believe that the tensor product provides the optimal approximation for a given node count.
I also make no claim that a plain Gauss-Hermite or Gauss-Legendre rule is necessarily the optimal way to do numerical integration. Accuracy can often be improved by transforming the integrand before a rule is applied, for instance by recentering and rescaling it around its mode, as in adaptive Gauss-Hermite quadrature \citep{naylor82,liu94}, or by splitting a bounded domain into subregions and applying a rule to each \citep{genz80}. But such improvements are not the focus of this paper: they can be used in conjunction with the rules provided here, also.
I am not the first to study this problem. Indeed, I will be making use of, combining, and building on methodology developed by others. A complete description of the methodology used, including references, can be found in \cref{app:methods}. The methodology used combines the following strategies. It attempts to find a (nodes, weights) pair that is a solution to moment equations of (an orthogonalized version of) the form \cref{eq:integral equality}. It eliminates nodes by exploiting symmetry, by trying to merge nodes, and by trying to find a new solution by reoptimizing after dropping nodes with small weights. I further tried imposing symmetry: imposing symmetry lightens the computational load, but can increase the minimal achievable node count.
I employed various warm-start strategies. First, I used a laddering approach by starting for a fixed dimension $d$ from a lower value of $p$ and adding nodes from there. I further tried starting from the tensor product of a $(d₁,p)$ and a $(d₂,p)$ solution to find a good $(d₁+d₂,p)$ solution. I tried taking a monotonic transformation from the Gauss-Hermite solution to obtain a starting point for the Gauss-Legendre solution and vice versa. Finally, I used existing rules as a warm start and improved on them from there.
Any rules found were polished using a high-precision optimizer.
The computations were conducted over a six-week period on a combination of the Open Science Pool, Penn State's ROAR cluster, my office and home desktops, and my work and personal laptops.
\begin{figure}[ht]
\centering
\begin{tikzpicture}
\begin{groupplot}[group style={group size=2 by 1, horizontal sep=2.3cm}, quad, scale only axis, width=0.34\linewidth, height=0.34\linewidth,
grid=none, colormap={qblue}{color=(qd2!45) color=(qd3) color=(qd5)}, colorbar,
colorbar style={width=2mm, tick label style={font=\scriptsize}, title={$\log_{10} w$}, title style={font=\scriptsize}, axis line style={gray!60}}]
\nextgroupplot[title={normal, $p=33$, $N=208$}, xmin=-8.2, xmax=8.2, ymin=-8.2, ymax=8.2]
\addplot[scatter, only marks, mark=*, mark size=1.3pt, scatter src=explicit, mark options={draw=gray!70, line width=0.15pt}] table[x=x, y=y, meta=logw] {fig/nodes_gh.dat};
\nextgroupplot[title={uniform, $p=77$, $N=1032$}, xmin=0, xmax=1, ymin=0, ymax=1]
\addplot[scatter, only marks, mark=*, mark size=0.9pt, scatter src=explicit, mark options={draw=gray!70, line width=0.1pt}] table[x=x, y=y, meta=logw] {fig/nodes_le.dat};
\end{groupplot}
\end{tikzpicture}
\caption{The largest planar rule of each table: nodes, shaded by the logarithm of their weight.}
\label{fig:nodes}
\end{figure}
As will become clear in \cref{sec:results}, the tables in this paper improve on the rules available in the literature for many cases and provide new ones for many cases for which no published rules (other than the tensor product) exist. For Gauss-Hermite, I develop entirely new or improved rules in 34 of 57 cases whereas for Gauss-Legendre the count is 62 of 85.
Although my primary objective was to improve the existing Gauss-Hermite results, my Gauss-Legendre rules cover more cases. This is simply due to the fact that finding a Gauss-Hermite solution is harder because of the presence of (node, weight) combinations that are extreme in their location (nodes) and magnitude (weight). This difference in the spread of the weights between Gauss-Hermite and Gauss-Legendre solutions is illustrated in \cref{fig:nodes}.
\section{Results}
\label{sec:results}
\begin{table}[p]
\centering
\setlength{\tabcolsep}{3.5pt}
\small
\begin{tabular}[t]{rrrrrrrc}
\toprule
$p$ & $N$ & $\rho$ & rel.\ err. & M\"o & sym & prev & src \\
\midrule
\multicolumn{8}{@{}l}{$d = 2$} \\
\rowcolor{gray!40}1 & 1 & 1.00 & 0 & 1 & 1 & 1 & --- \\
\rowcolor{gray!40}3 & 4 & 1.00 & 3.6e-70 & 4 & 4 & 4 & CO \\
\rowcolor{gray!40}5 & 7 & 0.88 & 2.7e-69 & 7 & 7 & 7 & CO \\
\rowcolor{gray!40}7 & 12 & 0.87 & 1.2e-69 & 12 & 12 & 12 & CO \\
\rowcolor{gray!15}9 & 18 & 0.85 & 6.8e-69 & 17 & 17 & 18 & HP \\
\rowcolor{gray!15}11 & 25 & 0.83 & 1.1e-79 & 24 & 25 & 25 & HP \\
\rowcolor{gray!15}13 & 34 & 0.83 & 8.9e-80 & 31 & 33 & 34 & CH \\
\rowcolor{gray!15}15 & 44 & 0.83 & 1.3e-79 & 40 & 44 & 44 & CO \\
17 & 55 & 0.82 & 6.4e-79 & 49 & 55 & 57 & SC \\
19 & 68 & 0.82 & 6.8e-78 & 60 & 68 & 71 & CO \\
21 & 82 & 0.82 & 3.0e-77 & 71 & 81 & 90 & CO \\
\rowcolor{gray!15}23 & 97 & 0.82 & 8.7e-74 & 84 & 97 & 97 & CO \\
25 & 114 & 0.82 & 6.5e-71 & 97 & 113 & 127 & CO \\
27 & 132 & 0.82 & 1.1e-69 & 112 & 132 & --- & --- \\
29 & 153 & 0.82 & 5.0e-78 & 127 & 151 & --- & --- \\
31 & 178 & 0.83 & 7.9e-69 & 144 & 172 & --- & --- \\
33 & 208 & 0.85 & 7.2e-69 & 161 & 193 & --- & --- \\
\midrule
\multicolumn{8}{@{}l}{$d = 3$} \\
\rowcolor{gray!40}1 & 1 & 1.00 & 0 & 1 & 1 & 1 & --- \\
\rowcolor{gray!40}3 & 6 & 0.91 & 5.4e-71 & 6 & 6 & 6 & CO \\
\rowcolor{gray!40}5 & 13 & 0.78 & 1.1e-69 & 13 & 13 & 13 & ST \\
\rowcolor{gray!15}7 & 27 & 0.75 & 6.0e-69 & 26 & 27 & 27 & CO \\
\rowcolor{gray!15}9 & 45 & 0.71 & 9.3e-69 & 43 & 45 & 45 & KO \\
\rowcolor{gray!15}11 & 77 & 0.71 & 1.0e-79 & 68 & 77 & 77 & CO \\
13 & 128 & 0.72 & 6.3e-80 & 99 & 119 & 137 & CO \\
15 & 184 & 0.71 & 5.2e-79 & 140 & 175 & --- & --- \\
17 & 264 & 0.71 & 2.6e-78 & 189 & 248 & --- & --- \\
19 & 354 & 0.71 & 1.4e-69 & 250 & 342 & --- & --- \\
21 & 476 & 0.71 & 1.8e-71 & 321 & 447 & --- & --- \\
23 & 597 & 0.70 & 3.0e-72 & 406 & 563 & --- & --- \\
25 & 776 & 0.71 & 5.1e-73 & 503 & 720 & --- & --- \\
27 & 966 & 0.71 & 2.6e-71 & 616 & 894 & --- & --- \\
29 & 1242 & 0.72 & 6.3e-72 & 743 & 1094 & --- & --- \\
31 & 1848 & 0.77 & 2.1e-71 & 888 & 1328 & --- & --- \\
33 & 2226 & 0.77 & 1.1e-69 & 1049 & 1588 & --- & --- \\
\bottomrule
\end{tabular}
\hfill
\begin{tabular}[t]{rrrrrrrc}
\toprule
$p$ & $N$ & $\rho$ & rel.\ err. & M\"o & sym & prev & src \\
\midrule
\multicolumn{8}{@{}l}{$d = 4$} \\
\rowcolor{gray!40}1 & 1 & 1.00 & 0 & 1 & 1 & 1 & --- \\
\rowcolor{gray!40}3 & 8 & 0.84 & 6.6e-70 & 8 & 8 & 8 & CO \\
\rowcolor{gray!15}5 & 22 & 0.72 & 1.8e-69 & 21 & 22 & 22 & CO \\
\rowcolor{gray!15}7 & 49 & 0.66 & 7.2e-69 & 48 & 49 & 49 & CO \\
9 & 116 & 0.66 & 7.7e-80 & 91 & 105 & --- & --- \\
11 & 193 & 0.62 & 8.3e-80 & 160 & 192 & --- & --- \\
13 & 414 & 0.64 & 6.7e-80 & 259 & 353 & --- & --- \\
15 & 577 & 0.61 & 7.2e-80 & 400 & 529 & --- & --- \\
17 & 1056 & 0.63 & 1.1e-77 & 589 & 856 & --- & --- \\
19 & 1505 & 0.62 & 1.5e-79 & 840 & 1281 & --- & --- \\
21 & 2318 & 0.63 & 1.9e-69 & 1161 & 1809 & --- & --- \\
23 & 3238 & 0.63 & 1.3e-71 & 1568 & 2512 & --- & --- \\
\midrule
\multicolumn{8}{@{}l}{$d = 5$} \\
\rowcolor{gray!40}1 & 1 & 1.00 & 0 & 1 & 1 & 1 & --- \\
\rowcolor{gray!40}3 & 10 & 0.79 & 2.2e-70 & 10 & 10 & 10 & CO \\
\rowcolor{gray!15}5 & 32 & 0.67 & 2.1e-69 & 31 & 32 & 32 & SS \\
\rowcolor{gray!15}7 & 83 & 0.61 & 8.2e-80 & 80 & 83 & 83 & ST \\
9 & 244 & 0.60 & 5.4e-80 & 171 & 212 & 395 & AD \\
11 & 485 & 0.57 & 9.7e-80 & 332 & 445 & --- & --- \\
13 & 1135 & 0.58 & 1.9e-79 & 591 & 846 & --- & --- \\
15 & 1767 & 0.56 & 5.5e-80 & 992 & 1496 & --- & --- \\
17 & 3986 & 0.58 & 1.4e-79 & 1581 & 2448 & --- & --- \\
19 & 7174 & 0.59 & 3.0e-79 & 2422 & 3978 & --- & --- \\
21 & 13199 & 0.61 & 9.8e-80 & 3583 & 6290 & --- & --- \\
\bottomrule
\end{tabular}
\normalsize
\caption{Positive-weight rules for the Gaussian weight $N(0, I_d)$. $N$: nodes in my best rule; $\rho = N^{1/d}/q$ with $q = (p+1)/2$ ($\rho = 1$ is the Gauss product grid, smaller is better); rel.\ err.: largest relative monomial error, $\max_{|a| \le p} |\sum_s w_s x_s^a - \mathrm{E}\,x^a| / \max(\sum_s w_s\,|x_s^a|, 1)$, of the rule's 80-digit file (rounded to double precision: at most $1.1\cdot 10^{-15}$); M\"o: M\"oller's lower bound; sym: the symmetric counting floor, the fewest nodes for which a rule invariant under a group $G$ has as many unknowns as it has $G$-invariant moment conditions, overall and among the conditions that only some orbit types see, minimized over the orbit types and over $G$ (no symmetry, $\pm$pairs, a sign change in each coordinate, the symmetry group $B_d$ of the cube and, for $d = 2$, the dihedral and cyclic groups), not counting the $d(d-1)/2$ unknowns that only rotate the rule where $G$ commutes with rotations --- a parameter count, not a proven bound; prev: best count in the literature, --- where none exists besides the product grid; src: where that count is recorded --- the original paper where I could identify it, otherwise a compilation (AD = \citealt{adurthi12}; CH = \citealt{coolio88b}; CO = \citealt{coolio}; HP = \citealt{haggy76,haggy77}; KO = \citealt{konyaev77}; SC = \citealt{degani05}; SS = \citealt{stroud63}; ST = \citealt{stroud71}). Light gray: no improvement on prev; dark gray: M\"oller's bound attained (proven minimal).}\label{tab:best-gh}
\end{table}
\begin{table}[p]
\centering
\setlength{\tabcolsep}{3.5pt}
\footnotesize
\begin{tabular}[t]{rrrrrrrc}
\toprule
$p$ & $N$ & $\rho$ & rel.\ err. & M\"o & pair & prev & src \\
\midrule
\multicolumn{8}{@{}l}{$d = 2$} \\
\rowcolor{gray!40}1 & 1 & 1.00 & 0 & 1 & 1 & 1 & --- \\
\rowcolor{gray!40}3 & 4 & 1.00 & 5.2e-75 & 4 & 3 & 4 & ST \\
\rowcolor{gray!40}5 & 7 & 0.88 & 1.4e-76 & 7 & 6 & 7 & ST \\
\rowcolor{gray!40}7 & 12 & 0.87 & 1.1e-74 & 12 & 11 & 12 & ST \\
\rowcolor{gray!40}9 & 17 & 0.82 & 1.4e-78 & 17 & 17 & 17 & MO \\
\rowcolor{gray!40}11 & 24 & 0.82 & 4.7e-73 & 24 & 24 & 24 & CH \\
\rowcolor{gray!15}13 & 33 & 0.82 & 4.1e-72 & 31 & 33 & 33 & CH \\
\rowcolor{gray!15}15 & 43 & 0.82 & 6.7e-74 & 40 & 43 & 43 & FS \\
\rowcolor{gray!15}17 & 54 & 0.82 & 1.5e-78 & 49 & 54 & 54 & FS \\
\rowcolor{gray!15}19 & 67 & 0.82 & 1.4e-71 & 60 & 67 & 67 & FS \\
\rowcolor{gray!15}21 & 81 & 0.82 & 1.4e-78 & 71 & 81 & 81 & FS \\
\rowcolor{gray!15}23 & 96 & 0.82 & 2.5e-78 & 84 & 96 & 96 & FS \\
\rowcolor{gray!15}25 & 113 & 0.82 & 1.9e-70 & 97 & 113 & 113 & FS \\
27 & 131 & 0.82 & 1.2e-70 & 112 & 131 & 132 & FS \\
29 & 150 & 0.82 & 1.3e-78 & 127 & 150 & 152 & XG \\
31 & 171 & 0.82 & 2.7e-72 & 144 & 171 & 172 & FS \\
33 & 194 & 0.82 & 2.7e-73 & 161 & 193 & 197 & FS \\
35 & 216 & 0.82 & 1.9e-77 & 180 & 216 & 220 & FS \\
37 & 242 & 0.82 & 5.3e-76 & 199 & 241 & 245 & FS \\
39 & 268 & 0.82 & 3.1e-73 & 220 & 267 & 274 & FS \\
41 & 294 & 0.82 & 3.8e-74 & 241 & 294 & 303 & FS \\
43 & 324 & 0.82 & 8.6e-73 & 264 & 323 & 331 & FS \\
45 & 356 & 0.82 & 8.9e-72 & 287 & 353 & 359 & FS \\
47 & 388 & 0.82 & 2.1e-71 & 312 & 384 & 396 & FS \\
49 & 422 & 0.82 & 2.5e-78 & 337 & 417 & 427 & FS \\
51 & 456 & 0.82 & 8.1e-73 & 364 & 451 & 462 & FS \\
53 & 492 & 0.82 & 7.1e-75 & 391 & 486 & 498 & FS \\
55 & 528 & 0.82 & 2.7e-70 & 420 & 523 & 536 & FS \\
57 & 570 & 0.82 & 2.0e-71 & 449 & 561 & 576 & DW \\
59 & 606 & 0.82 & 1.2e-72 & 480 & 600 & 613 & DW \\
61 & 642 & 0.82 & 2.0e-78 & 511 & 641 & 660 & DW \\
63 & 692 & 0.82 & 1.1e-74 & 544 & 683 & 709 & DW \\
65 & 732 & 0.82 & 5.7e-71 & 577 & 726 & 757 & DW \\
67 & 780 & 0.82 & 2.0e-73 & 612 & 771 & 805 & DW \\
69 & 826 & 0.82 & 1.0e-70 & 647 & 817 & 849 & DW \\
71 & 872 & 0.82 & 2.1e-73 & 684 & 864 & 904 & DW \\
73 & 922 & 0.82 & 1.5e-70 & 721 & 913 & 953 & DW \\
75 & 980 & 0.82 & 2.1e-70 & 760 & 963 & 1001 & DW \\
77 & 1032 & 0.82 & 6.9e-74 & 799 & 1014 & 1049 & DW* \\
\midrule
\multicolumn{8}{@{}l}{$d = 3$} \\
\rowcolor{gray!40}1 & 1 & 1.00 & 0 & 1 & 1 & 1 & --- \\
\rowcolor{gray!40}3 & 6 & 0.91 & 1.1e-76 & 6 & 4 & 6 & ST \\
\rowcolor{gray!40}5 & 13 & 0.78 & 4.9e-76 & 13 & 12 & 13 & ST \\
\rowcolor{gray!40}7 & 26 & 0.74 & 6.9e-76 & 26 & 26 & 26 & XG \\
\bottomrule
\end{tabular}
\hfill
\begin{tabular}[t]{rrrrrrrc}
\toprule
$p$ & $N$ & $\rho$ & rel.\ err. & M\"o & pair & prev & src \\
\midrule
\multicolumn{8}{@{}l}{$d = 3$ (cont.)} \\
9 & 48 & 0.73 & 6.9e-75 & 43 & 48 & 50 & XG \\
11 & 82 & 0.72 & 1.9e-74 & 68 & 81 & 84 & XG \\
13 & 128 & 0.72 & 5.7e-74 & 99 & 126 & 130 & XG \\
15 & 188 & 0.72 & 6.9e-74 & 140 & 186 & 190 & XG \\
17 & 266 & 0.71 & 4.2e-73 & 189 & 263 & 282 & DW \\
19 & 360 & 0.71 & 9.6e-73 & 250 & 358 & 369 & DW \\
21 & 476 & 0.71 & 1.5e-71 & 321 & 474 & 495 & DW \\
23 & 612 & 0.71 & 7.2e-70 & 406 & 612 & 617 & DW \\
25 & 776 & 0.71 & 2.2e-71 & 503 & 774 & 828 & DW \\
27 & 964 & 0.71 & 2.2e-75 & 616 & 963 & 984 & DW* \\
29 & 1184 & 0.71 & 9.3e-78 & 743 & 1180 & 1258 & DW* \\
31 & 1432 & 0.70 & 2.8e-76 & 888 & 1428 & 1478 & DW* \\
33 & 1714 & 0.70 & 6.7e-78 & 1049 & 1709 & 1787 & DW* \\
35 & 2028 & 0.70 & 2.6e-76 & 1230 & 2024 & 2102 & DW* \\
37 & 2380 & 0.70 & 1.0e-75 & 1429 & 2376 & 2506 & DW* \\
39 & 2770 & 0.70 & 2.7e-74 & 1650 & 2766 & 2856 & DW* \\
41 & 3200 & 0.70 & 7.9e-75 & 1891 & 3196 & 3338 & DW* \\
43 & 3704 & 0.70 & 4.9e-71 & 2156 & 3669 & 3870 & DW* \\
45 & 4308 & 0.71 & 2.3e-71 & 2443 & 4186 & 4414 & DW* \\
\midrule
\multicolumn{8}{@{}l}{$d = 4$} \\
\rowcolor{gray!40}1 & 1 & 1.00 & 0 & 1 & 1 & 1 & --- \\
\rowcolor{gray!40}3 & 8 & 0.84 & 2.2e-77 & 8 & 5 & 8 & ST \\
\rowcolor{gray!40}5 & 21 & 0.71 & 4.7e-76 & 21 & 19 & 21 & KK \\
7 & 54 & 0.68 & 2.4e-75 & 48 & 52 & 55 & KK \\
9 & 120 & 0.66 & 2.0e-74 & 91 & 118 & 138 & KK \\
11 & 234 & 0.65 & 5.4e-73 & 160 & 233 & 272 & FR \\
13 & 416 & 0.65 & 5.9e-72 & 259 & 415 & 512 & FR \\
15 & 690 & 0.64 & 5.9e-71 & 400 & 687 & 728 & FR \\
17 & 1078 & 0.64 & 3.9e-71 & 589 & 1074 & 1384 & FR \\
19 & 1612 & 0.63 & 2.7e-70 & 840 & 1606 & --- & --- \\
21 & 2322 & 0.63 & 5.9e-71 & 1161 & 2315 & --- & --- \\
23 & 3244 & 0.63 & 8.0e-76 & 1568 & 3235 & --- & --- \\
\midrule
\multicolumn{8}{@{}l}{$d = 5$} \\
\rowcolor{gray!40}1 & 1 & 1.00 & 0 & 1 & 1 & 1 & --- \\
\rowcolor{gray!40}3 & 10 & 0.79 & 4.1e-77 & 10 & 6 & 10 & ST \\
\rowcolor{gray!15}5 & 32 & 0.67 & 1.2e-75 & 31 & 30 & 32 & ST \\
7 & 100 & 0.63 & 2.9e-73 & 80 & 100 & --- & --- \\
9 & 266 & 0.61 & 4.4e-72 & 171 & 264 & --- & --- \\
11 & 602 & 0.60 & 2.2e-72 & 332 & 598 & --- & --- \\
13 & 1212 & 0.59 & 4.6e-73 & 591 & 1204 & --- & --- \\
15 & 2500 & 0.60 & 4.6e-73 & 992 & 2224 & --- & --- \\
17 & 3872 & 0.58 & 1.2e-71 & 1581 & 3840 & --- & --- \\
19 & 6826 & 0.58 & 1.4e-78 & 2422 & 6278 & --- & --- \\
21 & 10984 & 0.58 & 1.1e-78 & 3583 & 9820 & --- & --- \\
\bottomrule
\end{tabular}
\normalsize
\caption{Positive-weight rules for the uniform weight on $[0,1]^d$. $N$: nodes in my best rule; $\rho = N^{1/d}/q$ with $q = (p+1)/2$ ($\rho = 1$ is the Gauss product grid, smaller is better); rel.\ err.: largest relative monomial error, $\max_{|a| \le p} |\sum_s w_s x_s^a - \mathrm{E}\,x^a| / \max(\sum_s w_s\,|x_s^a|, 1)$, of the rule's 80-digit file (rounded to double precision: at most $3.7\cdot 10^{-17}$); M\"o: M\"oller's lower bound; pair: the counting floor for a $\pm$pair rule, $2\lceil M_{\mathrm{even}}/(d+1) \rceil$ with $M_{\mathrm{even}}$ the number of monomials of even total degree $\le p$ (one less where a node at the origin allows it) --- a parameter count for that ansatz, not a proven bound; prev: best count in the literature, --- where none exists besides the product grid; src: where that count is recorded --- the original paper where I could identify it, otherwise a compilation (CH = \citealt{coolio88b}; DW = \citealt{diablo}; FR = \citealt{frontin21}; FS = \citealt{fiesta}; KK = \citealt{cash}; MO = \citealt{molly76}; ST = \citealt{stroud71}; XG = \citealt{xiao10}); *: my rule was warm-started from theirs. Light gray: no improvement on prev; dark gray: M\"oller's bound attained (proven minimal). Every such file also meets the stopping criterion of \citealt{diablo}, $\|f - V^{\top} w\|_2 < 10^{-66}$, over the full Legendre basis (appendix).}\label{tab:best-le}
\end{table}
A description of the results for the Gauss-Hermite and Gauss-Legendre cases can be found in \cref{tab:best-gh,tab:best-le}, respectively. Each table segment corresponds to a different dimension $d$. Each row in a table segment corresponds to a different polynomial exactness degree $p$.
There are four columns with node counts in each table: $N$ is the number of nodes used by the proposed rule, Mö is Möller's lower bound \citep{molly79}, and prev refers to the node count of the best available rule prior to this paper. The Möller bound is attainable for some combinations of $d,p$ but for most cases with $d≥3$ it is unknown and doubtful that rules achieving the Möller bound exist. So, the fact that a rule does not attain the Möller lower bound does not imply that a better rule exists. The fourth node count column in each table (pair in \cref{tab:best-le} and sym in \cref{tab:best-gh}) is a counting floor rather than a bound: \cref{app:methods} sets out how each floor is counted, and what it does and does not imply. Again, these floors are not known lower bounds (and in some cases are even less than the Möller bounds), but are a better reflection of what is achievable than the Möller bounds if $d$ or $p$ is large.
The constraints imposed are not the same across papers. For instance, \citet[\citetalias{diablo}]{diablo} impose a greater degree of symmetry on the solution than some of the methods I use do, so it is not surprising that the node counts that I provide are lower than those in \citetalias{diablo}, except where \citetalias{diablo} attain the Möller bound. As mentioned in \cref{app:methods}, the rules provided by \citetalias{diablo} were moreover used to warm-start my algorithms.
In addition, some papers contain some positive-weight rules even though weights were not restricted to be positive. This is e.g.\ true for \citet{sandy} who, in contrast to most papers cited here, puts emphasis on even $p$ cases. None of the \citet{sandy} rules provide the lowest node counts for the (odd) $p$ values presented here.
Another column in the tables contains the value of $ρ$, which reflects the equivalent one-dimensional count in the sense that the number $ρ=0.60$ in \cref{tab:best-gh} for $d=5$, $p=9$ ($q=5$) means that the node count in that case is equivalent to that of the tensor product for $q=3$ ($3^5 = 243 ≈ 244$).
\begin{figure}[t]
\centering
\begin{tikzpicture}
\begin{groupplot}[group style={group size=2 by 1, horizontal sep=1.2cm}, quad, width=0.5\linewidth, height=0.4\linewidth,
ymin=0.5, ymax=0.92, xmin=5, xlabel={degree $p$}, restrict x to domain=5:99]
\nextgroupplot[title={normal}, ylabel={$\rho$}, ylabel style={rotate=-90, font=\normalsize}, xmax=38.5]
\addplot[qd2, line width=0.9pt, mark=*, mark size=1pt] table[x=p, y=rho] {fig/rho_gh_d2.dat};
\addplot[qd2, line width=0.6pt, densely dashed] table[x=p, y=floor] {fig/rho_gh_d2.dat};
\pgfplotsinvokeforeach{3,4,5}{
\addplot[qd#1, line width=0.9pt, mark=*, mark size=1pt, forget plot] table[x=p, y=rho] {fig/rho_gh_d#1.dat};
\addplot[qd#1, line width=0.6pt, densely dashed, forget plot] table[x=p, y=floor] {fig/rho_gh_d#1.dat};
}
\node[anchor=west, font=\scriptsize] at (axis cs:33.4,0.852) {$d=2$};
\node[anchor=west, font=\scriptsize] at (axis cs:33.4,0.765) {$d=3$};
\node[anchor=west, font=\scriptsize] at (axis cs:23.4,0.632) {$d=4$};
\node[anchor=west, font=\scriptsize] at (axis cs:21.4,0.598) {$d=5$};
\nextgroupplot[title={uniform}, xmax=93]
\pgfplotsinvokeforeach{2,3,4,5}{
\addplot[qd#1, line width=0.9pt, mark=*, mark size=0.8pt] table[x=p, y=rho] {fig/rho_le_d#1.dat};
\addplot[qd#1, line width=0.6pt, densely dashed] table[x=p, y=floor] {fig/rho_le_d#1.dat};
}
\node[anchor=west, font=\scriptsize] at (axis cs:77.5,0.824) {$d=2$};
\node[anchor=west, font=\scriptsize] at (axis cs:45.5,0.712) {$d=3$};
\node[anchor=west, font=\scriptsize] at (axis cs:23.5,0.632) {$d=4$};
\node[anchor=west, font=\scriptsize] at (axis cs:21.5,0.578) {$d=5$};
\end{groupplot}
\end{tikzpicture}
\caption{The node count of each rule for $p≥5$ in \cref{tab:best-gh,tab:best-le} as a fraction of the Gauss product grid's, per axis: $\rho = N^{1/d}/q$ with $q = (p+1)/2$ (solid, one marker per rule), so that $\rho = 1$ is the product grid, and the counting floor of the tables' sixth column on the same scale (dashed).}
\label{fig:rho}
\end{figure}
As \cref{fig:rho} shows, in each table the $ρ$ values for a fixed value of $d$ flatten out as $p$ increases. The fact that the $ρ$ values tick up for larger values of $p$ in \cref{tab:best-gh} may suggest that those rules are suboptimal. The value of $p$ at which I stopped producing rules for a given value of $d$ is thus chosen as the value of $p≥21$ at which the value of $ρ$ started to edge up.
Each table contains two more columns: one displaying the relative error of the result and one the source where I found the previous best result.\footnote{Cools' Encyclopaedia of Cubature Formulas is the recorded source for sixteen of the Gauss-Hermite counts in \cref{tab:best-gh}. Its overview tables are online; I requested access to the full collection and received no reply, so those counts are as the Encyclopaedia records them rather than checked against the original publications.}
Dark cell coloring means that rules attaining the Möller lower bound already exist and light cell coloring that I have not found rules improving on what exists, albeit that my solution can differ from the preexisting results.
That said, most of the entries have white backgrounds, reflecting the case in which my results either improve on what exists or are the first to provide a solution different from the tensor product at all. Here, as throughout the paper, only positive-weight rules count: sparse-grid rules \citep{smolyak63}, for the normal density typically built on the nested one-dimensional rules of \citet{genz96}, exist for every $d$ and $p$, but they generally have some negative weights.
\clearpage