Perfectly equidistributed Quasi-Monte Carlo sequences from Artin-Schreier polynomials
July 16, 2026
1 Introduction↩︎
Numerically estimating high-dimensional integrals benefits from Quasi-Monte Carlo sampling. The resulting estimator averages evaluations of the integrand at specific pseudo-random locations \[\int_{[0, 1)^s} f(x) \text{d}x \approx \frac{1}{n} \sum_{i=1}^n f(x_i)\,,\] where \(\{x_i\}_{i=1..n}\) uniformly covers the integration domain (we will consider an \(s\)-dimensional unit hypercube, for simplicity). Uniformity can be measured by a non-negative integer value \(t\), where smaller \(t\) indicates better uniformity. In this context, a \((t,m,s)\)-net is a set of \(n=b^m\) points in \(s\) dimensions, where if we regularly subdivide the integration domain in rectangles of side lengths \(\frac{1}{b^{k_i}}\) for \(k_i \in \mathbb{N}, i \in \{1 \dots s\}\), then each rectangle of area \(\frac{b^t}{n}\) has exactly \(b^t\) points inside (see figure [fig:t1]). This \(t\) value amounts to a stratification property of the point set. High uniformity is desirable as it improves the integration error in \(\mathcal{O}\left(C(s,b) b^t \frac{\left(\log n\right)^s}{n}\right)\) for integrands of bounded Hardy and Krause variation [1] where \(C\) only depends on \(s\) and \(b\). A \((t,s)\)-sequence is a sequence of sample points where suitable subsets of \(b^m\) sample points all satisfy the \((t,m,s)\)-net property – this allows for progressive estimation of an integral, refining the estimated integral value by adding more points while maintaining low integration error.
The Sobol’ method is widely used as it produces \((t,s)\)-sequences. This construction uses a set of irreducible polynomials \(\{p_i\}_{i=1 \dots s}\) in the Galois Field \(GF(b)\) and small initialization matrices in \(GF(b)\) to produce larger matrices, that, in turn, are used to produce points. The crux of Sobol’ sequences lies in determining these polynomials and initialization matrices to ensure low \(t\) values. It can be shown that by taking \(p_i(x) = x + i\), for \(i=0 \dots s\) (\(s < b\)), the resulting \(s\)-dimensional sequence has \(t=0\) and the matrices produced relate to powers of Pascal matrices [2]. While cases for which \(t=1\) (when base \(b=2\)) have been found [3], we are not aware of other conditions guaranteeing \(t=0\) (for more general \(b\)) in the context of Sobol’ sequences.
In this paper, we show that similarly to the case of affine polynomials, we can guarantee \(t=0\) by considering a set of polynomials \(p_i(x) = p(x) + c_i\) of degree \(e\) that only differ by a constant \(c_i\). In this case, we show that the Sobol’ constructed matrices also amount to powers of Pascal matrices, tensorized by an \(e \times e\) matrix, and that \(t=0\) as soon as the first few Sobol’ iterations yield \(t=0\) point sets. To ensure irreducibility of these polynomials and compatibility with Sobol’, we then use Artin-Schreier irreducible polynomials of the form \(p_i(x) = x^b - x + c_i\) in \(\mathrm{GF}(b)\) with prime \(b\) and \(c_i \in \mathrm{GF}(b)^\times\) to obtain \(b-1\) irreducible polynomials in base \(b\). We further find that initializing these Sobol’ recurrences with powers of Pascal matrices left and right-multiplied by diagonal matrices always ensures that \(t=0\), and provide a greedy solution to explore \(t=0\) candidate initializations. We finally conclude with numerical experiments to evaluate the uniformity of our \((b-1)\)-dimensional Artin-Schreier sequences as well as \((2b-1)\)-dimensional sequences obtained by concatenating coordinates obtained by the \(b\) linear polynomials and the \(b-1\) Artin-Schreier polynomials.
2 Background↩︎
Digital nets
To produce a well-distributed set of \(b^m\) sample points in \(s\) dimensions, the idea of digital nets is to first construct \(s\) matrices \(\{M_i\}_{i=1\dots s}\) of size \(m \times m\) in the Galois Field1 \(\mathrm{GF}(b)\). To construct the \(j-\)th sample, one decomposes the index \(j\) in base \(b\) to form a column vector of digits (e.g., in base \(b=2\), this is the binary representation of \(j\)). Then each matrix multiplies this column vector to obtain \(s\) new column vectors, that are reinterpreted as \(s\) coordinates in \([0, 1[\) by considering these vectors as base-\(b\) digits with weights \(b^{-1}, b^{-2}, \dots\). The uniformity of samples produced with this approach is dictated by the property of matrices \(\{M_i\}_{i=1\dots s}\). For digital nets, the quality factor \(t\) can be obtained directly from matrices [1], [4]. In particular, a sequence with \(t=0\) amounts to ensuring that if one builds a matrix by stacking the \(d_i\) first rows of matrix \(M_i\), with \(\sum_{i=1..s} d_i = m\), the \(m\times m\) compound matrix \(\mathcal{K}\) has \(\text{det}(\mathcal{K}) \neq 0\) (in \(\mathrm{GF}(b)\)) for all \(m\) and all \(\{d_i\}\).
Sobol’ construction
A well-established way of producing good matrices is through Sobol’ construction [5]. Given a primitive polynomial \(p\) of degree \(e\) in base \(b\), the original Sobol’ approach recursively builds an upper triangular matrix one column at a time, as a combination of the \(e\) previous columns. One matrix is produced for each dimension, and thus \(s\) dimensions require \(s\) polynomials and \(e \times e\) matrices to initialize the recursion. Numerical search over good \(e \times e\) matrices was notably performed by Joe and Kuo [6] for \(b=2\), and then Ostromoukhov et al. [7] for \(b=3\).
The core of our paper uses the construction of Faure and Lemieux [8] who obtained the same Sobol’ matrices by operating the recursive construction blockwise rather than columnwise. For a Sobol’ matrix \(M\) and polynomial \(p(x) = \sum_{i=0}^e a_i x^i\) of degree \(e\), their recurrence reads: \[\label{eq:faurelemieux} M_{i,j} = \left(M_{i,j-1} Q_p + M_{i-1,j-1}\right) F_p\tag{1}\] where \(i,j\) index the \(i^{\text{th}}\) block row and \(j^{\text{th}}\) block column of matrix \(M\), with blocks of size \(e \times e\), and \(Q_p\) and (inverse of) \(F_p\) are matrices2 given by [3], [8]:
\[\label{eq:q} Q_p = \left[ \begin{array}{ccccc} -a_0 & 0 & 0 & \dots & 0 \\ -a_1 & -a_0 & 0 & \dots & 0 \\ -a_2 & -a_1 & -a_0 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -a_{e-1} & -a_{e-2} & -a_{e-3} & \dots & -a_0 \\ \end{array} \right]\tag{2}\]
\[\label{eq:f} F^{-1}_p = \left[ \begin{array}{ccccc} 1 & -a_{e-1} & -a_{e-2} & \dots & -a_{1} \\ 0 & 1 & -a_{e-1} & \dots & -a_{2} \\ 0 & 0 & 1 & \dots & -a_{3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & 1 \\ \end{array} \right]\,.\tag{3}\]
Faure and Lemieux show that using the larger class of irreducible polynomials instead of primitive polynomials, with \(b\) a prime power, also yields \((t,s)\)-sequences.
3 Overview↩︎
Our paper investigates conditions that guarantee the highest quality \(t=0\) in the context of Sobol’ sequences in \(\mathrm{GF}(b)\). In base \(b\), one may find \(s\) irreducible degree-\(e\) polynomials that only differ by a constant (we call them consecutive polynomials), although they may not always exist. For instance, in base \(b=5\), we may find consecutive degree \(e=5\) polynomials \(x^5+4x+1\), \(x^5+4x+2\), \(x^5+4x+3\) and \(x^5+4x+4\) that are all irreducible. In section 4, we show that given 2 consecutive polynomials, the Sobol’ matrices produced are equivalent to powers of Pascal matrices tensorized by \(e \times e\) matrices. This relates our construction to Faure’s construction of \((0,s)\)-sequences that corresponds to the special case where \(e=1\) [2]. However, while in the case \(e=1\) the \((0,s)\)-sequence condition is automatically satisfied, for \(e>1\) this is not always the case. In section 5, we show that testing the \((0, m, s)\)-net conditions for the first few submatrices guarantees the \((0,s)\)-sequence property up to infinity. In section 6, we use Artin-Schreier theory to obtain \(s=b-1\) degree-\(b\) consecutive irreducible polynomials in \(\mathrm{GF}(b)\), and show that initializing the Sobol’ construction with powers of Pascal matrices always fulfills the conditions for \(t=0\). In section 7 we numerically experiment our sequences with a simple greedy algorithm to explore good initialization matrices, and evaluate discrepancies of the resulting point sets.
Settings
We consider a prime base \(b\) and \(p(x)\in\mathrm{GF}(b)[x]\) a polynomial of degree \(e\ge1\). We assume that for \(\{c_i\in\mathrm{GF}(b)\}_{i=1\dots s}\), the shifted polynomials \(p_i = p+c_i\) are distinct irreducible polynomials of degree \(e\). We denote \(M^{(c_i)}\) their corresponding Sobol’ matrix.
images/matrices/linearPolynomial_GF5_dim4.pdf (600,-100) (630,-92)10 (630,-74)9 (630,-58)8 (630,-42)7 (630,-24)6 (630,-8)5 (630,8)4 (630,26)3 (630,44)2 (630,62)1 (630,80)0 |
||||||||||||
Figure 2: Color-coded representation of Sobol’ matrices produced by consecutive polynomials \(p_i(x) = x + c_i\) (top row) and Artin-Schreier polynomials \(p_i(x) = x^b - x + c_i\) (bottom row) in \(\mathrm{GF}(5)\), \(\mathrm{GF}(7)\) and \(\mathrm{GF}(11)\).. a — image, b — image, c — image, d — image, e — image, f — image, g — image, h — image, i — image, j — image, k — image, l — image, m — image, n — image, o — image, p — image, q — image, r — image, s — image, t — image, u — image, v — image, w — image, x — image, y — image, z — image, { — image, | — image, } — image, ~ — image, — image, — image, — image, — image, — image, — image, — image, — image, — image, — image, — image, — image, — image
4 Consecutive polynomials produce Pascal powers↩︎
We first focus on the specific case of two-dimensional sequences, produced by Sobol’ recurrence with consecutive polynomials, producing matrices \(M^{(c_1)}\) and \(M^{(c_2)}\).
Characteristic matrices
For two-dimensional sequences, Ahmed et al. [9] define the characteristic matrix as \(K = M^{(c_2)}\left(M^{(c_1)}\right)^{-1}\). They show that this matrix uniquely characterizes a digital net: the 2D point set generated by \((M^{(c_1)},M^{(c_2)})\) is the same as the one given by \((\mathrm{Id},K)\), up to reordering.
We show that consecutive polynomials yield characteristic matrices in the form of tensorized powers of Pascal matrices:
Lemma 1. Given two consecutive polynomials over \(\mathrm{GF}(b)\) of degree \(e\) and their corresponding Sobol’ matrices \(M^{(c_1)}\) and \(M^{(c_2)}\), the characteristic matrix \(K = M^{(c_2)} \left(M^{(c_1)}\right)^{-1}\) restricted to a size \(em \times em\) (denoted \(K_{[em \times em]}\)) for \(m\geq 1\), is a Pascal matrix at the power \(\alpha=c_1-c_2\) tensorized with an initialization matrix \(K_{[e \times e]}\): \[K_{[em \times em]} = \mathcal{P}^{\alpha} \otimes K_{[e \times e]}\,,\] where \(\mathcal{P}^{\alpha}\) is an \(m \times m\) matrix of coefficients in \(\mathrm{GF}(b)\): \[\mathcal{P}^{\alpha}_{i,j} = \binom{j-1}{i-1}\alpha^{j-i}\,.\]
Equivalently, the block coefficients of the (infinite) matrix \(K\) are given by: \[\label{eq:rec} K_{i,j} = \alpha\,K_{i,j-1}+K_{i-1,j-1}, \qquad K_{i,1}=\delta_{i,1}\,K_{[e \times e]}.\tag{4}\]
Proof. We show that matrix \(L := \mathcal{P}^{\alpha} \otimes K_{[e \times e]}\) satisfies \(L M^{(c_1)}_{[em \times em]} = M^{(c_2)}_{[em \times em]}\). We exploit an identity relating powers of Pascal matrices [2]: \(\mathcal{P}^a_{i,j} = a^{j-i} \binom{j-1}{i-1}\), and rewrite it for \(e \times e\) block \((i,j)\) of matrix \(L\): \(L_{i,j+1} = \alpha L_{i,j}+L_{i-1,j}\). The matrices of equations 2 and 3 are such that \(F_{p+c_i} = F_p\) and \(Q_{p+c_i}=Q_p-c_i \,\text{Id}_{e\times e}\). Using these identities, \(\alpha = c_1-c_2\) and equation 1 , and denoting \(N = L M^{(c_1)}_{[em \times em]}\), we develop the block matrix-vector multiplication as \[\begin{align} N_{i,j} & =\sum_{k} L_{i,k}\,M^{(c_1)}_{k,j-1}\,(Q_p-c_1 \text{Id}_{e\times e})F_p + L_{i,k}\,M^{(c_1)}_{k-1,j-1}\,F_p\\ & = N_{i,j-1}(Q_p-c_1 \text{Id}_{e\times e})F_p + \sum_{k} L_{i,k+1} M^{(c_1)}_{k,j-1} F_p\\ &= N_{i,j-1}(Q_p-c_1 \text{Id}_{e\times e})F_p +\bigl(\alpha\,N_{i,j-1}+N_{i-1,j-1}\bigr)F_p\\ &=\bigl(N_{i,j-1}(Q_p-c_2 \text{Id}_{e\times e})+N_{i-1,j-1}\bigr)F_p, \end{align}\] Thus \(N\) satisfies the Faure and Lemieux recurrence of equation 1 for polynomial \(p+c_2\). The initialization is trivial, with \(N_{i,1}=L_{i,1} M^{(c_1)}_{[e \times e]} = \delta_{i,1} K_{[e \times e]} M^{(c_1)}_{[e \times e]} = \delta_{i,1} M^{(c_2)}_{[e \times e]}\), and thus \(N=M_{[em\times em]}^{c_2}\). \(L\) is thus the characteristic matrix of \(M^{(c_1)}\) and \(M^{(c_2)}\); in other words \(L = K_{[em \times em]}\) which concludes the proof. ◻
5 Pascal powers may give \((0,s)\)-sequences↩︎
We now move to the more general \(s\)-dimensional setting. The two-dimensional observation of Ahmed et al. [9] has been shown in the higher dimensional settings by Faure and Tezuka [10]. Given matrices \(M^{(c_1)}, M^{(c_2)}, \dots, M^{(c_{s})}\) producing a \((t,s)\)-sequence, right-multiplying these matrices by the same invertible matrix results in the same set of points, up to reordering, and thus they also produce a \((t,s)\)-sequence. One may choose right-multiplying by \(\left(M^{(c_1)}\right)^{-1}\). This brings the study of Sobol’ sequences produced by consecutive polynomials to the simpler study of tensorized powers of Pascal matrices.
We show that consecutive polynomials may lead to \((0,s)\)-sequences.
Theorem 1. Given \(s\) distinct consecutive polynomials of degree \(e\) over \(\mathrm{GF}(b)\), if the \(s-\)dimensional Sobol’ point set of \(b^m\) samples forms a \((0,m,s)\)-net for all \(m\le s\,(e-1)\) then it forms a \((0, s)-\)sequence in base \(b\).
Proof. For \(1\leq i\leq s\), with \(\alpha_i = c_1-c_{i}\), and \(K_{[m \times m]}^{(\alpha_i)}:=M_{[m \times m]}^{(c_i)}\left(M_{[m \times m]}^{(c_1)}\right)^{-1}\), the matrices \(\left\{K^{(\alpha_i)}_{[m \times m]} \right\}_{i=1\dots s}\) produce the same sequence as matrices \(\{M_{[m \times m]}^{(c_i)}\}_{i=1\dots s}\). In the following, we drop the \([m \times m]\) subscript for clarity. As noted in Sec. 2, the first \(b^m\) points form a \((0,m,s)\)-net if and only if every possible \(m \times m\) compound matrix \(\mathcal{K}_{\boldsymbol{d}}\), for all \({\boldsymbol{d}} = \{d_i\}_{i=1\dots s}\), formed by stacking the first \(d_i\) rows of \(K^{(\alpha_i)}\) has \(\text{det}(\mathcal{K}_{{\boldsymbol{d}}}) \neq 0\) (see figure 3). We assume that this is the case up to \(m \leq s (e-1)\) and argue by contradiction. Suppose that some compound matrix of total size \(m > s (e-1)\) is singular. Thus, for some row counts \({\boldsymbol{d}} = d_1,\dots,d_s\) with \(\sum_i d_i=m\), there exists a nonzero vector \(\mathbf{v}\in\mathrm{GF}(b)^m\) annihilated by the selected rows.
We group the entries of \(\mathbf{v}\) into blocks of length \(e\), \[\mathbf{v}=(\mathbf{v}_0,\mathbf{v}_1,\ldots),\qquad \mathbf{v}_\ell\in\mathrm{GF}(b)^e,\] padding the last block with zeros if necessary, and form the vector polynomial \[F(x)=\sum_{\ell\ge 0} \mathbf{v}_\ell x^\ell.\] Thus singularity of \(\mathcal{K}_{\boldsymbol{d}}\) is equivalent to the existence of a nonzero coefficient vector, or equivalently a nonzero vector polynomial \(F\), such that \(\mathcal{K}_{\boldsymbol{d}}\mathbf{v}=0\).
We note that powers of Pascal matrices act as polynomial shifts, i.e., for a scalar-valued polynomial \(p(x) = \sum_i a_i x^i\), with \(\mathbf{a}=(a_i)\), denoting \(\mathbf{a}' = \mathcal{P}^y \mathbf{a}\), then \(\sum_{i} a'_i x^i = p(x+y)\).
For our vector-valued polynomial \(F\) and tensorized Pascal matrices, the \(q\)-th block row of one matrix \(K^{(\alpha_i)} = \mathcal{P}^{\alpha_i}\otimes K^{(\alpha_i)}_{[e \times e]}\) hence computes \[K^{(\alpha_{i})}_{[e \times e]} \, \text{coeff}_{x^q} F(x+\alpha_i),\] where \(\text{coeff}_{x^q} F\) denotes the coefficient of \(x^q\) in \(F\).
Turning to the compound matrix \(\mathcal{K}_{\boldsymbol{d}}\), we write each row count as \(d_i=e q_i+r_i\), for \(0\le r_i<e\). The first \(e q_i\) rows of \(K^{(\alpha_i)}\) consist of \(q_i\) complete block rows. Since \(K^{(\alpha_i)}_{e \times e}\) is invertible, these rows force \[\text{coeff}_{x^q} F(x+\alpha_i)=0, \qquad q=0,\ldots,q_i-1.\] Equivalently, \(F(x)\) is divisible by \((x-\alpha_i)^{q_i}\). Since the \(\alpha_i\) are distinct, we can write for some vector polynomial \(H(x)\neq 0\): \[F(x)=\Phi(x)H(x), \qquad \Phi(x)=\prod_{i=1}^s (x-\alpha_i)^{q_i}.\]
Let \(Q=\sum_i q_i = \deg \Phi\), and \(R=\sum_i r_i\). Since \(\Phi\) is monic of degree \(Q\) and \(F\) has at most \(m=eQ+R\) scalar coefficients, \(H\) has at most \(R\) scalar coefficients. If \(R=0\), this implies \(H=0\), contradicting \(H\neq 0\). We may thus assume \(R>0\) and handle the remaining incomplete blocks of rows.
For a fixed \(i\), the remaining \(r_i\) selected rows of \(K^{(\alpha_i)}\) are precisely the first \(r_i\) rows inside the next block row \(q_i\). They therefore impose that the first \(r_i\) components of \(K^{(\alpha_i)}_{[e \times e]} \,\text{coeff}_{x^{q_i}}F(x+\alpha_i)\) vanish.
We factor \(\Phi\) by \((x-\alpha_i)^{q_i}\): \[\Phi(x)=(x-\alpha_i)^{q_i}\Phi_i(x), \qquad \Phi_i(x)=\prod_{j\ne i}(x-\alpha_j)^{q_j}\,,\] then \[F(x+\alpha_i) = x^{q_i}\Phi_i(x+\alpha_i)H(x+\alpha_i)\,,\] and therefore \[\text{coeff}_{x^{q_i}}F(x+\alpha_i) = \Phi_i(\alpha_i)H(\alpha_i)\,.\]
Since \(\Phi_i(\alpha_i)=\prod_{j\ne i}(\alpha_i-\alpha_j)^{q_j}\ne 0\), this is equivalent to the vanishing of the first \(r_i\) components of \(K^{(\alpha_i)}_{[e \times e]} H(\alpha_i)\) for every \(i\). In other words, the nonzero coefficient vector of \(H\), of length \(R=\sum_i r_i\le s(e-1)\), is annihilated by the first \(r_i\) rows of each \(K^{(\alpha_i)}\), i.e., by the compound matrix \(\mathcal{K}_{(r_1,\dots,r_s)}\), which the hypothesis forbids. Hence every compound matrix is non-singular and the sequence is a \((0,s)\)-sequence. ◻
6 Guaranteeing \(t=0\)↩︎
The conditions of Sec. 5 would suggest searching for \(s\) initialization matrices, each of size \(e \times e\), and identifying cases that satisfy conditions of theorem 1. While this works in practice for small values of bases \(b\), \(e\) and \(s\), it becomes intractable even for moderate values. Even focusing on invertible upper triangular initialization matrices, this would lead to \(b^{e(e-1) s/2} (b-1)^{e\, s}\) degrees of freedom, such that even 3-dimensional sequences in base \(b=5\) using degree \(e=5\) polynomials would require a search over \(10^{30}\) candidates.
When \(e=b\), we formulate the following condition to guarantee that sequences have \(t=0\).
Proposition 2. Considering \(s\) distinct consecutive polynomials of degree \(e=b\) over \(\mathrm{GF}(b)\), initializing their \(s\) Sobol’ matrices with \[M^{(c_i)}_{[e \times e]} = D \mathcal{P}^{c_i}_{[e\times e]} D^{-1}\] with \(D\) a diagonal invertible matrix, and \(\mathcal{P}_{[e\times e]}\) the \(e\times e\) Pascal matrix always produces \((0,s)\)-sequences.
To show proposition 2, we first exhibit the self-similar structure of powers of Pascal matrices.
Lemma 2 (Self-similarity of Pascal powers). Let \(b\) be prime and \(a\in\mathrm{GF}(b)\). Over \(\mathrm{GF}(b)\), \[\label{eq:tensor-fine} \mathcal{P}^{a}_{[bm\times bm]} =\mathcal{P}^{a}_{[m\times m]}\otimes\mathcal{P}^{a}_{[b\times b]} \qquad\text{for every }m\ge1,\tag{5}\] and, when \(m\) is a power of \(b\), \[\label{eq:tensor-coarse} \mathcal{P}^{a}_{[bm\times bm]} =\mathcal{P}^{a}_{[b\times b]}\otimes\mathcal{P}^{a}_{[m\times m]}.\tag{6}\]
Proof. We exploit the identity \(\mathcal{P}^a_{i,j} = a^{j-i} \binom{j-1}{i-1}\). For equation 5 , split off the last base-\(b\) digit of the (0-based) indices: \(i=\hat{\imath}\,b+i_0\), \(j=\hat{\jmath}\,b+j_0\) with \(0\le i_0,j_0<b\). Lucas’ congruence [11], [12] gives \(\binom{j}{i}\equiv\binom{\hat{\jmath}}{\hat{\imath}}\binom{j_0}{i_0}\pmod b\), and Fermat’s little theorem gives \(a^{\,j-i}=\bigl(a^{\,b}\bigr)^{\hat{\jmath}-\hat{\imath}}a^{\,j_0-i_0} =a^{\,\hat{\jmath}-\hat{\imath}}\,a^{\,j_0-i_0}\); the product of the two congruences is the \((i,j)\) entry of \(\mathcal{P}^{a}_{[m\times m]}\otimes\mathcal{P}^{a}_{[b\times b]}\). For 6 , split off the leading digit instead, \(i=I\,m+i'\), \(j=J\,m+j'\) with \(0\le I,J<b\): this is a single base-\(b\) digit exactly when \(m\) is a power of \(b\), and Fermat iterates to \(a^{\,m}=a\); the same two congruences conclude. ◻
We now turn to the proof of proposition 2.
Proof. As noted by Faure and Lemieux [8], a Sobol’ matrix obtained using an initialization \(M_{[e\times e]}\) can also be obtained by initializing the recurrence with an identity matrix and left-multiplying the result by a block-diagonal matrix with \(M_{[e\times e]}\) as the diagonal blocks. Denoting \(W_i = M^{(c_i)}_{[e \times e]} = D\,\mathcal{P}^{c_i}_{[e\times e]}D^{-1}\) the initialization matrix: \(M^{(c_i)}=\mathrm{diag}(W_i,W_i,\dots)\,M^{(c_i)}_{\mathrm{Id}}\). Using lemma 1 for the identity-initialized matrices, \[\begin{align} K^{(\alpha_i)} &=\mathrm{diag}(W_i,\dots)\,\bigl(\mathcal{P}^{\alpha_i}\otimes \mathrm{Id}_e\bigr)\, \mathrm{diag}(W_1,\dots)^{-1} \\ &=\mathcal{P}^{\alpha_i}\otimes W_iW_1^{-1}. \end{align}\] Pascal powers compose as the coordinate shifts they encode, \(\mathcal{P}^{a}\mathcal{P}^{a'}=\mathcal{P}^{a+a'}\), and conjugating by a geometric diagonal rescales the power, \(\mathrm{diag}(\lambda^u)_u\,\mathcal{P}^{a}\,\mathrm{diag}(\lambda^{-u})_u =\mathcal{P}^{a/\lambda}\) (since \(\lambda^{u}\binom{v}{u}a^{v-u}\lambda^{-v} =\binom{v}{u}(a/\lambda)^{v-u}\)); hence, with \(S=\mathrm{diag}\bigl((-1)^u\bigr)_u\) and \(D'=DS\), \[W_iW_1^{-1}=D\,\mathcal{P}^{\,c_i-c_1}_{[e\times e]}\,D^{-1} =D'\,\mathcal{P}^{\,\alpha_i}_{[e\times e]}\,D'^{-1}, ~~~ \alpha_i=c_1-c_i .\] Because \(e=b\), lemma 2 gives \(\mathcal{P}^{\alpha}\otimes\mathcal{P}^{\alpha}_{[b\times b]}=\mathcal{P}^{\alpha}\), and therefore \[K^{(\alpha_i)}=\Delta\;\mathcal{P}^{\alpha_i}\;\Delta^{-1}, \qquad \Delta=\mathrm{diag}(D',D',\dots),\] with the same diagonal \(\Delta\) for every \(i\). Diagonal conjugation preserves the nonsingularity of compound matrices: the left factor rescales each selected row by a nonzero scalar, and the right factor is a common non-singular triangular factor. The compound matrices of the family \(\{\mathcal{P}^{\alpha_i}\}\), with pairwise-distinct powers \(\alpha_i\) are all non-singular by Faure’s theorem [2]. All compound determinants are therefore nonzero, and the sequence is a \((0,s)\)-sequence. ◻
Artin-Schreier polynomials
While our proofs have not relied on polynomial irreducibility, enforcing the use of irreducible polynomials allows to maintain compatibility with Sobol’ sequences. By using \(s_1\) irreducible consecutive polynomials and \(s_2\) other polynomials commonly used for Sobol’ sequences [6], one may obtain samples of \(s_1+s_2\) dimensions with \(t=0\) guarantees in the \(s_1\) subspace, and preserving a bounded \(t\) value for the full \((s_1+s_2)\)-dimensional sequence. The question arises of whether \(s\) consecutive irreducible Sobol’ sequences always exist. The answer lies within Artin-Schreier theory [13]. This theory states that for prime \(b\), the \(b-1\) polynomials in \(\mathrm{GF}(b)\) of the form \(p_i(x) = x^b - x +c_i\) are always irreducible. Such polynomials perfectly fit our framework since proposition 2 requires \(e=b\), which also results from Artin-Schreier theory, and since \(s=b-1\) is the maximum number of possible different consecutive irreducible polynomials in \(\mathrm{GF}(b)\), apart from the \(b\) linear polynomials \(p_i(x) = x +c_i\). In practice, we use Artin-Schreier polynomials to define our sequences and for our numerical experiments. Illustrations of Sobol’ matrices produced by \(\{x+c_i\}\) and Artin-Schreier polynomials are shown in figure 2 for various bases.
7 Numerical experiments↩︎
7.1 Artin-Schreier sequences↩︎
In figure 4 (left), we evaluate a uniformity measure called generalized \(\ell^2\) discrepancy, which relates Monte Carlo integration error to point set uniformity [14] (lower is better) for our sequences with all \(4^4 = 256\) possible initialization \(D\) matrices in
\(\mathrm{GF}(5)\), using the tms_lib library [15]. For a chosen \(D\), we show
in figure 5 the pairwise 2D projections of point sets with \(b^m \approx 1200\) samples (for the nearest perfect power of \(b\)) corresponding to Sobol’ sequences
with Artin-Schreier polynomials in bases \(b=5, 7\) and \(11\) for up to 5 dimensions (the case \(b=3\) corresponds dimensions 13–14 of the \(\mathrm{GF}(3)\) sequence of Ostromoukhov et al. [7]). We also show other Sobol’-based sequences, such as the improvement of Faure
and Lemieux [8] using irreducible polynomials in \(\mathrm{GF}(2)\), the sequence of Bonneel et
al. [3] that uses polynomials of the form \(p\) and \(p^2+p+1\) to guarantee \(t=1\) in consecutive pairs of dimensions in \(\mathrm{GF}(2)\), and the sequence of Ostromoukhov et al. [7] that numerically optimizes quadruples of Sobol’ dimensions in \(\mathrm{GF}(3)\). In figure 6 (top), we compare the generalized \(\ell^2\) discrepancy of our sequences to these alternatives, as well as Joe and Kuo’s Sobol’ sequence [6], where discrepancy is computed after randomization by a base-\(b\) Owen scrambling [16] and averaged over 32 realizations.
7.2 Combined Sobol’-Artin-Schreier sequences↩︎
In this section we investigate the behavior of sequences obtained by combining the first \(b\) irreducible affine polynomials that produce powers of Pascal Sobol’ matrices (with no tensorization involved) and the \(b-1\) Artin-Schreier polynomials to form \((2b-1)\)-dimensional sequences. In practice, proposition 2 provides a way to construct \((b-1)^{b-1}\) different sequences (since multiplying \(D\) by a constant does not change the initialization matrix), all guaranteeing \(t=0\) in \(b-1\) dimensions, by exploring potential diagonal matrices \(D\). When combining these \(b-1\) dimensions with other Sobol’ dimensions, it can be desirable to optimize criteria over the combined \(2b-1\) dimensions. For instance, while taking \(D = \text{Id}_{b \times b}\) still produces a \(t=0\) sequence in \(b-1\) dimensions, the sequence obtained by combining the powers of Pascal matrices resulting from affine polynomials with the powers of Pascal matrices obtained by the initialization of proposition 2 with an identity diagonal matrix would result in the worst possible behavior among possible diagonal matrices in \(2b-1\) dimensions, and especially for sample counts \(n < b^b\) (see figure 4 for the combined \((2b-1)\)-dimensional discrepancy for all initializations in \(\mathrm{GF}(5)\)). While we rely on exhaustive search over these \((b-1)^{b-1}\) values for up to \(\mathrm{GF}(7)\), this can be intractable for larger \(b\). We alleviate this issue using a greedy approach. We first define our criterion as a lexicographic order of \(t\) values computed for each \(m\) as computed by the fast approach of Marion et al. [4], and progressively test diagonal elements \(D_{i,i}\) for increasing \(i\) of matrix \(D\), making the search only over \((b-1)^2\) values. The resulting \((2b-1)\)-dimensional discrepancy of our optimized \(D\) values for \(\mathrm{GF}(5)\), \(\mathrm{GF}(7)\) and \(\mathrm{GF}(11)\) can be seen in figure 6, compared to other Sobol’ samplers. Our resulting diagonal matrices contain the following values in the diagonal: \([1, 2, 3, 1, 4]\) in \(\mathrm{GF}(5)\), \([1, 2, 6, 5, 1, 2, 6]\) in \(\mathrm{GF}(7)\) and \([1, 2, 1, 1, 8, 7, 6, 9, 5, 1, 2]\) in \(\mathrm{GF}(11)\).
8 Conclusions↩︎
We showed that Sobol’ sequences in \(\mathrm{GF}(b)\) constructed with polynomials of degree \(e=b\) that only differ by a constant relate to Faure sequences [2] and correspond to tensorized powers of Pascal matrices. These can be initialized in such a way that they guarantee the highest quality \(t=0\). While irreducibility is not required per se when only using consecutive polynomials, one can ensure compatibility with Sobol’ sequences by using irreducible polynomials to guarantee that higher-dimensional sequences formed by concatenating Sobol’ sequences obtained by consecutive polynomials and other irreducible polynomials still produce \((t,s)\)-sequences. In this setting, choosing Artin-Schreier polynomials ensure that \(b-1\) irreducible polynomials can be found. We show that the resulting sequences are competitive in terms of discrepancy.
9 Materials and Methods↩︎
Disclosure of Delegation to Generative AI↩︎
The authors declare the use of generative AI in the research and writing process. According to the GAIDeT taxonomy (2025), the following tasks were delegated to GAI tools under full human supervision:
Research design
Proofreading and editing
The GAI tool used was: ChatGPT 5.5, Claude Opus 4.8. Responsibility for the final manuscript lies entirely with the authors. GAI tools are not listed as authors and do not bear responsibility for the final outcomes. Declaration submitted by: Nicolas Bonneel
Data, Materials and Software Availability↩︎
Our code for optimizing over initializations and comparing discrepancy graphs to other sequences is available at https://github.com/liris-origami/Artin-Schreier-sequences