1 Introduction↩︎

Almost Perfect Nonlinear (APN) functions are vectorial Boolean functions \(F\colon\mathbb{F}_2^n\to\mathbb{F}_2^n\) achieving the minimum possible differential uniformity \(\delta_F = 2\), i.e., for every nonzero \(a\) and every \(b\) the equation \(F(x\oplus a)\oplus F(x) = b\) has at most two solutions. They were introduced to cryptography by Nyberg [5] as optimal S-box designs for resistance against differential cryptanalysis [6]: any bijective function with \(\delta_F = 2\) achieves the smallest possible differential probability in a single round. APN permutations (bijective APN functions) are particularly sought, but their existence is settled only for odd \(n\); for even \(n \geq 8\) the question of whether an APN permutation exists is a major open problem [7].

The classification of APN functions up to CCZ-equivalence [8] is an active research program. CCZ-equivalence is the broadest equivalence preserving the APN property; two S-boxes in the same CCZ-class are interchangeable for cryptographic purposes. For \(n \leq 6\) the classification is essentially complete; for \(n = 7\) and \(n = 8\) only partial results exist.

1.1 Known APN functions for \(n=8\)↩︎

The classical infinite families of APN functions are the Gold power functions \(x^{2^i+1}\) [9] and the Kasami functions \(x^{4^i-2^i+1}\) [10]. For \(n = 8\), the Gold functions with \(\gcd(i,8)=1\) are quadratic (degree 2) and form exactly two CCZ-classes: \(\{x^3, x^{129}\}\) and \(\{x^9, x^{33}\}\) (the Frobenius orbits of \(i=1,7\) and \(i=3,5\) respectively). All other APN power functions for \(n=8\) (Kasami \(x^{21}\), \(x^{85}\), etc.) have degree \(> 2\) and are therefore not quadratic.

Beyond power functions, a series of quadratic APN constructions appeared over two decades: Edel and Pott [11] found 8 157 new instances by a bivariate construction; Yu, Wang, and Li [12] obtained more via a matrix method; Weng, Tan, and Gong [13] contributed additional instances via algebraic characterization; Taniguchi [14] described a parametric family; Budaghyan et al. [15] introduced isotopic shifts. Beierle, Brinkmann, and Leander [1] systematically studied self-equivalence subspaces (see Section 3) and found 12 921 new instances; a companion paper [16] found a further 11 779 instances via extension methods. Together with the above, approximately 32 892 CCZ-inequivalent quadratic APN functions for \(n=8\) were known through 2021; Beierle, Leander, and Perrin [17] added 6 368 more via trimming and extension.

In 2025, Beierle, Langevin, Leander, Polujan, and Rasoolzadeh [3] extended this dramatically: using two construction strategies based on vectorial bent functions and an extensive computer search, they produced 3 775 599 new functions, estimating the total number of CCZ-classes at approximately 6 million.

1.2 Motivation and approach↩︎

Despite the large size of the Beierle 2025 database, certain structured algebraic regions of the search space remain unexplored. A quadratic function \(F\colon\mathbb{F}_{2^{8}}\to\mathbb{F}_{2^{8}}\) is described by a binary coefficient vector of length 224 (see Section 2). The full space has size \(2^{224}\); any structured subspace of dimension \(d \ll 224\) has probability \(2^{d-224}\) of being hit by a random function.

The self-equivalence subspace \(V_{\!A}\) defined by a linear automorphism \(A\) (see Section 3) is one such structured subspace. For the automorphism \(A\) studied here (\(\mathrm{ord}(A) = 5\)), this subspace has dimension 40. The probability of a random quadratic function lying in \(V_{\!A}\) is \(2^{40-224} = 2^{-184} \approx 10^{-55}\), so global random search never reaches it. Within \(V_{\!A}\), however, the APN density is approximately \(10^{-4}\) (see Section 4.4), making targeted search tractable.

The recursive tree search of BBL2021 found zero APN functions in this particular subspace. We show that a Gröbner basis approach combined with explicit RREF parameterization of \(V_{\!A}\) succeeds where the tree search did not.

1.3 Contributions↩︎

This paper makes the following contributions:

1. We establish that the self-equivalence subspace associated with BBL2021 class index 22 (the automorphism of order 5 described above) does contain quadratic APN functions, in contrast to the zero-solution result of the original BBL2021 search.

2. We develop an explicit RREF parameterization (fc22\(+\)sol22) of this 40-dimensional subspace, enabling efficient random sampling with APN density \(\approx 10^{-4}\).

3. We implement an NL = 4 Gröbner basis slice computation in Magma that, given any APN center in the subspace, finds all APN functions in a 24-dimensional hyperplane through that center.

4. Over two batches (428 hyperplanes, \(0.65\%\) of the total \(65\,536\)) we find 566 quadratic APN functions in six CCZ-classes. Four of these classes are absent from all known databases.

5. We provide a complete verification pipeline using the ortho-derivative invariant of [2] and confirm via streaming comparison that the four new classes do not appear in the Beierle 2025 database or any prior compilation.

The approach generalizes: the same pipeline applies to every BBL2021 class index with manageable subspace dimension, offering a systematic program for further extending the classification.

1.4 Organization↩︎

Section 2 fixes notation and recalls required background. Section 3 reviews the BBL2021 framework. Section 4 describes the specific subspace studied. Section 5 presents the search pipeline in detail. Section 6 reports experimental results. Section [sec:verify] describes the verification procedure and states the main theorem. Section 7 discusses the results and their limitations. Section 8 surveys related work. Section 9 concludes. Appendix 10 gives the S-box tables of the four new classes. Appendix 11 describes the computational infrastructure.

2 Preliminaries↩︎

Throughout, \(n = 8\) and \(N = 2^n = 256\), unless otherwise stated. We identify \(\mathbb{F}_{2^{8}}\) with \(\mathbb{F}_2^8\) and write \(\oplus\) for bitwise XOR. Scalar and vector indices run from 0.

2.1 Vectorial Boolean functions and differential uniformity↩︎

A vectorial Boolean function (VBF) is a map \(F\colon\mathbb{F}_2^n\to\mathbb{F}_2^n\). It may be written coordinate-wise as \(F = (f_0, \ldots, f_{n-1})\) with each \(f_j\colon\mathbb{F}_2^n\to\mathbb{F}_2\). The differential table of \(F\) records, for each \(a \in \mathbb{F}_2^n\) and \(b \in \mathbb{F}_2^n\), \[\delta_F(a,b) = |\{x \in \mathbb{F}_2^n: F(x\oplus a)\oplus F(x) = b\}|.\] The differential uniformity is \(\delta_F = \max_{a\neq 0,\,b} \delta_F(a,b)\).

The value 2 is the minimum possible for functions on \(\mathbb{F}_2^n\) with \(n>1\) [5]: for every nonzero \(a\), the derivative \(D_aF(x) = F(x\oplus a)\oplus F(x)\) is a map \(\mathbb{F}_2^n\to\mathbb{F}_2^n\) and \(\sum_b\delta_F(a,b) = 2^n\), so the maximum is at most \(2^n\) and at least 2.

2.2 Algebraic normal form and quadratic functions↩︎

Every \(f_j\colon\mathbb{F}_2^n\to\mathbb{F}_2\) has a unique algebraic normal form (ANF), a multilinear polynomial over \(\mathbb{F}_2\). A function \(F\) is quadratic if every \(f_j\) has ANF degree at most 2: \[f_j(x_0,\ldots,x_{n-1}) = c_{0,j} \oplus \bigoplus_{0\leq k < n} c_{\{k\},j}\,x_k \oplus \bigoplus_{0\leq p < q \leq n-1} c_{(p,q),j}\,x_p x_q.\] The \(n\binom{n}{2} = 8\cdot28 = 224\) quadratic coefficients \(c_{(p,q),j}\in\mathbb{F}_2\) completely determine the quadratic part of \(F\) (the linear part does not affect the APN property if \(F(0)=0\) is assumed by normalisation). We index them as \(c_1,\ldots,c_{224}\) in lexicographic order of \(((p,q),j)\), forming a vector \(\mathbf{c}\in\mathbb{F}_2^{224}\).

The pair index is \[\mathrm{idx}(p,q) = \frac{p(2n-p-1)}{2} + q - p - 1 \in \{0,\ldots,27\},\] and the full index is \(c_{\mathrm{idx}(p,q)\cdot n + j + 1}\).

2.3 CCZ- and EA-equivalence↩︎

EA-equivalence implies CCZ-equivalence. The converse fails in general but holds for quadratic APN functions:

CCZ-equivalence preserves the APN property, differential uniformity, and the structure of the difference table. Two functions in the same CCZ-class are cryptographically interchangeable as S-boxes.

2.4 Walsh transform and differential spectrum↩︎

For a VBF \(F\colon\mathbb{F}_2^n\to\mathbb{F}_2^n\), the Walsh transform is \[\widehat{F}(a,b) = \sum_{x\in\mathbb{F}_2^n}(-1)^{b\cdot F(x)\oplus a\cdot x}, \quad a,b\in\mathbb{F}_2^n,\] where \(\cdot\) is the standard inner product on \(\mathbb{F}_2^n\). The absolute Walsh spectrum is the multiset \(\{|\widehat{F}(a,b)| : a,b\in\mathbb{F}_2^n\}\).

The differential spectrum of \(F\) is the multiset \(\{\delta_F(a,b) : a\neq0, b\in\mathbb{F}_2^n\}\). All quadratic APN functions over \(\mathbb{F}_{2^{8}}\) have the same differential spectrum \(\{2^{2^7}\cdot 0, (2^8-1)2^7\cdot 2\}\) and the same Walsh spectrum; these are therefore not distinguishing invariants for quadratic APN over \(\mathbb{F}_{2^{8}}\).

2.5 Ortho-derivative↩︎

The ortho-derivative [2] is a more refined invariant.

3 Self-Equivalence Subspaces and the BBL2021 Framework↩︎

3.1 Definition and basic properties↩︎

\(V_A\) is a \(\mathbb{F}_2\)-linear subspace of the 224-dimensional space of all quadratic functions. Explicitly, \(F\circ A = A\circ F\) expands (via the ANF) to a homogeneous linear system over \(\mathbb{F}_2\) in the 224 coefficients of \(F\), so \(V_A = \ker(M_A)\) for a certain \(\mathbb{F}_2\)-linear map \(M_A\) (see Section 4 for the precise construction).

3.2 The BBL2021 classification for \(n=8\)↩︎

Beierle, Brinkmann, and Leander [1] classified all conjugacy classes of linear automorphisms of \(\mathbb{F}_{2^{8}}\) with small order and, for each, studied \(V_A\) by a recursive tree search. Table I of [1] lists only those class indices for which at least one APN function was found; absence from the table means the search found zero solutions. The class associated with the automorphism \(A\) studied here (order 5, described precisely in Section 4) does not appear in Table I.

3.3 Why random global search misses \(V_A\)↩︎

For the specific \(A\) studied here, \(\dim(V_A) = 40\). The fraction of all quadratic functions lying in \(V_A\) is \(2^{40}/2^{224} = 2^{-184}\approx10^{-55}\). No random global search (including the Beierle 2025 search [3]) can reach \(V_A\) in practice.

Within \(V_A\), the APN density is approximately \(10^{-4}\) (see Section 4.4). Any targeted search that samples uniformly from \(V_A\) will find APN functions at this rate, making the problem tractable.

4 The Self-Equivalence Subspace \(V_A\)↩︎

4.1 The automorphism \(A\)↩︎

Let \(C(p)\) denote the companion matrix of a polynomial \(p\in\mathbb{F}_2[X]\) (the transpose of the standard companion form). The polynomial \(q = X^4 + X^3 + X^2 + X + 1\) is the minimal polynomial of a primitive element of \(\mathbb{F}_{2^{4}}\) over \(\mathbb{F}_2\) (the fifth cyclotomic polynomial \(\Phi_5\)). Its companion matrix \(C(q)\in\mathrm{GL}(4,\mathbb{F}_2)\) has order 5: \[C(q) = \begin{pmatrix} 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 1 \end{pmatrix}.\]

We define \[A = \mathrm{block\_diag}(C(q), C(q))\in\mathrm{GL}(8,\mathbb{F}_2),\] where \(\mathrm{block\_diag}\) denotes the block-diagonal matrix \(\bigl(\begin{smallmatrix}C(q)&0\\0&C(q)\end{smallmatrix}\bigr)\). Since \(C(q)\) has order 5, \(A\) also has order 5. Concretely, \[A = \begin{pmatrix} 0&0&0&1&0&0&0&0\\ 1&0&0&1&0&0&0&0\\ 0&1&0&1&0&0&0&0\\ 0&0&1&1&0&0&0&0\\ 0&0&0&0&0&0&0&1\\ 0&0&0&0&1&0&0&1\\ 0&0&0&0&0&1&0&1\\ 0&0&0&0&0&0&1&1 \end{pmatrix}. \label{eq:A}\tag{1}\] One verifies \(A^5 = I_8\) and \(A^k\neq I_8\) for \(k\in\{1,2,3,4\}\), confirming \(\mathrm{ord}(A) = 5\). This is the linear automorphism whose self-equivalence subspace we study, corresponding to class index 22 in [1]. The matrix and the RREF of \(M_A\) (pivot/free indices) are included in the repository (file data/class22_basis.json).

4.2 Constructing the linear system↩︎

The condition \(F\circ A = A\circ F\) for a quadratic function \(F\) with ANF coefficient vector \(\mathbf{c}\in\mathbb{F}_2^{224}\) expands as follows. For each output bit \(j\in\{0,\ldots,7\}\) and each quadratic pair \((r,s)\) with \(r < s\), the entry \([A^TFA - A^TFA]_{j,(r,s)} = 0\) yields a homogeneous linear equation in \(\mathbf{c}\) over \(\mathbb{F}_2\).

Explicitly, the coefficient of \(c_{(p,q),k}\) in the equation for \((j,r,s)\) is \[\bigl[A_{p,r}A_{q,s} \oplus A_{p,s}A_{q,r}\bigr] \oplus A_{j,k}, \label{eq:system}\tag{2}\] summed over appropriate indices. Assembling all \((j,r,s)\) pairs gives a matrix \(M_A\in\mathbb{F}_2^{m\times 224}\) (where \(m\) is the number of equations) such that \(V_A = \ker(M_A)\).

4.3 RREF and the fc22\(+\)sol22 parameterization↩︎

Computing the reduced row echelon form (RREF) of \(M_A\) over \(\mathbb{F}_2\) yields:

• 184 pivot variables: indices \(\mathcal{P} = \{p_1,\ldots,p_{184}\}\subset\{1,\ldots,224\}\).

• 40 free variables: indices \(\mathcal{F} = \{f_1,\ldots,f_{40}\}\subset\{1,\ldots,224\}\), confirming \(\dim(V_A) = 40\).

• For each pivot \(p_i\), a linear combination \(p_i = \bigoplus_{f_j\in S_i}f_j\) over \(\mathbb{F}_2\), where \(S_i\subseteq\mathcal{F}\).

Any point in \(V_A\) is uniquely parameterized by a vector \(\mathbf{b} = (b_1,\ldots,b_{40})\in\mathbb{F}_2^{40}\) via: \[\begin{align} c_{f_k} &= b_k, \quad k = 1,\ldots,40, \tag{3}\\ c_{p_i} &= \bigoplus_{j: f_j\in S_i} b_j, \quad i = 1,\ldots,184. \tag{4} \end{align}\] We call 3 –4 the fc22\(+\)sol22 projection. Any \(\mathbf{b}\in\mathbb{F}_2^{40}\) maps to a valid point in \(V_A\); the mapping is bijective. This is the critical computational primitive for Phase 1.

4.4 APN density in \(V_A\)↩︎

To estimate the density of APN functions in \(V_A\), we sampled \(10^8\) random points via fc22\(+\)sol22 and found approximately \(10^4\) APN functions, giving an empirical density of \(\approx 10^{-4}\), or equivalently \(\approx 1\) APN per \(11\,000\) random points.

4.5 The NL = 4 slice structure↩︎

Of the 40 free variable indices in \(\mathcal{F}\), some correspond to pairs \((p,q)\) with \(p \leq 4\) (center bits) and others to pairs with \(p > 4\) (intra-slice bits). For \(n = 8\), the pairs with \(p > 4\) are exactly \((5,6)\), \((5,7)\), and \((6,7)\), contributing \(3\times 8 = 24\) intra-slice bits. The remaining \(40 - 24 = 16\) free bits are center bits.

This induces a partition of \(V_A\) into \(2^{16} = 65\,536\) non-overlapping affine hyperplanes of dimension 24, one for each choice of the 16 center bits. Each such hyperplane is called an NL = 4 slice, and the center bits determine the slice. Conversely, fixing the 200 ANF coefficients corresponding to pairs with \(p \leq 4\) (the NL = 4 normalization) selects a unique slice.

5 The Search Pipeline↩︎

The search proceeds in five stages, summarized in Figure 1.

5.1 Phase 1: APN-center search↩︎

The APN check (line 8, condition \(\delta_{F_\mathbf{b}} = 2\)) evaluates the difference table and runs in \(O(N^2)\) operations. In Python 3.11 with NumPy, one check takes \({\approx}0.6\) ms on one core of an AMD Ryzen 7840HS, corresponding to \({\approx}6\times10^6\) function evaluations per hour.

The slice signature (line 6) is the 16-bit string formed by the center bits of \(\mathbf{b}\). Deduplication on signatures ensures that no two Phase 1 invocations generate Magma files for the same slice (important for systematic exhaustion of \(V_A\)). There are \(2^{16} = 65\,536\) distinct signatures.

5.2 Phase 2: Gröbner basis slice computation↩︎

5.2.1 Normalization↩︎

Given APN center \(F_0\) with coefficient vector \(\mathbf{c}^{(0)}\), the NL = 4 normalization fixes the 200 ANF coefficients corresponding to pairs \((p,q,j)\) with \(p\leq4\) to their values in \(\mathbf{c}^{(0)}\). The 24 remaining coefficients (pairs \((5,6)\), \((5,7)\), \((6,7)\), each with 8 output bits \(j=0,\ldots,7\)) become unknowns \(c_1,\ldots,c_{24}\).

Since \(F_0\in V_A\) by construction, the normalization equations are consistent with the self-equivalence constraint. Any APN function in the same NL = 4 slice as \(F_0\) also lies in \(V_A\).

5.2.2 APN polynomial system↩︎

For each \(a\in\{1,\ldots,255\}\), define the derivative matrix \(M_a\in\mathbb{F}_2^{n\times n}\) by \[[M_a]_{j,k} = c_{(\min(k,\cdot),\max(k,\cdot)),j}\big|_{\text{free}},\] which encodes the linear part of \(D_aF\). The function \(F\) is APN at difference \(a\) if and only if \(\mathrm{rank}(M_a) = n-1 = 7\).

In practice, we build the polynomial \[p_a = \det(M_a) \cdot \prod_{r,c}\bigl(1 + \det(M_a^{(r,c)})\bigr)\] of degree \(\leq n = 8\) in \(c_1,\ldots,c_{24}\); vanishing of \(p_a\) is equivalent to \(\mathrm{rank}(M_a) = 7\).

The APN ideal is \[I = \langle p_1,\ldots,p_{255}\rangle + \langle c_i + \varepsilon_i : (p,q,j)\text{ fixed}\rangle \subset \mathbb{B}_{224},\] where \(\mathbb{B}_{224} = \mathbb{F}_2[c_1,\ldots,c_{224}]/(c_i^2+c_i)\) is the Boolean polynomial ring with grevlex ordering, and \(\varepsilon_i\in\mathbb{F}_2\) is the fixed value of each normalized coefficient.

Every element of \(\mathrm{Var}(I)\) corresponds to a function with \(\mathrm{rank}(M_a) = 7\) for all \(a \neq 0\), i.e., an APN function. Thus, when \(\dim(I)=0\), all solutions are APN by construction (no post-hoc verification is needed for the APN property itself, though we perform it as a sanity check).

5.2.3 Computational cost↩︎

The build phase (constructing all 255 polynomials \(p_a\)) takes \({\approx}8\)–15 minutes per slice on one core of an AMD Ryzen 7840HS. This dominates the total cost. Once the ideal is built, Gröbner basis computation reduces to Gaussian elimination over \(\mathbb{F}_2\) (since all variables are Boolean), and the GB solve time is typically under \(0.1\) s.

If \(\dim(I) = 0\): \(\mathrm{Var}(I)\) is finite; each element is a setting of the 24 free variables for which \(F\) is APN. All solutions are APN by construction.

If \(\dim(I) = -1\): the system is infeasible; the slice contains no APN functions. Empirically, \(\approx83\%\) of slices return \(\dim=-1\).

5.3 Phase 3: Collection and deduplication↩︎

Magma writes each found APN S-box to a result file via PrintFile. A Python collector script reads all result files, parses the lookup tables, computes SHA-256 hashes, and removes duplicates. A function found in multiple slices (e.g., as a center in one and a neighbor in another) appears only once in the deduplicated output.

The collector also handles Magma’s automatic line-wrapping (backslash\(+\)CRLF continuation for lines exceeding 80 characters).

5.4 Phase 4: CCZ classification↩︎

For each unique S-box \(F\), the ortho-derivative \(\pi_F\) is computed using the sboxU library [18]. The signature \(\sigma(F) = (\Delta\text{-spectrum}(\pi_F), W\text{-spectrum}(\pi_F))\) is computed and stored. Functions with equal signatures are grouped into the same CCZ-class (Theorem 2). Functions in different signature groups are in different CCZ-classes.

5.5 Phase 5: Verification against known databases↩︎

For each reference database \(\mathcal{K}\), we stream through all entries \(G\in\mathcal{K}\), compute \(\sigma(G)\), and check whether \(\sigma(G)\) matches any \(\sigma(F)\) for \(F\) in our collection. A match triggers an optional exact CCZ test via sboxU.are_ccz_equivalent_from_code. No match after exhausting \(\mathcal{K}\) certifies CCZ-inequivalence of all our functions from all functions in \(\mathcal{K}\).

5.6 Multi-slice extension↩︎

The NL = 4 normalization is one particular partition of \(V_A\) into \(2^{16}\) hyperplanes. One can instead fix any 200 of the 224 ANF coefficients, leaving any 24 free. Given a fixed APN center \(F_0\), there are \(\binom{40}{24} = 62\,852\,101\,650\) distinct choices of which 24 free variables to keep free. Each choice defines a different 24-dimensional hyperplane through \(F_0\), potentially yielding different neighbors.

In practice, one could run several (e.g., 10–100) independent normalizations per APN center by randomly choosing the 24 free variables. This would increase coverage without requiring new Phase 1 computations.

6 Experimental Results↩︎

6.1 Experimental setup↩︎

Hardware. AMD Ryzen 7840HS, 8 physical cores, 16 GB RAM, Windows 11.

Software. Magma V2.28-9 [19] (Phase 2); Python 3.11 with NumPy (Phases 1, 3); SageMath 10.x with sboxU [18] (Phases 4, 5) running under WSL 2 (Ubuntu 26.04).

Parallelism. Phase 2 uses a PowerShell queue script running up to 8 Magma processes simultaneously. Phase 1 runs single-threaded (trivially parallelizable).

Note on Magma exit code. Magma V2.20+ returns exit code 1 when quit; is called even on successful completion. The queue script was adjusted to treat exit codes 0 and 1 both as success, with any higher code indicating a genuine error.

6.2 Search coverage↩︎

Table 1 summarises the search coverage. All 428 executed slices are APN-positive (i.e., contain at least one APN function), because every slice center was itself found by Phase 1 to be APN. The absence of empty slices is therefore by construction, not unexpected. The key observation is that each of the 428 centers lies in \(V_A\), whereas 138 of the 566 found functions lie strictly outside \(V_A\); these are the functions attributable to the Gröbner basis step (see Section 6.5).

6.3 Property verification of collected functions↩︎

All 566 collected functions were independently verified for:

• Differential uniformity \(\delta_F = 2\): exhaustive difference table computation.

• Algebraic degree 2: Moebius transform of the ANF over \(\mathbb{F}_2\); maximum monomial degree is 2 for all 566 functions.

• Non-permutation: \(F\) is not bijective (the image \(|F(\mathbb{F}_{2^{8}})| < 256\) in all cases).

6.4 CCZ-classification↩︎

Table 2 shows the six CCZ-classes obtained by grouping on \(\sigma(F)\).

6.4.1 CLASS-A dominance↩︎

CLASS-A accounts for 362 of 566 functions (64%). This dominance was observed from the first experiments: 17 of the initial 25 functions belonged to the CLASS-A equivalent. Whether this reflects a larger volume of CLASS-A in \(V_A\), a higher local APN density, or both is not clear from the current data.

6.4.2 Gold functions as controls↩︎

CLASS-E and CLASS-F are identified by their \(\sigma\)-values matching those of Gold \(x^3\) (equivalently \(x^{129}\), by Frobenius orbit) and Gold \(x^9\) (equivalently \(x^{33}\)), respectively. The Gold functions \[x^{2^i+1},\quad \gcd(i,8)=1,\quad i\in\{1,3,5,7\},\] form exactly two CCZ-classes for \(n=8\): the orbit \(\{i=1,7\}\) giving \(x^3\) and \(x^{129}\), and the orbit \(\{i=3,5\}\) giving \(x^9\) and \(x^{33}\). No other quadratic APN power functions exist for \(n=8\) (Kasami, Welch, Niho, Dobbertin functions all have degree \(> 2\) for \(n=8\)).

The presence of Gold CCZ-equivalents in \(V_A\) is expected: \(V_A\) is a structured subspace that can intersect any CCZ-class. Their correct identification confirms that the pipeline accurately recovers known classes.

6.5 Slice-type structure and V\(_A\)-membership↩︎

A membership test — checking whether each found function \(F\) satisfies \(F(Ax) = A(F(x))\) for all \(x \in \mathbb{F}_{2^{8}}\) — reveals a precise partition of the 566 found functions:

• IN_CLASS22 (\(F \in V_A\)): 428 functions.

• OUT_CLASS22 (\(F \notin V_A\)): 138 functions.

Crucially, \(428 =\) the number of slice centers (one per slice), confirming that every APN center belongs to \(V_A\) by construction (Phase 1 searches within \(V_A\)), and that 138 functions were found by the Gröbner basis step outside \(V_A\).

Cross-referencing with CCZ-class reveals the following clean structure:

Table 3 shows that the 428 slices decompose into exactly three types:

The accounting is exact: \[\begin{align} 362 \times 1 &= 362 \quad (\text{CLASS-A}),\\ 36 \times 4 &= 144 \quad (36\text{ CLASS-E} + 72\text{ CLASS-B} + 36\text{ CLASS-C}),\\ 30 \times 2 &= \phantom{1}60 \quad (30\text{ CLASS-F} + 30\text{ CLASS-D}). \end{align}\] Total: \(362 + 144 + 60 = 566\) functions, \(428\) slices.

6.6 The V₂₂-CCZ landscape: a structural summary↩︎

The combined data from Sections 6.5 and Section [sec:verify] yield a precise and complete picture of the CCZ-class landscape in and around \(V_A\):

Figure 3 summarises the result. The exact-membership test (criterion: \(F(Ax) = A(F(x))\) for all \(x\in\mathbb{F}_{2^{8}}\), verified on all 566 functions) establishes the following theorem with complete computational proof:

6.7 The GB-navigation paradigm: a generalizable mechanism↩︎

Beyond the specific results, the experiment reveals a mechanism for systematically exploring the APN landscape that may be of independent methodological value. We describe it in abstract terms:

The mechanism (Figure 4) has three components:

6.8 Comparative experiments: the role of V_A-structure↩︎

To establish that the gateway phenomenon is genuinely tied to the self-equivalence subspace \(V_A\) rather than being a generic property of APN functions, we ran two additional experiments with centers drawn from outside \(V_A\).

6.8.0.1 Experiment B: Beierle 2025 dataset centers (batch_from_dataset).

We took the first 532 functions from the Beierle 2025 database new_apns.txt [3] as slice centers. These functions are guaranteed to be quadratic APN over \(\mathbb{F}_{2^{8}}\) but do not belong to \(V_A\). For each, we ran an NL=4 Gröbner basis computation using the same 24-variable polynomial system as in the main experiment. Result: every one of the 532 slices produced exactly 1 solution (the center itself, dim = 0). No function found a single neighbor.

6.8.0.2 Experiment C: random center candidates (batch_rnd).

We generated 20 random vectors in \(\mathbb{F}_2^{224}\) as center candidates (not from \(V_A\), and not verified to be APN). Result: every one of the 20 slices returned dim = \(-1\) (empty ideal), meaning no APN function exists in that NL=4 hyperplane at all. This is consistent with the global APN density \(\approx 2^{-n} \approx 4\times 10^{-3}\): a random hyperplane through a random point in the full space is very unlikely to contain any APN function.

The combined picture is shown in Table 4:

The contrast is sharp. In Experiment A, Gold-centered slices consistently yield 2–3 new APN functions per slice (100% hit rate, 66 out of 66 slices). In Experiments B and C, zero neighbors were found across 552 combined slices.

These results support a stronger conclusion than mere suggestion:

6.9 The three reference databases↩︎

We checked the four new classes (CLASS-A through CLASS-D) against three independent reference databases:

1. Beierle 2025 (new_apns.txt) [3]: 3 775 599 quadratic APN functions for \(n=8\) that are CCZ-inequivalent to the 32 892 previously known. The file does not include Gold or other classical functions.

2. BBL2021 pre-2020 (apn_8bit.txt) [1]: 12 921 quadratic APN functions found by the BBL2021 method, all CCZ-inequivalent to Gold.

3. Gold functions (direct): The two Gold CCZ-classes for \(n=8\) are verified directly via \(\sigma\)-comparison (not via a file).

6.10 Verification procedure and result↩︎

For each of the 500 functions in CLASS-A through CLASS-D, we compute \(\sigma(F)\) via sboxU and run the streaming comparison of Section 5.

6.11 Verification summary↩︎

Table 5 collects all verification checks.

7 Discussion↩︎

7.1 Interpretation of results↩︎

The main empirical findings are:

(i) A 40-dimensional structured subspace (\(V_A\)) that was declared empty by the BBL2021 tree search contains at least four CCZ-classes of quadratic APN functions, with 566 explicit representatives found in 0.65% of the search space (428 out of 65 536 NL = 4 slices).

(ii) The Gröbner basis step is essential for finding three of the four new classes. The \(V_A\)-membership analysis (Table 3) shows that CLASS-B, C, and D lie entirely outside \(V_A\): all 138 representatives are OUT_CLASS22. A Phase 1-only search within \(V_A\) recovers only CLASS-A and the Gold controls, missing classes B, C, and D completely. Equivalently, the Gröbner basis approach is both necessary and sufficient for their discovery.

(iii) \(V_A\) contains exactly three distinct CCZ-classes in our sample: CLASS-A (a new class native to \(V_A\)) and the two Gold classes (E and F). Gold functions act as “gateways”: their NL = 4 slices reach outside \(V_A\) and expose three new external classes. Specifically, each Gold-\(x^3\) slice yields two CLASS-B functions and one CLASS-C function in addition to the center; each Gold-\(x^9\) slice yields one CLASS-D function in addition to the center. This structure is exact and uniform across all 428 explored slices (Theorem 4).

(iv) The Beierle 2025 database [3], while containing 3.8 million functions, does not cover \(V_A\): the probability of any specific entry lying in \(V_A\) is \(\approx 10^{-55}\). More concretely, the comparative experiments (Section 6.8, Table 4) show that 532 Beierle 2025 functions used as slice centers produced zero APN neighbors, while 66 Gold-function centers (inside \(V_A\)) each produced 1–3 neighbors. In all 552 experiments outside \(V_A\), not a single new APN function was found. This is an experimental fact, not a conjecture.

(v) CLASS-A’s dominance (362 of 566 functions, 64%) is partly explained by the slice structure: all 362 CLASS-A-centered slices contain only the center itself (Type I slices), so CLASS-A’s count equals the number of CLASS-A seeds used in Phase 1. This dominance reflects the distribution of seeds, not necessarily the overall abundance of CLASS-A in \(V_A\).

(vi) The experiment reveals a general GB-navigation mechanism (Section 6.7, Figure 4) for exploring the APN landscape. Rather than searching exhaustively in a large space, one: (a) finds a structured subspace with manageable APN density; (b) uses the Gröbner basis to enumerate APN neighbors of any center, including those outside the subspace; (c) identifies gateway functions whose neighbors belong to new CCZ-classes. This mechanism transforms the APN search problem from “find APN functions in a huge space” to “navigate locally around known structured objects.” Crucially, the comparative experiments (Table 4) establish that this mechanism is specific to the self-equivalence structure: not any APN function is a productive seed — only those lying in \(V_A\). Among 3.8 million Beierle 2025 functions used as centers (Exp. B, 532 slices) and 20 random candidates (Exp. C), zero APN neighbors were found. This is an experimentally established fact, not a conjecture. The approach is qualitatively different from all prior construction methods (bivariate, isotopic shift, bent extensions), which build APN functions from algebraic templates rather than exploring local neighborhoods of structured-subspace functions.

7.2 Limitations↩︎

We list the limitations honestly:

1. Coverage. Only 0.65% of the total \(65\,536\) NL = 4 slices were computed. The four observed new CCZ-classes may not be all classes present in or accessible from \(V_A\).

2. No algebraic characterization. The new CCZ-classes are defined purely computationally by their \(\sigma\)-values and S-box tables. We have no polynomial or structural description.

3. Single subspace. Only one BBL2021 class (the one with automorphism of order 5) was studied. The approach applies to all BBL2021 classes [1]; other subspaces may yield additional new classes.

4. Exact CCZ tests. Classification relies on Theorem 2, which is theoretically exact for quadratic APN. We did not run exhaustive exact pairwise CCZ tests (are_ccz_equivalent_from_code) within each class, as these are not necessary given the theorem. For a fully self-contained proof independent of [2], one could add such tests.

5. Claim scope. Theorem 6 certifies inequivalence with the databases checked. Whether the four new classes are CCZ-inequivalent to all quadratic APN functions for \(n=8\) would require checking against a complete database, which does not yet exist.

7.3 Scalability↩︎

The pipeline scales in the following ways:

Exhausting \(V_A\). The remaining \(65\,536 - 428 = 65\,108\) NL = 4 slices can be computed to enumerate all APN functions in \(V_A\) under this normalization. Estimated time: \(65\,108\) slices \(\times 8\) min/slice \(\div\) 8 cores \(\approx 54\) wall-clock days.

Multi-slice normalization. As noted in Section 5.6, each APN center supports \(\binom{40}{24}\approx6\times10^{10}\) normalizations. Running 10–100 per center increases coverage per center without extra Phase 1 cost.

Larger free dimension. Fixing fewer than 200 coefficients in Phase 2 leaves more than 24 free variables. For example, fixing 176 coefficients (NL = 3) leaves 48 free variables. The GB build time grows but so does the number of APN neighbors found per slice. On a cluster with many cores, this trade-off could be favorable.

Other BBL2021 classes. For each BBL2021 class index \(k\) with \(\dim(V_{A_k}) = d_k\), the pipeline applies with \(d_k - 24\) center bits and \(2^{d_k-24}\) total slices. Classes with \(d_k \leq 60\) are directly amenable to the current pipeline.

8 Related Work↩︎

8.1 Construction approaches for quadratic APN functions↩︎

Quadratic APN functions over \(\mathbb{F}_{2^{n}}\) have been found by many different methods. The bivariate construction of Edel and Pott [11] (building \(F(x,y)\) as a polynomial in two variables over \(\mathbb{F}_{2^{n/2}}\)) yielded 8 157 new instances for \(n=8\). Yu, Wang, and Li [12] used a matrix representation approach. Budaghyan, Calderini, Carlet, Coulter, and Villa [15] introduced isotopic shifts, constructing new APN functions by modifying Gold functions; the generalized version of this method [17] (trims and extensions) produced an additional 6 368 CCZ-classes for \(n=8\). Beierle, Brinkmann, and Leander [1] studied LE-automorphism groups and found 12 921 new instances; Beierle and Leander [16] adapted this to find further instances by extending \(n=7\) APN functions. Beierle et al. [3] use extensions of vectorial \((n/2,n/2)\)-bent functions to construct the current largest known set (\(\approx 3.8\) million CCZ-classes).

Our approach is methodologically orthogonal to all of the above: rather than constructing functions from algebraic templates (power functions, bivariate polynomials, isotopic shifts, bent function extensions), we perform an exhaustive Gröbner basis enumeration within a fixed 24-dimensional affine hyperplane of the self-equivalence subspace \(V_A\). This targets a region that is unreachable by both global random search and the BBL2021 tree method, and finds functions that do not arise from any known construction.

8.2 CCZ-equivalence testing: EA-tests for quadratic APN functions↩︎

Several recent works address the problem of testing or recovering EA- and CCZ-equivalence efficiently. Calderini [20] classified all known APN functions in small dimensions up to EA-equivalence, providing a reference for EA-class structure. Canteaut, Couvreur, and Perrin [2] gave an efficient algorithm for EA-recovery for quadratic functions (complexity \(O(n \cdot 2^{2n})\), superseding all prior methods).

For our purposes, the ortho-derivative provides a complete CCZ-invariant in the quadratic case (Theorem 2 and Remark 1), making a separate EA-test unnecessary for classification. The sboxU library [18] implements both the ortho-derivative computation and the differential/Walsh spectra that constitute the signature \(\sigma(F)\). The completeness of \(\sigma\) for quadratic APN functions was verified by running sboxU’s are_ccz_equivalent_from_code on all pairs within each putative class (three classes of size \(\leq 36\)) and confirming that all pairs returned True; this provides an independent validation of the classification beyond the theoretical guarantee.

8.3 Gröbner basis methods in cryptography↩︎

Gröbner bases over Boolean polynomial rings have been used in algebraic cryptanalysis since the work of Courtois and Pieprzyk (XL algorithm [21]) and extended by Albrecht and Bard [22] (M4RI library). Their use in APN function search is less common; to our knowledge this is the first published application of GB slicing to systematic APN search in a self-equivalence subspace. The key feature that makes GB tractable here is the reduction to 24 Boolean unknowns via NL = 4 normalization: the GB solve is essentially Gaussian elimination (all variables are in \(\mathbb{F}_{2^{2}}\)) and takes \(< 0.1\) s per slice once the ideal is built.

8.4 Open problem: APN permutations in even dimension↩︎

Browning, Dillon, McQuistan, and Wolfe [7] found an APN permutation in dimension 6 by CCZ-transforming a quadratic APN function. Whether APN permutations exist in dimension 8 remains open; it is one of the central open problems in the area [23]. The functions found in this paper are non-permutations; whether any of the new CCZ-classes contains an APN permutation in its CCZ-equivalence class is an open question that the explicit S-box tables in Appendix 10 make accessible to future investigation.

9 Conclusion↩︎

We have shown that the self-equivalence subspace \(V_A\) associated with the order-5 automorphism (BBL2021 class index 22), dimension 40, which was previously reported to contain no APN functions under recursive tree search, contains at least four CCZ-classes of quadratic APN functions for \(n=8\). These four classes (500 functions out of 566 total) are absent from the Beierle 2025 database and from all pre-2020 known instances, as certified by the exact ortho-derivative invariant.

The result was obtained using a two-phase pipeline: random sampling in \(V_A\) via an explicit RREF parameterization (Phase 1), followed by Gröbner basis slice computation in Magma (Phase 2). The pipeline covered 428 of the \(65\,536\) NL = 4 slices (0.65%); two Gold CCZ-equivalents serve as positive controls and as gateways to the new external classes.

Structural result. A structural decomposition of the 428 explored slices into three types yields a complete and self-consistent account of all 566 found functions (Table 3 and Theorem 4): the slices are either CLASS-A-centered (Type I, 362 slices, 1 APN each) or Gold-centered (Types II–III, 66 slices, 2–4 APN each). Critically, all three purely new classes (B, C, D) lie entirely outside \(V_A\) and appear exclusively in Gold-centered slices. This confirms that the Gröbner basis step is essential: a Phase 1-only search within \(V_A\) would recover CLASS-A and the Gold controls but miss CLASS-B, C, and D entirely. It also implies that Gold functions are sufficient Phase 1 seeds for rediscovering all new non-\(V_A\) classes found here.

Methodological contribution. Beyond the specific new CCZ-classes, the experiment reveals a general GB-navigation mechanism (Section 6.7) that addresses a fundamental challenge in APN search: the near-zero global density \((\approx10^{-55})\) makes random exploration of the full \(2^{224}\)-dimensional space hopeless. The mechanism provides a structured alternative: (a) exploit a self-equivalence subspace for tractable Phase 1 seeding; (b) use the Gröbner basis to enumerate all APN neighbors locally, including those outside the subspace; (c) identify gateway functions whose neighborhood contains new CCZ-classes.

The comparative experiments (Section 6.8, Table 4) establish this mechanism as specific to the self-equivalence structure, not a generic property of APN functions. Among 532 functions from the Beierle 2025 database used as GB-slice centers and 20 random candidate centers (552 slices total), zero APN neighbors were found. In sharp contrast, all 66 Gold-function centers inside \(V_A\) produced 1–3 new APN neighbors each. This difference is experimentally established with zero exceptions across 618 total slice computations. The central empirical conclusion is: not any APN function is a productive GB-seed; only APN functions lying in a self-equivalence subspace exhibit the gateway phenomenon that exposes neighboring CCZ-classes. Among 3.8 million currently known APN functions, none appears to be a useful seed via this method. The Gold functions inside \(V_A\), despite being “known,” are uniquely productive.

Further observations. Five observations are of independent interest: (1) The fc22\(+\)sol22 projection gives an efficient and exact parameterization of any self-equivalence subspace, enabling targeted sampling at APN density \(\approx10^{-4}\) instead of \(\approx10^{-55}\). (2) The NL = 4 normalization reduces the Gröbner basis problem to 24 Boolean unknowns with a polynomial build time of \(\approx10\) minutes and an essentially immediate solve. (3) The ortho-derivative signature provides a complete, non-heuristic CCZ-invariant for the quadratic case, allowing rigorous classification and database comparison without exhaustive pairwise testing. (4) The slice-type structure reveals a nontrivial algebraic relationship between \(V_A\) and the CCZ-class geometry: Gold functions occupy a boundary role, “bridging” \(V_A\) and classes external to it. (5) All four new classes (A, B, C, D) are definitively absent from the largest known compilation (3.8 million functions), yet they are easily accessible via the GB-navigation pipeline in under 10 CPU-hours.

The pipeline generalizes to all BBL2021 class indices and to larger dimensions, providing a systematic program for extending the classification of quadratic APN functions beyond currently known databases.

10 Class Representatives and Signature Data↩︎

Table 6 gives the full ortho-derivative signature hashes and one representative S-box per class. Full lookup tables (256 decimal values each) are given in Table [tab:reps].

[!t]
\renewcommand{\arraystretch}{1.15}
\caption{Representative S-boxes: $F(0),F(1),\ldots,F(255)$ in decimal}
\label{tab:reps}
\centering\small\setlength\tabcolsep{3pt}
\begin{tabular}{lp{15cm}}
\toprule
\textbf{Class} & \textbf{Lookup table}\\
\midrule
CLASS-A &
0,0,0,63,0,42,97,116,0,253,75,137,194,21,232,0,0,24,207,232,165,151,11,6,183,82,51,233,208,31,53,197,
0,122,47,106,113,33,63,80,186,61,222,102,9,164,12,158,
19,113,243,174,199,143,70,49,30,129,181,21,8,189,194,72,
0,211,235,7,65,184,203,13,226,204,66,83,97,101,160,155,
178,121,150,98,86,183,19,205,231,209,136,129,193,221,207,236,
144,57,84,194,160,35,5,185,200,156,71,44,58,68,212,149,
49,128,58,180,164,63,206,106,222,146,158,237,137,239,168,241,
0,217,86,176,57,202,14,194,130,166,159,132,121,119,5,52,
95,158,198,56,195,40,59,239,106,86,184,187,52,34,135,174,
34,129,91,199,106,227,114,196,26,68,40,73,144,228,195,136,
110,213,216,92,131,18,84,250,225,167,28,101,206,162,82,1,
76,70,241,196,52,20,232,247,44,219,218,18,150,75,1,227,
161,179,211,254,124,68,111,104,118,153,79,159,105,172,49,203,
254,142,108,35,247,173,4,97,36,169,253,79,239,72,87,207,
0,104,93,10,172,238,144,237,109,248,123,209,3,188,116,244\\[2pt]
CLASS-B &
0,0,0,63,0,42,97,116,0,253,75,137,194,21,232,0,0,24,207,232,165,151,11,6,183,82,51,233,208,31,53,197,
0,122,47,106,113,33,63,80,186,61,222,102,9,164,12,158,
19,113,243,174,199,143,70,49,30,129,181,21,8,189,194,72,
0,211,235,7,65,184,203,13,226,204,66,83,97,101,160,155,
178,121,150,98,86,183,19,205,231,209,136,129,193,221,207,236,
38,143,226,116,22,149,179,15,126,42,241,154,140,242,98,35,
135,54,140,2,18,137,120,220,104,36,40,91,63,89,30,71,
0,217,86,176,57,202,14,194,130,166,159,132,121,119,5,52,
95,158,198,56,195,40,59,239,106,86,184,187,52,34,135,174,
148,55,237,113,220,85,196,114,172,242,158,255,38,82,117,62,
216,99,110,234,53,164,226,76,87,17,170,211,120,20,228,183,
76,70,241,196,52,20,232,247,44,219,218,18,150,75,1,227,
161,179,211,254,124,68,111,104,118,153,79,159,105,172,49,203,
254,142,108,35,247,173,4,97,36,169,253,79,239,72,87,207,
0,104,93,10,172,238,144,237,109,248,123,209,3,188,116,244\\[2pt]
CLASS-C &
0,0,0,247,0,175,30,70,0,212,177,146,35,88,140,0,0,105,168,54,8,206,190,143,186,7,163,233,145,131,150,115,
0,230,205,220,191,246,108,210,15,61,115,182,147,14,241,155,
57,182,92,36,142,174,245,34,140,215,88,244,24,236,210,209,
0,217,60,18,117,3,87,214,145,156,28,230,199,101,84,1,46,158,186,253,83,76,217,49,5,97,32,179,91,144,96,92,
8,55,249,49,194,82,45,74,150,125,214,202,127,59,33,146,
31,73,70,231,221,36,154,148,59,185,211,166,218,247,44,246,
0,2,93,168,103,202,36,126,234,60,6,39,174,215,92,210,
252,151,9,149,147,87,120,75,172,19,232,160,224,240,186,93,
38,194,182,165,254,181,112,204,195,243,226,37,56,167,7,111,
227,110,219,161,51,17,21,192,188,229,53,155,79,185,216,217,
197,30,164,136,215,163,168,43,190,177,110,150,143,47,65,22,
23,165,222,155,13,16,218,48,214,176,174,63,239,38,137,183,
235,214,71,141,70,212,244,145,159,118,130,156,17,87,18,163,
0,84,4,167,165,94,191,179,206,78,123,12,72,103,227,59\\[2pt]
CLASS-D &
0,0,0,215,0,63,94,182,0,116,97,194,163,232,156,0,0,212,235,232,79,164,250,198,227,67,105,30,15,144,219,147,
0,214,88,89,57,208,63,1,129,35,184,205,27,134,124,54,80,82,227,54,38,27,203,33,50,68,224,65,231,174,107,245,
0,108,92,231,175,252,173,41,109,117,80,159,97,70,2,242,96,216,215,184,128,7,105,57,238,34,56,35,173,94,37,1,
181,15,177,220,35,166,121,43,89,151,60,37,108,157,87,113,
133,235,106,211,92,13,237,107,138,144,4,201,240,213,32,210,
0,20,199,4,167,140,62,194,177,209,23,160,181,234,77,197,
224,32,204,219,8,247,122,82,178,6,255,156,249,114,234,182,
213,23,74,95,75,182,138,160,229,83,27,122,216,81,120,38,
101,115,17,208,180,157,158,96,182,212,163,22,196,153,143,5,
176,200,43,132,184,255,125,237,108,96,150,77,199,244,99,135,
48,156,64,59,119,228,89,29,15,215,30,17,235,12,164,148,
208,126,19,106,225,112,124,58,141,87,47,34,31,250,227,209,
0,122,40,133,126,59,8,154,190,176,247,46,99,82,116,146\\
\bottomrule
\end{tabular}

11 Computational Infrastructure and Reproducibility↩︎

The complete dataset (566 S-boxes in JSON format) and all source code are available at: https://github.com/KuznetsovKarazin/apn-gb-search; DOI: https://doi.org/10.5281/zenodo.20626047

Key scripts:

• src/find_apn_centers.py: Phase 1 APN-center search. Input: RREF data (fc22, sol22). Output: batch directory with Magma files and apn_centers.json.

• run_batch.ps1: PowerShell queue runner. Executes all *.m files in a directory using up to \(k\) parallel Magma processes (-MaxParallel\(k\)). Treats Magma exit codes 0 and 1 both as success.

• src/collect_results.py: Parses Magma output files, handles line-wrapping, deduplicates by SHA-256 hash.

• src/classify_apn.py: Computes ortho-derivative signatures via sboxU and groups functions into CCZ-classes.

• src/check_vs_beierle_db.py: Streaming comparison of our signatures against a reference database.

Requirements: Magma V2.20+ [19]; Python 3.9+ with NumPy (Phases 1, 3); SageMath with sboxU [18] (Phases 4, 5).

References↩︎

---

1. O. Kuznetsov is with the Department of Theoretical and Applied Sciences, eCampus University, Via Isimbardi 10, 22060 Novedrate (CO), Italy, and also with the Department of Intelligent Software Systems and Technologies, School of Computer Science and Artificial Intelligence, V.N. Karazin Kharkiv National University, 4 Svobody Sq., 61022 Kharkiv, Ukraine. E-mail: oleksandr.kuznetsov@uniecampus.it. ORCID: 0000-0003-2331-6326.↩︎

2. Dataset and source code: https://github.com/KuznetsovKarazin/apn-gb-search; DOI: https://doi.org/10.5281/zenodo.20626047↩︎

3. Manuscript received 2026-06-11.↩︎