Enumerating Discrete Resonant Rossby/Drift Wave Triads and Their Application in Information Security

Abstract

We propose a new parametrization of the resonant Rossby/drift wave triads to develop an algorithm to enumerate all resonant triads in a given grid of wavenumbers. To arrive at such a parametrization, we have employed tools from arithmetic/algebraic geometry to project resonant triads on a certain class of conics. Further, we extend the newly developed algorithm for the enumeration of quasi-resonant triads and experimentally show that the said algorithm is robust to design the network of quasi-resonances. From the experimental results, we observed that the new algorithm enumerates all triads in low computation time when compared with the existing methods. Finally, we apply this work to information security by constructing a total order on the enumerated resonant triads to design a substitution box (S-box) generator. Via extensive analyses over several indicators (nonlinearity, algebraic complexity, linear and differential approximation probabilities, strict avalanche criteria, and bit independence criterion) we show that the newly developed S-box outperforms the S-boxes constructed by most of the existing schemes.

1. Introduction

This paper combines our recent results from two seemingly unrelated scientific areas: geophysical waves on the one hand [1,2,3,4], and information security on the other hand [5,6,7,8], to produce an interdisciplinary application, whereby a calculation of exact discrete wave resonances in a geophysical system leads to the solution of a number-theory problem, which in turn we use to develop an algorithm to construct a cryptographically robust substitution box (S-box) generator.

Below we present brief reviews of the geophysical aspects (Section 1.1) and of the information security aspects (Section 1.2), then we summarize the contributions of this paper (Section 1.3), and finally we provide a plan for the remaining sections.

1.1. A Brief Review of Rossby/Drift Wave Triads

Nonlinear wave resonances play an important role in a variety of systems, ranging from fusion reactors to weather prediction to water waves. The review by Horton and Hasegawa [9] presents, with historical and scientific perspectives, the analogy between Rossby waves (pertaining to atmospheric and oceanic physics) and drift waves (pertaining to plasma physics). For simplicity, we will restrict our discussion to the atmospheric side. Resonant and quasi-resonant Rossby/drift wave triads play a vital role in atmospheric dynamics [10,11]. Petoukhov et al. [12] linked Rossby waves with several extreme weather phenomena via quasi-resonances. Coumou et al. [13] explored a connection between Rossby waves and global warming. The enumeration of all resonant Rossby wave triads is a practical problem and recently has gained great importance. In [14,15], generic algorithms were used to enumerate resonant Rossby wave triads in a given box. Bustamante and Hayat used elliptic surfaces to classify resonant and quasi-resonant triads [1]. Kopp [16] used methods from projective geometry to obtain almost all resonant triads in a box of wavenumber size 5000.

A Rossby wave [12,13,17,18,19,20,21,22,23,24] will be defined in this paper as a parametrized solution of the linearized form of the Charney–Hasegawa–Mima equation (CHME) [1], a partial differential equation that expresses the conservation of potential vorticity [25]. In other words, a Rossby wave represents a plane-wave solution of the linearized version of CHME. There are many such solutions to the linearized CHME, and these are important for the study of resonances. Mathematically, in the context of the CHME with periodic boundary conditions, a Rossby/drift wave is determined by a wavevector $\mathbf{k}\in {\mathbb{Z}}^2$. The first component of this vector represents the number of peaks of the wave along the zonal direction, and the second component represents the number of peaks along the meridional direction. In a special case, the nonlinear interaction of two waves of different wavevectors enumerates a third wave with a new wavevector, thus forming a ‘triad’ of wavevectors such that the interaction of any two waves in the triad produces the third one. In such a case, the nonlinear interaction is limited to only three waves under the constraint that no further waves may be produced. Such a triple of waves is known as a resonant triad.

The equations determining a resonant triad in CHME are Diophantine equations, namely, equations for integers. Solutions to these equations are of importance in reduced models of atmosphere and plasmas, as they give the wavevectors involved in three-wave interactions. The enumeration of all resonant triads is a practical problem. Over several centuries, classical methods were developed by Fermat, Euler, Lagrange, and Minkowski to classify solutions to some Diophantine equations [26], but the CHME resonant triads lead to new Diophantine equations whose analytical and computational enumeration problem has received great attention. Bustamante and Hayat [1] proposed a new method to enumerate numerically all resonant triads. The newly developed method relies on the transformation of the wavevectors to a new set of variables, converting the Diophantine equations for CHME resonant triads into a simpler set of equations, solvable by Fermat’s Xmas theorem. They extended the algorithm to include the enumeration of quasi-resonant triads, and the extended method was found practical. Kopp [16] provided a new parametrization to the resonant triads and found almost all the resonant triads in a box of wavenumber size 5000 using his parametrization. Subsequently, Hayat et al. [4] explicitly parametrized the resonant wavevectors by two rational parameters and proposed a new method of cubic complexity to enumerate all irreducible triads in a specific region.

1.2. A Brief Review of Substitution Box

Transfer of useful information via internet is usual in today’s modern life. In cryptography, a substitution box (or S-box) is a nonlinear vectorial Boolean function with m input and n output bits. Here, m and n are two positive integers. For data encryption, mainly two techniques are used: block cipher and stream cipher [27]. A block cipher encrypts data in blocks of fixed length, whereas in a stream cipher one bit is encrypted in one go. In a block cipher, an S-box is used as a fundamental nonlinear part to achieve a higher level of confusion in the data [28]. Furthermore, the security of the S-box dependent cryptosystems mainly relies on the strength of the S-box.

Using various mathematical methods, many S-box generators have been developed in recent years [7,29,30] for possible applications in image encryption algorithms using S-boxes as nonlinear components [31,32]. Cryptographers are widely using chaotic maps to develop new algorithms for S-box generation [32,33,34,35,36,37,38,39]. Optimization techniques are adopted to design highly nonlinear S-boxes in [40,41,42]. Other domains of mathematics, such as graph theory [43], cubic polynomial mappings [44], and linear trigonometric transformations [45], are also used for the construction of S-boxes. In [46], a new S-box construction algorithm is developed using chaotic maps, symmetric groups, and Mobius transformations to generate a highly secure S-box.

Elliptic curves are another important structure to construct secure S-boxes [6,47,48]. In [5], Azam et al. proposed an efficient S-box generator over elliptic curves to construct dynamic S-boxes. Saleh and Abbas [49] designed an S-box generator using the points over an elliptic curve to construct highly secure and key-dependent S-boxes. Ullah et al. [50] designed efficient S-boxes and pseudo-random number generators over elliptic curves. Hayat et al. [7] presented an S-box generation technique using an elliptic curve over a finite ring of integers. The designed S-boxes are further evaluated for the encryption of images. Recently, Murtaza et al. [51] introduced a dynamic S-box generator using elliptic curves and binary sequences to design robust and dynamic S-boxes for block ciphers.

1.3. Our Contributions

Motivated by the above work, we propose a new parametrization of the resonant triads to develop a new algorithm for the enumeration of all resonant triads in a specific box. The proposed work is capable of enumerating all triads with cubic time complexity. In addition, a new total order is defined using the proposed parametrization. The newly developed order and parametrization are employed to design a novel substitution box (S-box) generator.

Our main results are the following:

(1)

A new geometric parametrization with low time complexity to enumerate all resonant triads, our parametrization takes only 2 days to enumerate 472 irreducible triads in a grid of size 5000, whereas parametrization in [4] took 10.5 days on the same machine;

(2)

The new parametrization to triads enumerates 472 irreducible resonant triads in a grid of size 5000 when compared with the algorithm in [16], which enumerates 463 irreducible resonant triads;

(3)

An extended method to enumerate quasi-resonant triads as well;

(4)

As an application, we define a total order on the set of triads to develop an efficient S-box generator;

(5)

Our analysis shows that our generator is capable of constructing a cryptographically strong S-box.

The remaining parts of this paper are organized as follows: Section 2 explains the basic concepts of CHME. Section 3 discusses the new parametrization for the enumeration of triads. Section 4 contains the proposed ordering, S-box generator, and their detailed analyses. Section 5 presents the conclusion.

2. The Charney–Hasegawa–Mima Equation: Periodic Boundary Conditions and Arbitrary Aspect Ratio f

The discussion in this section is based on [1] and references therein. We consider the partial differential equation known as the barotropic vorticity equation in the $\beta$-plane approximation, also known as the Charney–Hasegawa–Mima equation (CHME):

$$\frac{\partial }{\partial t}(▵\psi -F\psi )+\beta \frac{\partial \psi }{\partial x_1}+[\psi ,▵\psi ]=0,$$

(1)

where, in the atmospheric context, $\psi :{\mathbb{R}}^2\times [0,T)\to \mathbb{R}$ is the streamfunction, $\beta >0,F\ge 0$ are constants, and

$$▵\psi =\frac{{\partial }^2\psi }{\partial x_1^2}+\frac{{\partial }^2\psi }{\partial x_2^2},[A,B]=\frac{\partial A}{\partial x_1}\frac{\partial B}{\partial x_2}-\frac{\partial A}{\partial x_2}\frac{\partial B}{\partial x_1}.$$

The sum of the first and second term represent the linear part of (1); the last term $[\psi ,▵\psi ]$ is the nonlinear part of (1). From here on, we consider the case of periodic boundary conditions $(x_1,x_2)\in [0,2\pi )\times [0,2\pi f^{-1})$, where $f>0$ is the so-called aspect ratio. To begin with, we consider the case of unit aspect ratio, $f=1$. The general solution of (1) is not known, as it generically displays spatio-temporal chaos and turbulence. However, a number of particular exact solutions, known as Rossby waves, are available in the form of a family of cosine functions

$$\psi ={\psi }_{k,\ell }(x_1,x_2,t)=cos\left(k{x}_1+\ell x_2-\omega t\right),$$

(2)

for the angular frequency

$$\omega =\omega \left(k,\ell \right)=\frac{-\beta k}{k^2+{\ell }^2+F}.$$

(3)

Notice that, for these solutions, both the linear part and the nonlinear part of (1) vanish independently. The vector $(k,\ell )$ is a wavevector, whereas the integers k and ℓ are known as the zonal and meridional wavenumbers, respectively.

When considering interacting Rossby waves with different wavevectors, nonlinear terms do not vanish. In fact, the dynamics leads to triad interactions, and in weakly nonlinear regimes an important role is played by the so-called resonant triads [52]. A resonant triad is a triple $(k_1,{\ell }_1),(k_2,{\ell }_2),(k_3,{\ell }_3)$ of Rossby waves satisfying the set of equations

$$\begin{array}{cc}\hfill k_2=& k_3-k_1,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\ell }_2=& {\ell }_3-{\ell }_1,\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill {\omega }_2=& {\omega }_3-{\omega }_1,\hfill \end{array}$$

(6)

for ${\omega }_i=\omega (k_i,{\ell }_i),i=1,2,3$. If for a small positive number $\delta$, condition (6) on angular frequencies in the above set of equations is replaced by

$$|{\omega }_1+{\omega }_2-{\omega }_3|\le \delta ,$$

(7)

then the aforesaid triple becomes a quasi-resonant triad, and the positive number $\delta$ is called a detuning level.

In the special case where $F=0$, the parameterized solutions for the resonant triads were first obtained in [1] and subsequently developed in [16], giving

$$\begin{array}{cc}\hfill \frac{k_1}{k_3}=& (u^2+t^2)(u^2+t^2-2u)/(1-2u),\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \frac{{\ell }_3}{k_3}=& \left((2u-1)u+(u^2+t^2)(u^2+t^2-2u)\right)/\left((1-2u)t\right),\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \frac{{\ell }_1}{k_3}=& (u^2+t^2)\left((2u-1)+(u^2+t^2-2u)u\right)/\left((1-2u)t\right).\hfill \end{array}$$

(10)

An independent and particular case is discussed in [53] as well. In [1], the resonant triad is explicitly transformed to a triple $(x,y,d)$ with rational components as

$$\frac{k_1}{k_3}=\frac{x}{y^2+d^2},\frac{{\ell }_1}{k_3}=\left(\frac{x}{y}\right)\left(1-\frac{d}{y^2+d^2}\right),\frac{{\ell }_3}{k_3}=\frac{d-1}{y}.$$

(11)

Then $(x,y,d)\in {\mathbb{Q}}^3$ is inversely mapped to a triad via

$$x=\frac{k_3(k_1^2+{\ell }_1^2)}{k_1(k_3^2+{\ell }_3^2)},y=\frac{k_3(k_3{\ell }_1-k_1{\ell }_3)}{k_1(k_3^2+{\ell }_3^2)},d=\frac{k_3(k_3{k}_1+{\ell }_1{\ell }_3)}{k_1(k_3^2+{\ell }_3^2)}.$$

(12)

Let us consider now the case of general aspect ratio f. It is enough to consider aspect ratios whose squares are rational, so we can write $f=f_1\sqrt{f_2}$ for rational $f_1$ and square-free integer $f_2$. As shown in [1], in this case relations (1)–(7) hold true except for Equation (3), which must be replaced with

$$\omega =\omega \left(k,\ell \right)=\frac{-\beta k}{k^2+f^2{\ell }^2}.$$

(13)

For such a choice of f, the mappings analogous to Equation (11) take the form

$$\frac{k_1}{k_3}=\frac{x}{f^2{y}^2+d^2},\frac{{\ell }_1}{k_3}=\left(\frac{x}{f^{2}y}\right)\left(1-\frac{d}{f^2{y}^2+d^2}\right),\frac{{\ell }_3}{k_3}=\frac{d-1}{f^{2}y},$$

(14)

with the inverse mappings as follows:

$$x=\frac{k_3(k_1^2+f^2{\ell }_1^2)}{k_1(k_3^2+f^2{\ell }_3^2)},y=\frac{k_3(k_3{\ell }_1-k_1{\ell }_3)}{k_1(k_3^2+f^2{\ell }_3^2)},d=\frac{k_3(k_3{k}_1+f^2{\ell }_1{\ell }_3)}{k_1(k_3^2+f^2{\ell }_3^2)}.$$

(15)

3. New Parametrization of the Elliptic Surface of Resonant Discrete Rossby Wave Triads

In the notation of Section 2, the following equation of an elliptic surface was derived in [1] for resonant Rossby wave triads, namely those triads satisfying Equations (4)–(6):

$$f^2{y}^2=x^3-2d{x}^2+2dx-d^2,$$

(16)

where $x,y$, and d are defined in Equation (15), and f is the aspect ratio of the system. We fix the variable x to equal some constant $\xi$. It is direct to see that the surface (16) can be rewritten as the following equation of an ellipse:

$$\frac{f^2{y}^2}{r^2}+\frac{{(d+a)}^2}{r^2}=1,$$

(17)

where $a={\xi }^2-\xi$ and $r=\sqrt{{\xi }^4-{\xi }^3+{\xi }^2}$. Now transform (17) into a polar form by substituting

$$y=\frac{r}{f}cos\theta ,d+a=rsin\theta$$

(18)

for a parameter $\theta \in [0,2\pi )$. In summary, we have

$$x=\xi ,y=\frac{\sqrt{{\xi }^4-{\xi }^3+{\xi }^2}}{f}cos\theta ,d=\sqrt{{\xi }^4-{\xi }^3+{\xi }^2}sin\theta -{\xi }^2+\xi .$$

(19)

From (18) and (19), the inverse relation is obtained as

$$\xi =x,\theta ={tan}^{-1}\left(\frac{d+x^2-x}{fy}\right).$$

(20)

Now using the transformation (19) and the fundamental identity ${cos}^2\theta +{sin}^2\theta =1$, we get

$$f^2{y}^2+d^2=2{\xi }^4-3{\xi }^3+2{\xi }^2-2({\xi }^2-\xi )\sqrt{{\xi }^4-{\xi }^3+{\xi }^2}sin\theta .$$

(21)

Thus, using (14) and (19), we have the following three equations defining the triad wavenumbers in terms of the new parameters $\xi$ and $\theta$:

$$\begin{array}{cc}\hfill \frac{k_1}{k_3}=& \frac{1}{2{\xi }^3-3{\xi }^2+2\xi -2(\xi -1)\sqrt{{\xi }^4-{\xi }^3+{\xi }^2}sin\theta },\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill \frac{{\ell }_1}{k_3}=& \frac{\xi }{f\sqrt{{\xi }^4-{\xi }^3+{\xi }^2}cos\theta }\left(1-\frac{\sqrt{{\xi }^4-{\xi }^3+{\xi }^2}sin\theta -{\xi }^2+\xi }{2{\xi }^4-3{\xi }^3+2{\xi }^2-2({\xi }^2-\xi )\sqrt{{\xi }^4-{\xi }^3+{\xi }^2}sin\theta }\right),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \frac{{\ell }_3}{k_3}=& \frac{\sqrt{{\xi }^4-{\xi }^3+{\xi }^2}sin\theta -{\xi }^2+\xi -1}{f\sqrt{{\xi }^4-{\xi }^3+{\xi }^2}cos\theta }.\hfill \end{array}$$

(24)

Now from (15) and (20), $\xi$ and $\theta$ can be written as

$$\xi =\frac{k_3(k_1^2+f^2{\ell }_1^2)}{k_1(k_3^2+f^2{\ell }_3^2)},$$

(25)

$$\theta ={tan}^{-1}\left[\left((k_3-k_1)k_1+({\ell }_3-{\ell }_1)f^2{\ell }_1+\left(\frac{k_3{(k_1^2+f^2{\ell }_1^2)}^2}{k_1(k_3^2+f^2{\ell }_3^2)}\right)\right)\left(\frac{1}{f(k_3{\ell }_1-k_1{\ell }_3)}\right)\right].$$

(26)

Hence, (22)–(24) and the inverse (25) and (26) represent the explicit parametrization of the resonant triad in terms of the parameters $\xi$ and $\theta$. Using this new parametrization, we enumerate all irreducible resonant triads in a specific box such that $\left|k\right|,|\ell |\le \mathrm{L}$. For this, choose two rational numbers $a^{\prime },b^{\prime }$ and an integer e to design a set of rational numbers ${}_{a^{\prime }}{\mathrm{A}}_{b^{\prime }}$ of size n with end points $a^{\prime },b^{\prime }$ such that $a_1=a^{\prime },a_n=b^{\prime }$, and $a_i=a_{i-1}+h$ for $h=\frac{b^{\prime }-a^{\prime }}{e}$ and $1<i<n$. Further, select a subset ${}_{c^{\prime }}{\mathrm{B}}_{d^{\prime }}\subseteq [0^{\circ },{360}^{\circ }]$ such that $d^{\prime }=c^{\prime }+tg$ for some integer t and a rational number g. The step by step technique for the enumeration of resonant triads is explained in Algorithm 1.

Algorithm 1: Enumeration of resonant triads

Input: Two sets ${}_{a^{\prime }}{\mathrm{A}}_{b^{\prime }},{}_{c^{\prime }}{\mathrm{B}}_{d^{\prime }}$, a box size $\mathrm{L}$ and an aspect ratio f. Output: A set $\mathrm{T}$ of resonant triads.       /* $\mathrm{T}$ is a set of resonant triads and initialized as an empty set. Furthermore, △ represents an arbitrary triad and $k_1^{{}^{\prime }},{\ell }_3^{{}^{\prime }},{\ell }_1^{{}^{\prime }}$ are the right hand sides of (22)–(24), respectively. Moreover, $\lfloor ·\rceil$ is an integer function.       $\mathrm{T}:=\varnothing$;       for $\xi \in {}_{a^{\prime }}{\mathrm{A}}_{b^{\prime }}$ do           for $\theta \in {}_{c^{\prime }}{\mathrm{B}}_{d^{\prime }}$ do               Compute $k_1^{{}^{\prime }},$${\ell }_3^{{}^{\prime }}$ and ${\ell }_1^{{}^{\prime }}$ for $\xi ,\theta$ by using (22)–(24).               for $k_3\in [1,\mathrm{L}]$ do                  $k_1:=\lfloor k_1^{\prime }·k_3\rceil ,{\ell }_3:=\lfloor {\ell }_3^{\prime }·k_3\rceil$ and ${\ell }_1:=\lfloor {\ell }_1^{\prime }·k_3\rceil$;                  $k_2:=k_3-k_1,{\ell }_2:={\ell }_3-{\ell }_1$ and ${\omega }_i:=k_i/(k_i^2+f^2{\ell }_i^2),i:=1,2,3;$                  if ${\omega }_3:={\omega }_2+{\omega }_1$ and $|k_i|,|{\ell }_i|<\mathrm{L},i:=1,2,3;$ then                      $\mathrm{T}:=\mathrm{T}\cup \{△\}$;           Output $\mathrm{T}$ as a set of triads.

We borrow an idea from physical considerations: to fully understand the dynamics of a system of Rossby waves, it is necessary to understand the behavior of quasi-resonant triads. Therefore, to further investigate the newly designed parametrization, we enumerate quasi-resonant triads using (4)–(6), where (6) is replaced by (7). In practice, we apply the same parametrization (22)–(24) but we approximate the output wavenumbers to numbers within a box of size L. This leads to a mismatch in the frequency resonance condition. The value of the newly introduced “detuning parameter” $\delta$ from Equation (7) determines the number of quasi-resonant triads obtained. Our proposed parametrization directly computes the quasi-resonant triads. For this purpose, we select the parameters $f=1,a^{\prime }=1.025,b^{\prime }=1.204,c^{\prime }=0,d^{\prime }=360,e=716,g=0.0125,\mathrm{L}=100$ and several choices of $\delta$: $2\times {10}^{-7},4\times {10}^{-7},8\times {10}^{-7},2\times {10}^{-6},4\times {10}^{-6},8\times {10}^{-6},2\times {10}^{-5}$. The wavevectors of the quasi-resonant triads enumerated for the chosen parameters are illustrated in Figure 1.

A bar graph analysis for all the computed unique and irreducible quasi-resonant triads is given in Figure 2.

It is evident from Figure 1 and Figure 2 that the proposed method is capable of enumerating a large number of triads by relaxing the condition on the angular frequencies. In Example 1, we explain how the designed algorithm maps a triad on the surface point and vice versa.

Example 1.

To map the triad on the surface and hence on the conic, let us choose the triad $(1,-11)$, $(8,34)$, $(9,23)$, then calculate $\xi =\frac{9}{5}$ and augmented angle $\theta ={180}^{\circ }+{\frac{1807}{36}}^{\circ }$ by (25) and (26), respectively. For $\xi =\frac{9}{5}$, we have $a=\frac{36}{25}$ and $r=\frac{\sqrt{4941}}{25}$. Hence from (19), it follows that $x=\frac{9}{5},y=-\frac{9}{5}$, and $d=-\frac{18}{5}$. Now to map the surface point back to triad, take the point $(x,y,d)=(\frac{9}{5},-\frac{9}{5},-\frac{18}{5})$ on the elliptic surface, and by (20) we have $\xi =\frac{9}{5}$ and augmented angle $\theta ={180}^{\circ }+{\frac{1807}{36}}^{\circ }$, so that (11) gives that $\frac{k_1}{k_3}=\frac{1}{9},\frac{{\ell }_1}{k_3}=-\frac{11}{9},\frac{{\ell }_3}{k_3}=-\frac{23}{9}$. From this, we can write that $k_3=9,k_1=1,{\ell }_1=-11$, and ${\ell }_3=23$. Using these values, we compute $k_2=k_3-k_1=8,{\ell }_2={\ell }_3-{\ell }_1=34$. Hence, we enumerate the same triad $(1,-11)$, $(8,34)$, $(9,23)$.

It was verified in [4] using a brute-force method that there is a total of 472 irreducible resonant triads in a grid of size 5000. Previously, Kopp [16] introduced a fast parametrization to enumerate irreducible triads. Using this parametrization and a simple search, he could enumerate only 463 out of the 472 irreducible resonant triads in a grid of size 5000. For this reason, the authors in [4] developed another parametrization. This parametrization could enumerate all 472 irreducible resonant triads in a grid of size 5000, but it took 10.5 days to compute this on a 16-core-machine. To overcome the shortcomings of the aforementioned algorithms, in this paper we introduced parametrization (22)–(24) to irreducible triads. A major advantage of our parametrization over the existing ones is that we can enumerate all 472 irreducible resonant triads in a given grid of size 5000 in only 2 days on a 16-core-machine.

Table 1. Number of irreducible triads for aspect ratio $f=1$.

Methods	Time	No. of Triads		System Specifications
		Optimal	Obtained
Proposed	2 days	472	472	MATLAB, 16-core-machine
Ref. [4]	10.5 days	472	472	Mathematica, 16-core-machine
Ref. [16]	-	472	463	Mathematica, MacBook 2.9 GHz dual-core i7

indicates that, for aspect ratio $f=1$, the new method to enumerate triads is more efficient than the existing algorithms.

For different values of aspect ratio f including the standard case ($f=1$), we have computed the number of irreducible resonant triads in a grid of size 100 using a brute-force method. The results are shown in

Table 2. Number of irreducible triads against different values of aspect ratio f.

f	1	${\displaystyle \frac{1}{2}}$	$\sqrt{2}$	${\displaystyle \frac{\sqrt{2}}{3}}$	$\sqrt{3}$	${\displaystyle \frac{\sqrt{3}}{2}}$	${\displaystyle \frac{\sqrt{3}}{3}}$	${\displaystyle \frac{\sqrt{3}}{4}}$	${\displaystyle \frac{\sqrt{3}}{5}}$	$2\sqrt{3}$	$3\sqrt{3}$	$4\sqrt{3}$	${\displaystyle \frac{\sqrt{6}}{2}}$	${\displaystyle \frac{\sqrt{5}}{2}}$	${\displaystyle \frac{\sqrt{7}}{2}}$
Triads	12	7	6	3	337	231	281	65	63	97	64	30	4	3	7

. Notice that the case $f=\sqrt{3}$ gives an extremely large number of resonant triads.

The plot in Figure 3 contains the largest clusters that occur within a box of size L = 100 and aspect ratio $f=\sqrt{3}$. We have excluded the repeated “mirrored” clusters. Also, we have excluded special triads with $k_3=0,|k_1|=|k_2|,$ because they are a kind of degenerate case from the viewpoint of the dynamical system. Remarkably, the Figure 3 shows resonant triads connected by two common modes, a feature that enhances the turbulent behaviour of the dynamical system [3]. The isolated triads provide roughly 50% of the overall number of resonant triads, and this fact holds true for increasing box sizes. However, within a specific box size, the distribution of cluster sizes is nontrivial.

4. An Application in Cryptography

To get the desired security in an S-box based cryptosystem, an S-box should be capable of creating enough confusion and diffusion. Total orders play an important role in achieving the aforementioned purpose. For example, Azam et al. [5] defined three different orderings on the points of elliptic curves to construct secure S- boxes. Similarly, quasi-resonant triads are ordered in [8] to develop an image encryption scheme. Motivated by [8], we develop a new ordering ${\prec }_{\theta }$ on resonant triads with respect to the new parametrization. Let △ and ${△}^{\prime }$ be two arbitrary triads (e.g., $△=(k_i,{\ell }_i)$ and ${△}^{\prime }=(k_i^{\prime },{\ell }_i^{\prime })$ for $i=1,2,3$), then the ordering ${\prec }_{\theta }$ is defined by

$$△{\prec }_{\theta }{△}^{\prime }\iff \left\{\begin{array}{c}\mathrm{either}\xi <{\xi }^{\prime },\mathrm{or}\hfill \\ \xi ={\xi }^{\prime }\mathrm{and}\theta \le {\theta }^{\prime }.\hfill \end{array}\right.$$

We prove that ${\prec }_{\theta }$ is a total order.

Lemma 1.

If $\mathrm{T}$ represents the set of irreducible resonant triads, then ${\prec }_{\theta }$ is a total order on $\mathrm{T}$.

Proof. 

We need to show that ${\prec }_{\theta }$ is reflexive, antisymmetric, and transitive. The reflexive property follows from the fact $\xi =\xi$ and $\theta =\theta$ always.

Now, suppose that $△{\prec }_{\theta }{△}^{\prime }$ and ${△}^{\prime }{\prec }_{\theta }△$ then $\xi <{\xi }^{\prime }$ and ${\xi }^{\prime }<\xi$. Therefore, it can be concluded that $\theta \le {\theta }^{\prime }$ and ${\theta }^{\prime }\le \theta$. Hence $\theta ={\theta }^{\prime }.$ Consequently, from (22)–(24), it follows that

$$\frac{k_1}{k_3}=\frac{k_1^{\prime }}{k_3^{\prime }},\frac{{\ell }_3}{k_3}=\frac{{\ell }_3^{\prime }}{k_3^{\prime }},\frac{{\ell }_1}{k_3}=\frac{{\ell }_1^{\prime }}{k_3^{\prime }},$$

(27)

which is possible if $k_1=c_1{k}_1^{\prime },k_3=c_1{k}_3^{\prime },{\ell }_3=c_2{\ell }_3^{\prime },k_3=c_2{k}_3^{\prime }$, and ${\ell }_1=c_3{\ell }_1^{\prime },k_3=c_3{k}_3^{\prime }$ for some integers $c_1,c_2$, and $c_3$, but $k_3=c_1{k}_3^{\prime }=c_2{k}_3^{\prime }=c_3{k}_3^{\prime }$ implies that $c_1=c_2=c_3=c$ for some integer c. Moreover, ${\ell }_1=c{\ell }_1^{\prime }$ and ${\ell }_3=c{\ell }_3^{\prime }$. From (4) and (5), it is evident that $k_2=k_3-k_1$ and ${\ell }_2={\ell }_3-{\ell }_1$, and it follows that $k_2=c{k}_2^{\prime }$ and ${\ell }_2=c{\ell }_2^{\prime }$. That is, we have $(k_i,{\ell }_i)=(c{k}_i^{\prime },c{\ell }_i^{\prime })$ for $i=1,2,3.$ Thus $△=c{△}^{\prime }$, but all triads in $\mathrm{T}$ are irreducible. Thus $c=1$, and ${\prec }_{\theta }$ is antisymmetric.

For transitivity, let $△{\prec }_{\theta }{△}^{\prime }$ and ${△}^{\prime }{\prec }_{\theta }{△}^{″}$. Then one possibility is $\xi <{\xi }^{\prime }$ and ${\xi }^{\prime }\le {\xi }^{″}$, which implies that $\xi <{\xi }^{″}$. The second possibility is $\xi ={\xi }^{\prime }$ and ${\xi }^{\prime }<{\xi }^{″}$, which also implies that $\xi <{\xi }^{″}$. The last case is $\xi ={\xi }^{\prime }={\xi }^{″}$, which gives that $\theta \le {\theta }^{\prime }\le {\theta }^{″}$ and hence $\theta \le {\theta }^{″}$. Thus, in all possible cases, $△{\prec }_{\theta }{△}^{″}$. Consequently, ${\prec }_{\theta }$ is transitive and hence a total order. □

Based on the ordering ${\prec }_{\theta }$, an S-box generator is introduced in the following section.

4.1. The Proposed S-box Generator

Suppose we want to construct an S-box over the set $[0,2^m-1]$ for some positive integer m. Consider the following steps:

Step 1.

Choose $\mathrm{L}$ as a grid size and generate two sets ${}_{a^{\prime }}{\mathrm{A}}_{b^{\prime }}$ and ${}_{c^{\prime }}{\mathrm{B}}_{d^{\prime }}$ as required by Algorithm 1 with the constraint that the parameters $a^{\prime },b^{\prime },c^{\prime },d^{\prime },e,g$, and t are chosen in such a way that we can enumerate exactly $u=\lceil \frac{2^m}{6}\rceil$ triads.

Step 2.

After enumerating all u triads, take the absolute of all wavevectors. Then, arrange the triads using the ordering ${\prec }_{\theta }$ to obtain a matrix ${\mathrm{M}}_{u\times 6}$, where the ith row represents the ith triad.

Step 3.

Select an integer $\ell \le |{\mathsf{\Lambda }}_{\mathrm{M}}|$ to apply the modulo ℓ operator on the matrix $\mathrm{M}$ in order to obtain the matrix ${\mathrm{M}}_{\ell }$. Here, ${\mathsf{\Lambda }}_{\mathrm{M}}$ denotes the greatest value in $\mathrm{M}$.

Step 4.

Neglect ${\mathrm{M}}_{\ell }\left(i\right)$ for $i>2^m$, and define a mapping $\varphi :[0,2^m-1]\to [0,2^m-1]$ such that $\varphi \left(n\right)=i-1$, where i represents the index of the nth least value of ${\mathrm{M}}_{\ell }$ in linear ordering.

It is noted that if ${\mathrm{M}}_{\ell }\left(i\right)$ is the rth least value and ${\mathrm{M}}_{\ell }\left(i\right)={\mathrm{M}}_{\ell }\left(j\right)$ for $i<j$ then ${\mathrm{M}}_{\ell }\left(j\right)$ is considered to be the $(r+1)$th least value of ${\mathrm{M}}_{\ell }$. For parameters $a^{\prime }=1.0667,b^{\prime }=1.1248,c^{\prime }=14.92,d^{\prime }=254.70,e=581,g=$ 23,978, $f=1,\mathrm{L}=5000$, and $\ell =859$, the S-box constructed by the proposed algorithm is shown in

Table 3. The S-box constructed by the proposed method.

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

. The proposed algorithm is important in the sense that an S-box is guaranteed for each value of ℓ. Consequently, the total number of S-boxes for the chosen parameters is equal to the number of values of $\ell .$ The time complexity of the proposed S-box is given by the following:

Lemma 2.

Suppose that all the inputs $\mathrm{L}$, $a^{\prime },b^{\prime },c^{\prime },d^{\prime },e,g$, f, m, and integer $\ell \le |{\mathsf{\Lambda }}_{\mathrm{M}}|$ are known and the size of the set ${}_{c^{\prime }}{\mathrm{B}}_{d^{\prime }}$ is η. Then, the time complexity for the proposed S-box is $\mathcal{O}\left(\max (n\eta \mathrm{L},2^{2m})\right)$.

Proof. 

For given inputs the enumeration of u triads the S-box generator needs $\mathcal{O}\left(n\eta \mathrm{L}\right)$ time. Further, the arrangements of triads according to the ordering ${\prec }_{\theta }$ takes $\mathcal{O}\left(\left(u\right)\log \left(u\right)\right)$ time. The time taken by the execution of nested two loops is $\mathcal{O}\left(2^{2m}\right)$. As $2^{2m}>\left(u\right)\log \left(u\right)$, therefore, the time complexity of the proposed S-box algorithm is $\mathcal{O}\left(\max \{n\eta \mathrm{L},2^{2m}\}\right)$. □

4.2. Security Analysis

To test the cryptographic strength of an S-box obtained by the new technique, several standard security tests are applied, and the results are compared with some well-known schemes. The analysis is given as follows.

4.2.1. Linear Attacks

Nonlinearity (NL): The idea of NL was first proposed in [54], which determines the ability of an S-box to create confusion and diffusion in a plaintext. For an $8\times 8$ S-box $\sigma$ over $\mathrm{GF}\left(2^8\right)$, the NL is defined as

$$\mathrm{NL}\left(\sigma \right)=\underset{x,y,z}{\min }\{\alpha \in \mathrm{GF}\left(2^8\right)\mid x·\sigma \left(\alpha \right)\ne y·\alpha \oplus z\},$$

where $x\in \mathrm{GF}\left(2^8\right)\backslash \left\{0\right\}$, $y\in \mathrm{GF}\left(2^8\right)$, $z\in \mathrm{GF}\left(2\right)$, and operation “.” is an inner product over $\mathrm{GF}\left(2\right)$. An S-box has high resistance against linear attacks if it has high NL. However, Meier and Staffelbach [55] demonstrated that an S-box with high NL may lack some other cryptographic features. Therefore, it is crucial to design an S-box that has the best NL and passes other security performance tests. The NLs of all component functions of the proposed S-box are given in Figure 4. The minimum NL of the proposed S-box is 106, which is sufficient to resist powerful linear attacks.

Algebraic complexity (AC): For an S-box, the concept of a linear polynomial was first proposed in [56]. The AC represents the number of non-zero terms in a linear polynomial of an S-box. An S-box can have a maximum AC value of 255. For the proposed S-box, the AC is 254, which is very near to the optimal value. This shows that the proposed S-box is very strong against algebraic attacks. In

Table 4. The coefficients of the linear polynomial for the proposed S-box.

0	6	42	166	207	41	47	29	248	71	175	56	89	156	93	154
241	241	180	24	18	184	167	158	9	166	24	197	4	169	90	189
85	126	114	51	130	26	66	202	86	218	16	69	73	203	194	155
151	217	132	51	41	126	162	153	33	144	181	249	81	247	218	68
130	220	111	70	184	240	112	99	115	56	247	188	245	20	172	57
207	30	209	22	68	44	140	7	24	94	151	135	247	161	0	229
119	52	6	99	233	126	100	52	52	121	123	192	172	4	131	57
251	18	119	193	98	39	51	202	26	195	123	230	134	251	137	70
191	209	105	173	146	39	35	233	72	32	90	58	177	194	177	88
139	47	119	82	101	46	234	180	18	60	18	75	16	28	226	155
38	217	248	78	108	197	36	225	149	233	144	153	123	83	30	152
222	94	17	211	25	54	20	83	217	143	121	44	45	248	236	244
130	76	104	161	8	106	117	57	93	46	247	29	47	131	26	210
132	73	200	219	154	199	19	94	109	110	142	220	90	199	31	18
116	212	147	189	50	252	88	163	116	137	213	98	125	31	74	49
82	227	83	131	238	5	60	170	176	245	221	174	180	51	154	2

, the coefficients of the linear polynomial are shown.

Linear approximation probability (LAP): Mitsui [57] was the first to introduce the LAP test. The LAP is determined by the number of bits of plaintext and ciphertext that overlap. The formula that computes the LAP is

$$\mathrm{LAP}\left(\sigma \right)={\displaystyle \frac{1}{2^8}}\underset{x,y}{\max }|\#\{\alpha \in \mathrm{GF}\left(2^8\right)\mid x·\alpha =y·\sigma \left(\alpha \right)\}-2^7|$$

where $x\in \mathrm{GF}\left(2^8\right)$, $y\in \mathrm{GF}\left(2^8\right)\backslash \left\{0\right\}$. An S-box has high resistance against linear attacks if it has low LAP. The proposed S-box has the LAP value 0.133, which is close to the optimal value and shows that the proposed S-box is secure against linear attacks.

4.2.2. Differential Attacks

Differential approximation probability (DAP): In 1991, the concept of the DAP was introduced in [58]. The DAP is used to test the resistance of an S-box against differential attacks. For an S-box $\sigma$ of size $8\times 8$, the DAP is measured by

$$\mathrm{DAP}\left(\sigma \right)=\frac{1}{2^8}\underset{△x,△y}{\max }\#\left\{x\in \mathrm{GF}\left(2^8\right)\mid \sigma \left(x\right)\oplus \sigma (x+△x)=△y\right\},$$

where $△x\in \mathrm{GF}\left(2^8\right)\backslash \left\{0\right\},△y\in \mathrm{GF}\left(2^8\right)$, and the operation ⊕ is bit-wise addition over $\mathrm{GF}\left(2^8\right)$. Cryptographically, an S-box has high security against differential attacks if its DAP value is close to zero. The DAP of the proposed S-box is 0.039, which is very low. Consequently, the proposed S-box has high security against approximation attacks.

4.2.3. Analysis of Boolean Functions

Strict avalanche criteria (SAC): The SAC test [59] is a basic criterion to check the ability of an S-box to create diffusion in a plaintext. The SAC calculates the changing effect of output bits when a single input bit has been changed. The values for SAC of an $8\times 8$ S-box $\sigma$ is computed by matrix A($\sigma$)$=\left[m_{ij}\right]$ for

$$m_{ij}=\frac{1}{2^8}\sum _{x\in \mathrm{GF}\left(2^8\right)}w\left({\sigma }_i(x\oplus {\alpha }_j)\oplus {\sigma }_i\left(x\right)\right),$$

where $w\left(y\right)$ is the hamming weight of $y\in \mathrm{GF}\left(2^8\right)$, ${\sigma }_i$ and ${\sigma }_k$ are ith and kth Boolean functions of $\sigma$, respectively, and $1\le i,j,k\le 8$. The ideal value for SAC is $0.5$. If the calculated value is closer to 0.5, then it means that an S-box fulfills the SAC criterion. The dependence matrix of the proposed S-box is shown in

Table 5. Dependence matrix with entries $m_{ij}$.

0.424	0.563	0.578	0.531	0.500	0.516	0.500	0.531
0.484	0.578	0.469	0.469	0.500	0.5	0.516	0.453
0.422	0.547	0.500	0.484	0.563	0.547	0.516	0.422
0.469	0.547	0.500	0.563	0.500	0.469	0.500	0.516
0.500	0.547	0.500	0.531	0.563	0.484	0.500	0.500
0.500	0.531	0.516	0.531	0.578	0.547	0.578	0.531
0.563	0.500	0.438	0.500	0.547	0.469	0.469	0.500
0.516	0.453	0.531	0.453	0.453	0.563	0.531	0.531

, where the minimum and maximum values of the SAC are $0.422$ and $0.578$, respectively. Based on these values, we can say that the proposed S-box meets the SAC criterion.

Bit independence criterion (BIC): The BIC test [59] is used to determine how independent a pair of output bits is when one input bit is inverted. The diffusion-creating ability of an S-box is also determined by the BIC criterion. It is found by computing the dependence matrix $B\left(\sigma \right)=\left[d_{ij}\right]$, where $d_{ij}$ is calculated by

$$d_{ij}=\frac{1}{2^8}\sum _{\underset{1\le k\le 8,k\ne i}{x\in \mathrm{GF}\left(2^8\right)}}w\left({\sigma }_i(x\oplus {\alpha }_j)\oplus {\sigma }_i\left(x\right)\oplus {\sigma }_k(x\oplus {\alpha }_j)\oplus {\sigma }_k\left(x\right)\right).$$

The requirement of the BIC analysis is that all values of $d_{ij}$ should be approximately equal to 0.5. It can be observed that each $d_{ij}$ ranges between $0.473$ and $0.531$. This means that the S-box satisfies the BIC criterion. Furthermore,

Table 6. Dependence matrix with entries $d_{ij}$.

—	0.488	0.484	0.488	0.477	0.525	0.531	0.523
0.488	—	0.5	0.490	0.508	0.504	0.516	0.473
0.484	0.5	—	0.486	0.522	0.506	0.510	0.488
0.488	0.490	0.486	—	0.512	0.504	0.494	0.498
0.477	0.508	0.522	0.512	—	0.498	0.522	0.5
0.525	0.504	0.506	0.504	0.498	—	0.486	0.527
0.531	0.516	0.510	0.494	0.522	0.486	—	0.514
0.523	0.473	0.488	0.498	0.5	0.527	0.514	—

indicates that the values of each element $d_{ij}$ of the correlation matrix of $x_i\oplus x_j$ for all input $x_i\in \mathrm{GF}\left(2^8\right)$$(i\ne j,$ where $i,j=1,2,\dots ,8)$ of the given S-box are all close to 0.5, which shows that the S-box meets the BIC. To determine the BIC results for NL, we could calculate NL for each output bit of ($x_i\oplus x_j$) for all input $x_i\in \mathrm{GF}\left(2^8\right)$, where $i,j=1,2,\dots ,8$. A bar chart of BIC-NL for the proposed S-box is shown in Figure 5, which reveals that the newly designed S-box satisfies the BIC criterion.

4.2.4. Alteration in S-boxes

To have a sufficient cryptographic strength, the S-box construction technique should be capable of constructing a number of variant S-boxes [5]. This is because many cryptosystems require more than one secure S-box. We took a fixed set of u triads enumerated by the parameters noted in Section 4.1. Since for each value of ℓ, an S-box is guaranteed, we picked all the corresponding S-boxes for 1319 randomly chosen values of ℓ and computed the NL of each S-box. We found that the total number of distinct S-boxes is 1256, which is $95\%$ of all the S-boxes. Further, the distribution of values of ℓ is shown in Figure 6a, and the behavior of the NL of the constructed S-boxes is illustrated by Figure 6b. More explicitly, the jth value in Figure 6b represents the minimum NL of the S-box generated by the randomly chosen jth value of the parameter ℓ in Figure 6a. The fluctuation in the NL values is evident from Figure 6b, which implies a variation in the associated S-boxes.

Furthermore, Figure 6b shows that the minimum NL for most of the S-boxes oscillates between 90 and 100. However, there exists a large number of S-boxes with NL greater than 100. Thus, the above arguments explain that the proposed method is not only capable of constructing a number of distinct S- boxes but also has the capability to construct highly nonlinear S-boxes.

4.3. Performance Comparison

In this part, we compare our newly designed S-box with some other S-boxes constructed by different methods, including chaotic maps and elliptic curves [5,6,7,31,33,34,35,36,37,38,39,40,41,42,43,44,47,48,60,61]. The NL comparison of the proposed S-box with other S-boxes in

Table 7. Comparison of S-boxes across a range of methods (rows) against our method based on triads. The comparison is made over nine different indicators (columns). For each indicator (column), we present in boldface the methods that outperform our method.

S-boxes	Method	NL			SAC		BIC			DAP
		NL	AC	LAP	Min	Max	Min	Max	NL
Ref. [33]	Chaos	100	255	0.129	0.422	0.594	0.477	0.525	98	0.039
Ref. [34]		96	254	0.023	0.391	0.625	0.477	0.531	92	0.050
Ref. [35]		103	255	0.132	0.398	0.570	0.472	0.535	96	0.039
Ref. [36]		100	255	0.152	0.391	0.586	0.468	0.537	100	0.039
Ref. [37]		98	254	0.133	0.422	0.609	0.477	0.535	94	0.054
Ref. [31]		104	255	0.133	0.359	0.609	0.457	0.535	96	0.039
Ref. [38]		104	254	0.133	0.438	0.641	0.475	0.547	98	0.039
Ref. [39]		96	255	0.125	0.422	0.609	0.471	0.547	96	0.047
Ref. [46]	Other	112	253	0.063	0.453	0.563	0.480	0.525	112	0.016
Ref. [44]		104	253	0.141	0.406	0.594	0.461	0.522	98	0.054
Ref. [40]		102	254	0.133	0.359	0.562	0.467	0.535	98	0.039
Ref. [41]		102	254	0.148	0.375	0.609	0.470	0.521	100	0.039
Ref. [60]		104	255	0.109	0.391	0.593	0.454	0.49	102	0.047
Ref. [42]		104	255	0.094	0.391	0.578	0.476	0.529	103	0.023
Ref. [61]		96	255	0.125	0.359	0.609	0.477	0.541	98	0.039
Ref. [43]		104	254	0.148	0.438	0.578	0.482	0.543	96	0.047
Ref. [5]	Elliptic	106	255	0.148	0.406	0.641	0.471	0.537	98	0.039
Ref. [6]	Curves	106	253	0.188	0.406	0.609	0.465	0.527	98	0.039
Ref. [47]		106	255	0.148	0.438	0.594	0.465	0.545	98	0.039
Ref. [48]		104	255	0.145	0.391	0.625	0.471	0.531	98	0.039
Ref. [7]		106	255	0.148	0.406	0.625	0.471	0.539	96	0.047
Our Method	Triads	106	254	0.133	0.422	0.578	0.473	0.531	100	0.039

shows that our S-box has better NL than the S-boxes in papers [31,33,34,35,36,37,38,39,40,41,42,43,44,48,60,61] and, consequently, it has better resistance and security against linear attacks as compared to the S-box with lower NL. The comparison reveals that the newly designed S-box has lower LAP than the S-boxes constructd by the methods [5,6,7,36,41,43,44,47,48]. Similarly, the DAP results of the proposed S-box is better than the DAP results of the S-boxes in [7,34,37,39,43,44,60]. Hence, our S-box is more secure against approximation attacks. The SAC values of the newly constructed S-box ranges from $0.421$ to $0.578$, which are very close to the ideal value $0.5$. From Table 7, it is clear that the proposed S-box has better SAC results as compared to that of S-boxes in [5,6,7,31,34,35,36,37,38,39,40,41,42,48,60,61]. The S-boxes in Table 7 have satisfactory BIC values, but the proposed S-box has better BIC-NL than that of the S-boxes in [5,6,7,31,33,34,35,37,38,39,40,43,44,47,48,61]. The AC of the proposed S-box is higher than the AC of S-boxes in [6,44] and comparable with other schemes. There are two methods that outperform our method in more than four indicators: Ref. [42], which outperforms our method in four indicators, while in two indicators (AC and max BIC), it outperforms our method by a very small amount. Ref. [46] outperforms our method in all indicators except algebraic complexity.

This performance comparison of our proposed scheme based on triads against already established methods reveals that the proposed scheme has the capability to design highly secure S-boxes when compared with other schemes with different underlying mathematical structures such as chaotic maps, elliptic curves, and some other algebraic structures.

5. Conclusions

In this paper, we have developed a new geometric method to enumerate all distinct resonant triads in a given box of wavenumbers. Considering aspect ratio $f=1$, we computed all triads in a specific region and observed that the new method is very efficient for obtaining all triads in the grid of wavenumber size 5000 when compared with the state of the art [4,16]. The new method can also be used to obtain quasi-resonant triads. For aspect ratios $f\ne 1$, we proved numerically that the proposed method is capable of enumerating quasi-resonant triads in any region. As an application, we have defined a new total order on the set of triads to subsequently design an S-box generator via a novel algorithm whose time complexity we prove mathematically.We show that the proposed S- box generator is useful in cryptosystems using a number of S-boxes. Analysis of the S-box generator shows that our generator is capable of enumerating a highly secure S- box: considering nine key indicators based on nonlinearity, algebraic complexity, linear and differential approximation probabilities, strict avalanche criteria, and bit independence criterion, we compared our method based on triads against 20 state-of-the-art methods based on elliptic curves, chaos, and other algebraic methods. We outrank 18 of these methods [5,6,7,31,34,35,36,37,38,39,40,41,43,44,47,48,60,61] in the majority of the nine key indicators, and we basically draw with the methods from Refs. [33,42].

Our future goal is to improve our parametrization so as to further reduce the time complexity of the resonant triad search algorithm. Moreover, we will extend the current scheme of the S-box enumeration to construct a random number generator.