A Framework for Fine-Grained Nonlinearity Optimization of Boolean and Vectorial Boolean Functions

Boolean functions and vectorial Boolean functions (S-boxes) are widely used cryptographic primitives for achieving cryptanalytic resistance of modern block or stream ciphers. In the aspect of information security, one of the most desirable characteristics a given S-box should possess is a high nonlinearity. In this paper, we project the nonlinearity optimization problem to the domain of binary integer programming. Then, we demonstrate how this interconnection could be successfully exploited by SAT solvers. The provided toolbox could serve in cases, where the designer’s goal is to increase (or intentionally decrease) the nonlinearity of a given S-box by applying as minimum changes as possible. For example, we demonstrate how the Skipjack S-box, developed by the U.S. National Security Agency (NSA), and the Kuznyechik S-box, developed by the Russian Federation’s standardization agency, could be optimized to a higher nonlinearity by tweaking, respectively, just 4 and 12 bits (out of 2048). In the end, we show that bijective (8,8) S-boxes, the eight coordinates of which possess the currently known optimal nonlinearity value of 116, do exist.


I. INTRODUCTION
S-boxes are the main building blocks of the nonlinearity layer of modern block or stream ciphers, like AES [1], SNOW [2], TWOFISH [3], WHIRLPOOL [4], PICARO [5], PRINCE [6]. A weak S-box, in a cryptographic perspective, can be exploited by various attacks such as linear cryptanalysis [7], differential cryptanalysis [8], boomerang attack [9], algebraic attacks [10], and others [11]. One of the most desirable characteristics a given S-box should possess is a high nonlinearity value. An S-box with optimal or near-optimal nonlinearity value can be constructed by using the finite field inversion method as shown in [12]. However, such an S-box would be closely related to an algebraic structure. Thus, as a proactive countermeasure to possible algebraic attacks, new ways of generation or optimization of pseudo-randomly generated S-boxes are proposed. The published methods are manifold. However, heuristic optimization of a given S-box The associate editor coordinating the review of this manuscript and approving it for publication was Jiafeng Xie. is a resource-consuming task. In general, the size of the search space of bijective S-boxes with dimension (n, n) is 2 n !. For dimension (8,8), which is a common choice among the majority of popular S-boxes, the size of the search space is approximately 2 1684 . Furthermore, during the search for a better S-box, each candidate is repeatedly evaluated by some complex cost function. For example, considering state-of-theart optimization algorithms [13]- [17], the Walsh-Hadamard transform matrix elements are regularly collected during each iteration. However, such calculations add significant overhead to the optimization routine.
In this work, we find an interconnection between the S-box nonlinearity optimization problem and binary integer programming. Thus, we create a lightweight optimization routine, which does not cause any significant computational burden. Moreover, the proposed toolbox could be used as a proof of infeasibility.
Another major drawback of the state-of-the-art heuristic techniques is their aggressiveness on the initial S-box. Hence, in most cases, it is difficult to link the resulted S-box with the initial S-box. In fact, it is difficult to prove that such a link exists in the first place. However, in this work, we present a fine-grained optimization routine, which allow us to optimize the nonlinearity value of a given S-box with minimum changes as possible. From the designer perspective, this property is particularly beneficial, since we could focus the optimization routine on the weak components of a given S-box, without degrading the remaining ones. We demonstrate the effectiveness of the proposed algorithm by increasing the nonlinearity of the Skipjack S-box, developed by NSA, and Kuznyechik S-box, developed by Russian Federation's standardization agency, by tweaking respectively 4 and 12 (out of 2048) bits only.
During our research, we collected more than 300 papers, which emphasize on the average coordinate nonlinearity value, or ACNV, of bijective S-boxes, instead of the nonlinearity value of the weakest component. In our previous work [18], we show that such constructions should be used with considerable caution. The currently known maximum nonlinearity value for 8-variable balanced Boolean functions is 116 [19]. Furthermore, as shown in [20], the nonlinearity value of 8-variable balanced Boolean functions is upper bounded by 120, which means that the maximum theoretical ACNV of (8,8) bijective S-boxes is less or equal to 118.0. If a bijective S-box with ACNV greater than 116 is found, at least one of its eight coordinates will possess a nonlinearity value of 118, which will finally give an answer to the long-standing problem of the maximum possible nonlinearity value for 8-variable balanced Boolean functions. However, there is an academic skepticism that 8-variable balanced Boolean functions with nonlinearity value 118 really exist. Having this in mind, one open question to be answered is: • Does bijective (8,8) S-box with an ACNV value of 116 really exist? In this paper, by using SAT solving techniques, we show that bijective (8,8) S-boxes with ACNV value of 116 really exist. However, despite our attempts, we were not able to find an 8-variable balanced Boolean function with a nonlinearity of 118.

II. PRELIMINARIES
In this section, we present some standard definitions and S-box characteristics. More detailed overview could be found, for example, in [21] and [22].

Definition 2 (Algebraic Normal Form -ANF f ):
The algebraic normal form of an n-variable Boolean function f (x), denoted by ANF f , is given by the following equation: ANF f = a 0 ⊕ a 1 x 1 ⊕ a 2 x 2 ⊕ · · · ⊕ a n x n ⊕ a 1,2 x 1 x 2 ⊕ · · · ⊕ a 1,2,··· ,n x 1 x 2 · · · x n , where the coefficients a belongs to B.
Definition 3 (Hamming Distance): The Hamming distance between two n-variable Boolean functions f (x) and g(x), denoted by d H (f , g), represents the number of differing elements in the corresponding positions of their truth tables.
Definition 4 (Linear Boolean Function): Any n-variable Boolean function of the form: where w, x ∈ B n , is called a linear function.
Definition 5 (Affine Boolean Function): Any n-variable Boolean function of the form: where w 0 ∈ B and w, x ∈ B n , is called an affine function.
Definition 6 (Walsh-Hadamard Transform -WHT): For an n-variable Boolean function f (x), represented by its polarity table [f (x)], the Walsh-Hadamard transform, or WHT, Table -S-Box: An n-binary input to m-binary output mapping S : B n ⇔ B m , which assigns some y = (y 1 , y 2 , · · · , y m ) ∈ B m by S(x) = y to each x = (x 1 , x 2 , · · · , x n ) ∈ B n , is called an (n, m) substitution table (S-box) and is denoted by S(n, m). Definition 8 (Bijective S-box): An S-box S(n, m) is said to be bijective, if it maps each input x ∈ B n to a distinct output y = S(x) ∈ B m and all possible 2 m outputs are present.
where w and v are arranged lexicographically in B n and respectively B m .
Definition 13 (Linear Approximation Table -LAT): The linear approximation table of an S-box S(n, m), denoted by LAT S or S LAT , is a table with 2 n rows and 2 m columns, which entries are given by: where Y is a consequent linear combination of coordinates of the current S-box and X is the consequent linear function with length n. Definition 14 (S-box Nonlinearity): The nonlinearity of an S-box S(n, m), denoted by S NL , is defined as: where {w i } is the set of all elements in LAT, except the uppermost left one.

Definition 15 (S-box ACNV):
The average coordinate nonlinearity value, or ACNV, of a given S-box S, is the average value of all nonlinearities of coordinates of S. Table -DLUT): Each S-box is uniquely defined by its LUT. Translating each row of the LUT as a decimal number uniquely defines the same S-box as a decimal look-up table (DLUT). A bijective S-box S(n, n) has exactly n coordinates. Each coordinate is defined by its truth table, which consists of 2 n elements from B. So, we have a total of n2 n elements from B which uniquely define S. Following this observation, we have a total of 2 n2 n possible choices of S-boxes. But not all of them are bijective. To further restrict our choices to bijective S-boxes only, we need to restrict the set of possible S-boxes with the following observation -all elements in the DLUT of S should be distinct. This reduces the possible choices of S-boxes from 2 n2 n to 2 n !.

III. DESCENSION THROUGH COUPLINGS
In this section, we introduce the concept of couplings, coordinate decomposition, degree of descendibility, S-box coordinate extended linear approximation table (CELAT), as well as some useful properties and inner relationships. For convenience, let us denote as f (n) i the integer extracted from n, by flipping its i-th bit of its binary representation. Obviously, Tweaking a bijective S-box S by flipping just one bit in its corresponding Lookup Table (LUT) will convert S to a non-bijective S-box.
Proof: [Proof of Lemma 1] We take an arbitrary bijective S-box S(n, n) and its corresponding Look-up S LUT and Decimal Look-up S DLUT tables. We pick the flipped bit to be somewhere inside the row with index i of S LUT . The resulted Look-up Table will be denoted as S LUT . We will prove that the S-box S is not bijective.
, then the resulted Decimal Look-up Table of where d i is the decimal integer, which corresponds to the bit-concatenation of all the bits from the i-th row of S LUT . However, from the bijectivity property it follows that ∀i : i = j → d i = d j . Furthermore, by definition (see Lemma 8), S DLUT is in fact a permutation of all the 2 n integers in the interval [0, 2 n − 1]. Since S is with dimensions (n, n), each element of S DLUT is represented by exactly 2 n bits. Having this in mind, the number of possible distinct values of d i is 2 n (d i with first bit flipped, d i with second bit flipped, . . . , or d i with last bit flipped). Since the binary representations of all those distinct values consist of exactly n bits, their decimal representations are values less or equal than 2 n − 1. Therefore, no matter which bit of d i is flipped, d i will collide with exactly one d j , for some j = i. Hence, S DLUT will hold two different elements, d i and d j , with equal values, and therefore, S is a non-bijective S-box.
It is not possible to get a bijective S-box by modifying (flipping) a single bit of the S LUT of other bijective S-box. However, as shown in the next Lemma, the minimum count of bits we need to change in the LUT of a random bijective S-box to get a new bijective S-box is 2.

Lemma 2 (Couplings Lemma):
The smallest nonzero number of bits from the LUT of a random bijective S-box needed to be modified to obtain another bijective S-box is 2.
Proof: [Proof of Lemma 2] Let us take a bijective S-box S(n, n) and define the DLUT of S as an array We recall that S DLUT is in fact a permutation of all the 2 n integers in the interval [0, 2 n − 1]. Without loss of generality, we pick the first flipped bit to be somewhere inside the row with index i of S LUT . Let us denote the resulted Look-up Table, when this bit is flipped, as S LUT . As shown in the previous lemma, the S-box S , which corresponds to the Look-up Table S LUT , is not bijective. However, we will show that we could always flip another distinct (non-trivial) bit, which could transform the S-box S to some bijective S-box S , where S = S.
Using the notations introduced throughout the proof of the previous lemma, we have where d i is the decimal integer, which corresponds to the bit-concatenation of all the bits from the i-th row of S LUT . Following the same observation made in Lemma 1, no matter which bit of d i is flipped, d i will collide with exactly one d j , for some j = i, and S DLUT will hold two different elements, d i and d j , with equal values However, all the remaining 2 n −2 elements from S DLUT , i.e. all the elements in S DLUT with d j and d i excluded, differ from each other. Since d j = d i , and d i = d i , we could highlight two reasons for the non-bijectivity of the S-box S : Having this in mind, if we could modify d j to d i , by using just a single flip, we could convert S DLUT to S DLUT , where the S-box S , which corresponds to the Decimal Look-up Table S DLUT , is bijective. It is trivial to be shown that this modification is possible. Let us recall that d i is created by flipping a single bit in d i on position, for example, x. Therefore, since d j = d i , flipping the bit on position x in d j will convert d j back to d i , and we will have the DLUT S DLUT , s.t: Since the elements in the S-box S , which corresponds to the Decimal Look-up Table S DLUT , are now a permutation of all the 2 n integers in the interval [0, 2 n − 1], S is bijective. Furthermore, the permutation S DLUT is exactly one transposition away from the permutation S DLUT .
We have shown that if we start from a random bijective S-box S it is possible to construct another bijective S-box S by flipping exactly two bits from the LUT of S. It appears that the first flip can be on a random element in the LUT, but the second flip is uniquely determined by the first one. We define each such pair as coupling.
Definition 17 (Couplings): Let us take a bijective S-box S(n, n) and its corresponding DLUT We define as a coupling each set

Corollary 1 (Properties of Couplings Pivot Sets):
Considering the definitions of the Couplings Pivot Sets on bijective S-box S(n,n), the following properties hold true: Definition 19 (Coordinate Decomposition): Let S be an (n, n) bijective S-box. We take a random element with coordinates (x, y) of its corresponding linear approximation table S LAT . We denote the binary representation of y as: The coordinate decomposition of an element with coordinates (x, y), denoted by x,y , is the set:

Definition 20 (Nonlinearity Bottleneck Snapshot -NBS):
We define the nonlinearity bottleneck snapshot S NBS of a bijective S-box S(n, n) as a set of tuples holding all coordinates of the elements of S LAT , which are holding down the nonlinerity value S NL of S, i.e.
For a given bijective S-box S(n, n), we say that coordinate j is descendible, if the following properties hold:

S c i = S • c i
When we have a list of couplings {c 1 , c 2 , · · · , c i }, which we want to use for transformation of S in this exact order, we will use the following expression: In the case when c a = c b , since they belong to the same coupling pivot set, it follows that c a ∩ c b = ∅, which concludes the proof.
Corollary 3: For a given bijective S-box S(n, n), for any couplings c j , which belongs to the same couplings pivot set Proof: For a given bijective S-box S(n, n) and some coupling c = d x , f (d x ) j , we denote as S c LAT the transformed LAT of S caused by c. Let us take some element e x,y from the LAT of S before the coupling transformation. We have: , for some linear function L q and some linear combination in binary representation of the coordinates of S: b = b 1 b 2 · · · b 2 n . If j / ∈ x,y , e x,y is not affected after applying the coupling. However, if j ∈ x,y , we know that exactly two of the bits of the linear combination b are flipped. We denote them as b s and b t . Let us denote the element on position (x, y) on the newly created LAT as e x,y .
Since the expression ±1 ± 1 is equal to one of the three possible values: -2, 0, and 2, the proof concludes.

Definition 26 (S-box Coordinate Extended LAT -CELAT):
For a given bijective S-box S(n, n), and a given coordinate i, we can define the one-dimensional linear approximation table of S as: Furthermore, we denote all the couplings in the couplings pivot set {S i } as c 1 , c 2 , · · · , c 2 n−1 . We have: Following the same concept used in the construction of one-dimensional LAT of S, we can define one-dimensional CTM, i.e.: Finally, we define S-box i-th Coordinate Extended LAT S i CELAT as the following table: CELAT has 2 n−1 + 1 rows and 2 2n columns. For example, let us consider an S-box S(2, 2) with S DLUT = [0, 2, 1, 3]. For n = 2 the S 2 CELAT has 2 2 − 1 = 3 rows and 2 2 * 2 = 16 columns. Considering coordinate 2, we have:

IV. TRANSLATING AN S-BOX NL OPTIMIZATION PROBLEM TO A BINARY INTEGER PROGRAMMING MODEL
In this chapter, the transition from the S-box optimization problem to a system of linear constraints is presented. We are following the same notations as in [23]. where the data consists of the row vector c with size n, (m, n) matrix A, and column vectors b and x with respective sizes of m and n. The columnn vector x contains the binary variables to be optimized. We say that the set S is the set of fesible solutions, i.e.: Definition 29 Binary Integer Programming -Feasibility or SAT Problem: A feasibility binary integer program is a problem of the form: where the data consists of (m, n)-matrix A and column vectors b and x with respective sizes of m and n. The column vector x contains the binary variables to be optimized. We say that the set S is the set of feasible solutions, i.e.: In the context of the feasibility problem we are looking for just one element in the set S, not the optimal one.
For an (n, n) S-box S, we denote 2 n−1 by r and 2 2n by m. Let us construct its CELAT using coordinate i i.e: We want to apply some coupling transformations subset P = p 1 , p 2 , · · · , p k which belongs to the pivot coupling set From corollary 3 it follows that: Then, we can construct the following system of equations: q 1 = l 1 + c 11 x 1 + c 21 x 2 + · · · + c r1 x r q 2 = l 2 + c 12 x 1 + c 22 x 2 + · · · + c r2 x r · · · q m = l m + c 1m x 1 + c 2m x 2 + · · · + c rm x r (4) where x = (x 1 , x 2 , · · · , x r ) ∈ B r , and x t = 1 iff p t ∈ P. We have S NL = 2 n−1 − max m j=1 abs(l j ). If coordinate i is VOLUME 9, 2021 descendable, we can construct the following binary integer programming feasibility problem: where A is a column vector with 2 n−1 +1 elements, each equal to 2 n−1 − S NL − 2, while B is a column vector with 2 n−1 + 1 elements, each equal to S NL −2 n−1 +2. Let us denote the SAT problem descending on coordinate i in equation 5 as S,i . This is NP-hard problem with a total of 2 n−1 binary variables and 2 n + 2 restrictions. However, we can further divide the problem to an union of subproblems, i.e.: where each subproblem d S,i is modelled using the following restrictions: Solving any of the subproblems will yield a solution to the original problem.
are always satisfied. Proof: For a subproblem d S,i we have the following restrictions: l 1 + c 11 x 1 + c 21 x 2 + · · · + c r1 x r ≤ 2 n−1 − S NL − 2 l 1 + c 11 x 1 + c 21 x 2 + · · · + c r1 x r ≥ S NL − 2 n−1 + 2 l 2 + c 12 x 1 + c 22 x 2 + · · · + c r2 x r ≤ 2 n−1 − S NL − 2 l 2 + c 12 x 1 + c 22 x 2 + · · · + c r2 x r ≥ S NL − 2 n−1 + 2 From lemma 6, we know that the possible values of the elements c ij are −2, 0 or 2. Hence: If for some l j the following inequalities hold: then and on the other hand l 1 + c 11 x 1 + c 21 x 2 + · · · + c r1 x r ≥ l 1 + min(c 11 x 1 + c 21 x 2 + · · · + c r1 x r ) which completes the proof. Definition 30 (CELAT with radius R): For a given bijective S-box S(n, n), and a given coordinate i, we have: We define as S i,R CELAT a matrix constructed of those columns of S i CELAT with first element ρ, for which the following inequalities hold: Hence, a given suproblem d S,i could be further reduced and launched on S i,d CELAT , instead of its corresponding full (unreduced) version S i CELAT .

V. ALGORITHM
By using automata notation, Figure 1 presents the distinct steps of the optimization process. State S is the initial state of the automata. In this phase, we initialize and process the input. We make some additional checks about the properties of the S-box. For example, we check the bijectivity property of the S-box. We further analyze and extract the descendable coordinates (if such exist). For example, if we first choose the coordinate j to descent into, we further calculate the corresponding matrix S j CELAT . Then, the adjacent state Re, by further processing the generated matrix, and by using some radius R, generates the matrix S j,R CELAT . Finally, state BIP is translating the problem to a binary integer programming feasibility problem. In case a feasible solution or proof that the problem is infeasible is found, the result is reported back to state Re, the major role of which is to orchestrate the behavior of the optimization routine -increasing the radius, changing the descendible coordinate, or give up.
Few important properties of the automata should be emphasized. For a given S-box S, if S = 1, then at least one descendible coordinate, for example j, does exist. Thus, if a feasible solution of R S,j is found, the nonlinearity of S could be increased by activating exactly R couplings. Therefore, we could not only optimize the nonlinearity of S, but dictate the impact of our changes to the original S-box as wellincreasing the value of R will increase the total count of flipped bits in S.
In those cases where S > 1, the aforementioned algorithm, as we will later demonstrate, is still applicable. We just pick a random coordinate j ∈ S , instead of a descendible one. As a consequence, finding a solution to the R S,j problem will not increase the nonlinearity of S. However, it will decrease the value of S by 1. Thus, by repeating the reduction phase, we would eventually reduce the initial problem to a problem having S = 1, for some S-box S , yielded by the optimization routine, or the composition of optimization routines, performed on S.

VI. EXPERIMENTAL RESULTS
We have implemented the algorithm by using Python, for states S and Re, and the Gurobi SAT Solver [24], for the BIP state itself. We analyzed two famous S-boxes: the Skipjack S-box, developed by the U.S. National Security Agency (NSA) [25], which we will denote as Sk, and the Kuznyechik S-box, standardized by the Russian Federation's standardization agency [26], which we will denote as Kk. On the other hand, if we require a higher nonlinearity, combined with a greater impact of the structure of Sk, we could significantly increase the value of R. Indeed, using a radius value of 10, the translated BIP model of 10 Sk,4 consists of 25125 rows, 128 columns, and 1621980 nonzeros. Despite the greater model, after 14.542 seconds, a solution was found, yielding an S-box with nonlinearity 102, constructed from Sk by flipping exactly 20 bits.  (1,3,5,7), (1,2,4,7), (0,2,5,6,7), (0,3,5,6,7), (3,4,5), (0,2,4,6,7), (1,2,6), (1,2,4,6), (2,4,5,6,7), (4,5,6,7)} and Kk = 2. Since the degree of descendibility is greater than 1, more precisely Kk = {2, 5}, we pick a random coordinate from Kk . 1 Kk,5 and 2 Kk,5 are reported back as infeasible. However, 3 Kk,5 , consisting of 319 rows, 128 columns and 20904 nonzeros, is feasible. Again, the solver took less than a second to find a solution, i.e.  (1,2,4,6), (1,2,4,7)} and Kk = 1. As expected, the value of is decreased by 1. Thus, we continue the optimization process by the descendible coordinate with index 2. First, the infeasibility of the two models 1 Kk,2 and 2 Kk,2 are proved. However, the BIP model of 3 Kk,2 , consisting of 379 rows, 128 columns, and 25344 nonzeros, yielded a solution for less than a second. Indeed, the found coupling set {(95, 127), (207, 239), (108, 157)} increased the nonlinearity of Kk from 100 to 102. Hence, we have shown how Kk could be optimized to an S-box with higher nonlinearity by tweaking just 12 bits. VOLUME 9, 2021 C. ACNV PROBLEM The ACNV optimization problem could be represented as a special, and significantly lighter, in terms of computational burden, case of S i,R CELAT , where S denotes the initial S-box and i denotes the coordinate of S to be optimized. Since our goal is ACNV optimization only, we could significantly reduce the size of S i,R CELAT . Let us denote as S i,R ACNV , the matrix formed by the matrix S i,R CELAT , with columns corresponding to linear combinations of coordinates of S removed. Indeed, this is a significant reduction of the feasibility model. For example, if S is of dimension (n, n), S i,R CELAT has 2 n−1 + 1 rows and 2 2n columns, while S i,R ACNV has 2 n−1 + 1 rows and 2 n columns. We further denote the corresponding feasibility problem corresponding to S i,R ACNV as S,i . As usual, we could divide the problem S,i to a union of subproblems d S,i . We want to emphasize, that collecting all published results regarding the ACNV problem is a rigid and tedious task. However, during our extensive research, in Table 1, we have collected the most significant of them. It appears, that the majority of papers focused on optimizing the ACNV value of (8,8) S-boxes, were skeptical about the existence of (8,8) S-box with ACNV greater than 112. Indeed, as stated in [70], ''It is seen that the average value obtained from the proposed method is close to 112 (theoretical optimal value for an 8 × 8 S-box)''. However, as we have already shown in [18], bijective (8,8) S-boxes with ACNV up to 114 are easy to construct.
We have initiated the optimization routine on a bijective S-box, for simplicity denoted as S, from [18], possessing the highest, currently known, ACNV of 114.5. It is composed of 6 coordinates with a nonlinearity value of 114 and 2 coordinates with a nonlinearity value of 116. We will outline a possible trace of improvement, which led to an S-box with an ACNV of 116.0.
• We launched ≤9 S,4 . After around 153 seconds a feasible solution, with exactly 9 couplings, was found. We activated the couplings to get S 1 . ACNV was lifted to 114.75.
• We launched ≤8 S 1 ,1 . After around 16 seconds a feasible solution, with exactly 8 couplings, was found. We activated the couplings to get S 2 . ACNV was lifted to 115.
• Finally, we launched ≤9 S 5 ,5 and a feasible solution was found after 69 seconds. We activated the 9 couplings to get S 6 . ACNV was lifted to 116. We present S 6 as a hexadecimal format: a9, 7b, 99, 49 The overall nonlinearity of S 6 is 92. Significant efforts were made to reach higher ACNV -reaching higher ACNV would reveal a balanced Boolean function having a nonlinearity of 118. Unfortunately, all instances x S 6 ,y , for ∀x, y : x ≤ 21, y ∈ [1,8], were proofed infeasible. We want to emphasize that the searching routine is highly efficient. For example, 17 S 6 ,y , for some y, was proved infeasible for 1065 seconds or approximately 18 minutes. Since the given search space size is equal to 128 17 , or approximately 2 69 , this results in checking simultaneously roughly 2 59 distinct cases per second.

VII. CONCLUSION
In this work, we have projected the S-box nonlinearity optimization problem to the domain of binary integer programming. The provided framework is not only capable of increasing the nonlinearity value of a given S-box, but also achieving that with a minimum impact on the initial structure. We have demonstrated how the Skipjack S-box and the Kuznyechik S-box could be both optimized to a higher nonlinearity by tweaking, respectively, just 4 and 12 bits. We further revealed that bijective (8,8) S-boxes, the eight coordinates of which possess the currently known optimal nonlinearity value of 116, do exist.