## Abstract

The role of Boolean functions is prominent in several areas including cryptography, sequences, and coding theory. Therefore, various methods for the construction of Boolean functions with desired properties are of direct interest. New motivations on the role of Boolean functions in cryptography with attendant new properties have emerged over the years. There are still many combinations of design criteria left unexplored and in this matter evolutionary computation can play a distinct role. This article concentrates on two scenarios for the use of Boolean functions in cryptography. The first uses Boolean functions as the source of the nonlinearity in filter and combiner generators. Although relatively well explored using evolutionary algorithms, it still presents an interesting goal in terms of the practical sizes of Boolean functions. The second scenario appeared rather recently where the objective is to find Boolean functions that have various orders of the correlation immunity and minimal Hamming weight. In both these scenarios we see that evolutionary algorithms are able to find high-quality solutions where genetic programming performs the best.

## 1 Introduction

In cryptography, one standard division is into symmetric key cryptography and public key cryptography (Diffie and Hellman, 1976; Paar and Pelzl, 2010). Going one step further, symmetric key cryptography can be again divided into block ciphers and stream ciphers. A common trait for all these ciphers is that they are designed in accordance with the number of cryptographic criteria they need to fulfill. These varying criteria enable ciphers to resist various cryptanalysis attacks. Some of the most common attacks on block ciphers are differential (Biham and Shamir, 1991) and linear (Matsui and Yamagishi, 1993) cryptanalysis. While on stream ciphers, attacks are the Berlekamp--Massey (Massey, 1969), fast correlation (Meier and Staffelbach, 1988), algebraic (Courtois and Meier, 2003), and fast algebraic (Courtois, 2003) cryptanalysis, and additionally on combiner generators, the correlation cryptanalysis (Siegenthaler, 1985).

One well researched source of nonlinearity in ciphers is Boolean functions. In block ciphers, the nonlinearity often comes from Substitution Boxes (S-Boxes) which are actually a number of Boolean functions (hence also the name vectorial Boolean functions). On the other hand, in stream ciphers that nonlinearity comes from the Boolean functions. Both of those scenarios, while not the only ones, show us the important role of Boolean functions in cryptography. Finding Boolean functions fitting all the criteria and analyzing the best possible trade-offs between these criteria are still crucial questions today.

Historically, Boolean functions were predominantly used in conjunction with Linear Feedback Shift Registers (LFSRs). Two models that were most often used were filter generators and combiner generators. In a combiner generator, several LFSRs are used in parallel and their output is the input for a Boolean function. On the other hand, in a filter generator, the output is obtained by a nonlinear combination of a number of positions in one longer LFSR (see, e.g., Carlet, 2010) because Boolean functions need to be balanced, with high nonlinearity, large algebraic degree, large algebraic immunity, large fast algebraic immunity, and large correlation immunity (in the case of combiner generators).

We said that a Boolean function needs to be balanced (among other criteria) to be suitable for cryptography. Indeed this is true, but only when we consider the role of Boolean functions in filter and combiner generators. However, recently one more application emerged where we are actually interested in Boolean functions that have the minimal Hamming weight and are therefore as far as possible from being balanced. Such Boolean functions can be used to help resist side-channel attacks. These attacks do not rely on the security of the underlying algorithm, but rather on the implementation of the algorithm in a device (Mangard et al., 2007). One class of countermeasures against side-channel attacks are masking schemes. In masking schemes one randomizes the intermediate values that are processed by the cryptographic device. One obvious drawback of such an approach is the masking overhead which can be substantial in embedded devices or smart cards. Correlation immune Boolean functions can reduce the masking overhead either by applying a leakage squeezing method (Carlet et al., 2012; Carlet and Guilley, 2014) or with Rotating S-Box masking (Carlet and Guilley, 2013).

In order to obtain adequate Boolean functions, there exists a number of construction methods. Such methods can be roughly divided into algebraic constructions, random search, heuristics, and combinations of those methods (Picek, Jakobovic et al., 2015). In this article, we examine one branch of heuristics, more precisely Evolutionary Algorithms (EAs) in order to evolve Boolean functions in accordance with the aforementioned properties. It is worth noting that EAs can be used either as the primary or the secondary construction method as will be shown in the next section. In primary constructions one obtains new functions without using known ones. In secondary constructions, one uses already known Boolean functions to construct new ones (either with different properties or sizes) (Carlet, 2010). Furthermore, we experiment with Boolean function sizes that are of practical importance. However, such practical sizes also mean that the search space is very large. For a Boolean function with inputs, there are in total possible Boolean functions.

We experiment with evolutionary algorithms that use three different representations: genetic algorithms (GAs) with binary encoding, genetic programming (GP) with tree representation, and finally, Cartesian genetic programming (CGP) with graph representation. Besides the aforesaid single-objective algorithms, we experiment with a multiobjective approach. We emphasize that here we present only two possible applications of EAs for cryptography. There are many other possibilities, for example, evolution of S-Boxes (Picek, Ege et al., 2014), finding search space parameters that can cause faults (Carpi et al., 2013), design of pseudorandom number generators (Lamenca-Martinez et al., 2006), and design of block ciphers (Hernandez-Castro et al., 2006), to name a few.

The remainder of this paper is organized as follows: Section 2 gives relevant information about Boolean functions, their representations and properties. In Section 3, we present a survey of related work. Next, in Section 4, we give our motivation for this research. In Section 4.1, we present the problem of finding minimal Hamming weight Boolean functions and their usage in masking schemes, and then in Section 4.2, we discuss the necessary properties of a Boolean function to be used in combiner or filter generators. Section 5 gives our experimental setup, the algorithms we consider, as well as the fitness functions for two design problems. Next, in Section 6, we give results for all algorithms and fitness functions as well as a short discussion and future work. Finally, in Section 7, we conclude the article.

## 2 Preliminaries

Let be positive integers, that is, . The set of all -tuples of elements in the field is denoted as where is the Galois field with two elements. The inner product of two vectors and is denoted by and equals . Here, “” represents addition modulo two (bitwise XOR). The Hamming weight (HW) of a vector , where , is denoted . It is the number of nonzero positions in the vector, that is, . The support of a vector , where , is the set of indices corresponding to the nonzero positions in the vector .

An -function is any mapping *F* from to . If equals 1, the function is called a Boolean function.

A Boolean function on can be uniquely represented by a truth table (TT), which is a vector that contains the function values of , ordered lexicographically, that is, .

The *support* () of a Boolean function is the set containing the nonzero positions in the truth table (TT) representation, that is, . The HW of a Boolean function is denoted by ; it is the cardinality of its support.

*balanced*if it takes the value 1 exactly the same number of times as the value 0 when the input ranges over . For the Walsh--Hadamard transform, a Boolean function is balanced if (Preneel et al., 1991):

*nonlinearity*of . The nonlinearity of a Boolean function can be expressed in terms of the Walsh--Hadamard coefficients as (Carlet, 2010):

*bias of nonlinearity*if it has the same output as its best affine approximation with probability of (Massey, 1995):

Boolean function is *-resilient* if it is balanced and with correlation immunity of order (Siegenthaler, 2006).

*algebraic degree*() of a Boolean function

*f*is defined as the number of variables in the largest product term of the function’s ANF having a nonzero coefficient (MacWilliams and Sloane, 1977): Here, is defined by the Möbius inversion principle.

The *algebraic immunity* () of a Boolean function is the lowest degree of a nonzero function from into for which or where and are Boolean functions. Here, is the Hadamard product of and , whose support is the intersection of the supports of and . A function such that is called an annihilator of (Carlet, 2010; Meier et al., 2004).

*fast algebraic immunity*(): For further information about Boolean functions, their properties and roles in cryptography, we refer interested readers to Carlet (2010) and MacWilliams and Sloane (1977).

## 3 Related Work

There exists a number of works that examine Boolean functions in cryptography and their generation with Evolutionary Computation (EC) techniques. Here, we try to enumerate them as exhaustively as possible.

As far as the authors know, the first application of genetic algorithms (GAs) to the evolution of cryptographically suitable Boolean functions was done by Millan et al. (1997) where they experimented with GA to evolve Boolean functions with high nonlinearity.

In Clark's 1998 thesis, he presents several applications of optimization techniques in the field of cryptology. One of the applications is the evolution of Boolean functions with high nonlinearity using GA and hill climbing techniques.

Millan et al. (1998) use GA to evolve Boolean functions that have high nonlinearity. In conjunction with the GA, they use hill climbing together with a resetting step in order to find Boolean functions with even higher nonlinearity and sizes of up to 12 inputs. They find balanced Boolean functions with 8 inputs that have nonlinearity 112 and correlation immunity equal to one. Furthermore, when using GA with hill climbing, they find eight input balanced Boolean functions with nonlinearity equal to 116. Millan et al. (1999) use variations of the hill climbing method in order to find Boolean functions that have high nonlinearity and low autocorrelation.

Clark and Jacob (2000) experiment with two-stage optimization to generate Boolean functions. They use a combination of simulated annealing (SA) and hill climbing with a cost function motivated by Parseval theorem in order to find functions with high nonlinearity and low autocorrelation.

Clark et al. (2002) use simulated annealing to generate Boolean functions with cryptographically relevant properties. In their work, they consider balanced function with high nonlinearity and with the correlation immunity property less than or equal to two.

Kavut and Yücel (2003) develop improved cost functions for a search that combines SA and hill climbing. With that approach, the authors are able to find some functions of eight and nine inputs that have a combination of nonlinearity and auto-correlation values previously unattained. They also experiment with a three-stage optimization method that combines SA and two hill climbing algorithms with different objectives.

Clark et al. (2003) experiment with SA in order to design Boolean functions by spectral inversion. They observe that many cryptographic properties of interest are defined in terms of Walsh--Hadamard transform values. Therefore, they work in the spectral domain where the cost function punishes those solutions that are not valid Boolean functions.

Burnett et al. (2004) present two heuristic methods where the goal of the first method is to generate balanced Boolean functions with high nonlinearity and low autocorrelation. The second method aims at generating resilient functions with high nonlinearity and algebraic degree that maximizes the Siegenthaler inequality.

Millan et al. (2004) propose a new adaptive strategy for the local search algorithm for generation of Boolean functions with high nonlinearity. Additionally, they introduce the notion of the graph of affine equivalence classes of Boolean functions.

Burnett (2005) in her thesis uses three heuristic techniques to evolve Boolean functions. The first method aimed to evolve balanced functions with high nonlinearity. The second method is used to find balanced Boolean functions with high nonlinearity that are correlation immune. The last method is used to find balanced functions with high nonlinearity and propagation characteristics different from zero. Furthermore, she experiments with the evolution of S-Boxes.

Aguirre et al. (2007) use a multiobjective random bit climber to search for balanced Boolean functions of size up to eight inputs that have high nonlinearity. Results indicate that the multiobjective approach is highly efficient when generating Boolean functions that have high nonlinearity.

Izbenko et al. (2008) use a modified hill climbing algorithm to transform bent functions to balanced Boolean functions with high nonlinearity.

McLaughlin and Clark (2013), on the other hand, use SA to generate Boolean functions that have optimal values of algebraic immunity, fast algebraic resistance, and algebraic degree. In their work, they experiment with Boolean functions of sizes up to 16 inputs.

Picek et al. (2013) experiment with GA and GP to find Boolean functions that possess several optimal properties. As far as the authors know, this is the first application of GP to the evolution of cryptographically suitable Boolean functions.

Hrbacek and Dvorak (2014) experiment with CGP to evolve bent Boolean functions of size up to 16 inputs. The authors investigate several configurations of algorithms in order to speed up the evolution process. Since they do not limit the number of generations, they succeed in finding bent function in each run for sizes between 6 and 16 inputs.

Several EAs are used by Picek, Batina et al. (2014) to evolve Boolean functions that have better side-channel resistance. This paper presents the first application of optimization techniques to Boolean functions with improved side-channel resistance. With the goal of finding maximal nonlinearity values of Boolean functions Picek, Marchiori et al. (2014) experiment with several EAs. Furthermore, they combine optimization techniques with algebraic constructions in order to improve the search. Although they are unable to find a balanced Boolean function with nonlinearity equal to 118, they present several possible avenues to follow when looking for highly nonlinear balanced Boolean functions. Picek, Jakobovic et al. (2015) compare the effectiveness of the CGP and GP approaches when looking for highly nonlinear balanced Boolean functions of eight inputs. Picek, Carlet et al. (2015) investigate several EAs in order to evolve Boolean functions with different values of the correlation immunity property. In the same paper, the authors also discuss the problem of finding correlation immune functions with minimal Hamming weight, but they experiment with only one size of Boolean functions. More extensive investigation on finding correlation immune Boolean functions with minimal Hamming weight and different sizes is conducted by Picek, Guilley et al. (2015).

Mariot and Leporati (2015b) use Particle Swarm Optimization (PSO) to find Boolean functions with good trade-offs of cryptographic properties for dimensions up to 12. The same authors use GAs where the genotype consists of the Walsh--Hadamard values in order to evolve semibent (plateaued) Boolean functions (Mariot and Leporati, 2015a).

Picek et al. (2016) conduct a detailed analysis of the efficiency of a number of evolutionary algorithms and fitness functions for Boolean functions with 8 inputs.

## 4 Motivation and Contributions

In this section, we present the reasoning behind the choice of the two scenarios explored in this article as well as our contributions.

### 4.1 Boolean Functions and Masking

#### 4.1.1 Masking with Codewords

Some applications manipulate sensitive data, such as cryptographic keys. Obviously, such data should remain secret. But skillful attackers might try to probe bits within a processor or a memory; they often succeed unless countermeasures are implemented.

In a view to protect secrets from probing attempts, it is customary to implement a countermeasure known as masking. It consists in changing randomly the representation of any so-called “sensitive” data, which depends on the key, to deceive the attacker. For example, if each bit , of a key is masked with a random bit , then an attacker can probe only . However, provided is uniformly distributed, the knowledge of does not disclose any information on bit , unless of course the attacker can also probe separately , in which case a higher-order masking would be necessary.

Now, for implementation reasons, it is often impractical to mask each bit individually. Indeed, the generation of unpredictable random numbers is costly. Hence, limiting the number of required random bits is to be welcomed. Furthermore, masking each bit of a vector of length would correspond to choosing the global mask in the whole set of possible masks, which may be too costly. Specifically, in the example of the AES cipher, some key bytes are mixed with plaintext bytes and enter a Substitution Box (S-Box). Generating all the 256 S-Boxes suitable for all possible masks in is expensive (especially for embedded systems). Hence, by restricting the number of possible masks, the overhead incurred by the countermeasure becomes more affordable.

Let us consider the simple example of masking one byte (). This can be achieved by using two complementary masks, such as and . An attacker who measures one bit of the masked byte cannot derive any information about the corresponding unmasked bit. If we define by the Boolean function whose support is , then plays its masking role as it is balanced: any bit can take value 0 and 1 with equal probability. In general, any Boolean function which is is a valid masking.

Let us now consider a stronger attacker who is able to probe two bits simultaneously. In this case, some information can be recovered. Typically, if the two masked bits are equal, then so are the two unmasked bits. So, by testing all pairs of bits, the attacker can recover the whole key (precisely: the whole key or its complement). It happens that, against such “second-order attacker,” it would be desirable that the masks be the support of a Boolean function. But clearly, this support must have a cardinality strictly greater than 2.

Hence, masking can be summarized as the problem of finding Boolean functions whose support is the masks’ set, with the two following constraints:

It should have small Hamming weight, for implementation reasons, and

it should have high correlation immunity to resist an attacker with multiple () probes.

Clearly, there is a trade-off which motivates the research for low Hamming weight, high correlation immunity Boolean functions.

#### 4.1.2 Example of Masking Suitable for the AES

An interesting example concerns the AES block cipher. It uses, up to composition by linear permutations, the same S-Box 16 times in each cipher round. Thus fast implementations (e.g., in hardware) instantiate 16 identical copies of the S-Box. If each S-Box is masked with a different mask and mapped randomly to each byte of the state, then the masking has no overhead. The question is thus to find masks as codewords of a code of length , dimension , and maximum dual distance. The solution is reached by a auto-dual code (this is the indicator of an 8-input Boolean function of Hamming weight 16) (Carlet and Guilley, 2013). Recall, a Boolean function is if and only if its support is a (not necessarily linear) code of dual distance strictly larger than (see, e.g., Carlet, 2010).

#### 4.1.3 State-of-the-Art and Open Issues

Some work on the topic of *low Hamming weight high correlation immunity* Boolean functions has been summarized by Hedayat et al. (1999) in their book. However, in this book, some entries (minimum Hamming weight for a given pair ) are expressed as nontight bounds. Recently, the exact value for entries corresponding to , and have been obtained by Carlet and Guilley (2014, Table 2, page 66), using a computer exhaustive search with a Satisfiability Modulo Tool (SMT).

Let be the lowest weight of Boolean functions with -bit inputs. Some remarkable values can be computed thanks to mathematical results:

When , then the function must be nonzero and have algebraic degree at most , that is, have even Hamming weight, and the indicator of the repetition code (i.e., of the pair of the all-0 and all-1 words) reaches this bound. Hence, .

According to the Siegenthaler bound or using the inverse Fourier transform formula, the only nonzero Boolean function is the constant function 1. Thus, .

A corollary of the Fon-Der-Flaas theorem states that for all and (see Carlet and Guilley, 2014, Corollary 4.7, page 61).

But for , there is no known mathematical procedure to compute .

Thus, the method to find in those cases consists of (i) deriving mathematically some lower bounds, and (ii) finding some codes for which the lower bound is tight.

Regarding the requirement (i) of the method, examples of lower bounds are as follows (Carlet and Guilley, 2014, Lemma 4.2, page 59):

for and , ;

for and , ;

for and , .

Now, condition (ii) can consist of selecting Boolean functions as indicators of binary codes. One solution consists of looking into databases, such as in Grassl (2007). However, codes in this database are linear, whereas some functions of minimal Hamming weight are indicators of nonlinear codes. This is the case, for example, for and (for which ) and for and (for which ). Obviously, the cardinality of support of a linear code is a power of two, which is not the case of 12 or 24.

Thus, it is required to find codes, possibly nonlinear, by unconventional methods. Indeed, as said, conventional constructions (Bose Ray--Chaudhuri Hocquenghem codes, Reed--Muller codes, Quadratic Residue codes, etc. (MacWilliams and Sloane, 1977)) favor linear codes, and furthermore do not have size minimality of codes (corresponding to Hamming weight minimality of Boolean functions) as a constraint. Thus, Boolean Satisfiability Problem (SAT) and SMT solvers are an option, but fail for problems with more than a few thousand variables (i.e., when ).

Now, it is well known that large S-Boxes make the design of block ciphers easier against cryptanalysis (Heys and Tavares, 1994). In addition, the fraction of good S-Boxes increases dramatically with the number of input variables (Youssef and Tavares, 1995). Even data bitwidths () are interesting from an implementation standpoint (for instance, Seberry et al. (1993) study a 12 bit input S-Box), whereas with odd data bitwidth () it is easier to attain S-Box of high quality (Detombe and Tavares, 1992). In software, the memory has become widely available. For instance, a table with entries is now not a practical issue. In hardware, the integration capability has also drastically improved. Besides, generic power-efficient (Bertoni et al., 2004) and area-efficient (Canright, 2005) methods have been developed, which makes the area-overhead in hardware more acceptable. Thus, it is interesting to study values of up to the reasonable value (for which the S-Box will have entries, which is large, but perfectly manageable in multimillion gate circuits).

This is a motivation for the choice of the first experimental scenario. For instance, the exact values that are of practical interest, remained unknown until this research. To conclude, our *first* optimization problem is to evolve Boolean functions that have certain order of correlation immunity and minimal Hamming weight. Furthermore, we do this for Boolean function of sizes between 11 and 13 inputs in an effort to find the exact remaining unknown values, but also for sizes 14 to 16 inputs to additionally investigate whether EAs can handle larger sizes.

### 4.2 On Combiner and Filter Generators

Recall from Section 1 that symmetric key cryptography is mainly divided into block and stream algorithms. In stream ciphers the encryption is made bitwise, through the addition, most often mod 2, of a *keystream* of the same size as the plaintext; output is by a pseudorandom generator (PRG) parameterized by a secret key. The resulting ciphertext can be decrypted by the same bitwise addition of the keystream, which gives back the plaintext. Stream ciphers are meant to be used on lighter devices than block ciphers and are supposed to be faster.

This is a difficult challenge for stream ciphers, since their security is dependent on the single choice of the PRG (they do not have the advantage of involving several rounds like block ciphers). In addition, nowadays the situation is still more difficult because modern block ciphers like the AES are very fast. We describe now briefly the criteria for the design of a stream cipher; more can be found in Carlet (2013). For reasons of efficiency, the PRG contains a linear part, initialized by means of a secret key and possibly a binary string sent more often over a public channel, called an initial vector, or initial value (IV). The keystream is made of the bits output by the PRG during a sufficient number of clock cycles, the linear part being updated at each clock cycle by a bijective linear function, and its content being nonlinearly combined in some way to produce the output bit. As observed already by Claude Shannon in the 1940s (Shannon, 1949), every cipher can be translated into a system of equations whose coefficients are deduced from the public data and from the observed data (the attacker being supposed to have access to a part of the keystream and needing to recover the rest of it); the unknowns of the equations can be the secret key bits or (in the case of stream ciphers) the initialization of the PRG. Note that the description of the PRG itself is supposed to be public, only its initialization being kept secret, since it is usually believed that making a part of the description secret, which requires the use of a part of the secret key for this description, reduces too much the possible size of the PRG and gives less robust systems than using the whole key-IV for generating the initialization. These equations need to have many unknowns and to be nonlinear for ensuring resistance to the possible cryptanalyses using this method. This nonlinearity is often ensured thanks to a *Boolean function* which, in the so-called *filter model*, takes as input bits chosen at fixed positions in the linear part of the PRG at its current state ( being much smaller than the size of the linear part, for reasons of speed). In the simplest version of PRG, the function itself outputs the current bit of the keystream and the choice of the function must ensure resistance to all known attacks (in practice, additional devices such as memory cells are added, but this simplest model of PRG must have been proven already resistant to all known attacks). If a function is weak against one of the existing attacks, it does not have practical interest, even if it is very strong against all the other attacks. Note that mainstream products in the industry are based on other stream ciphers (e.g., the eStream finalists). However, in the governmental market, in very light weight cryptography, or in “internal cryptography” dedicated to the inner protection of a chip, the filter model is often preferred, because its design has been extensively analyzed and has a low gate count.

#### 4.2.1 Algebraic Attacks on Stream Ciphers

As observed by Meier et al. (2004), since implies , that is, (recall that the addition of such Boolean functions is mod 2), the existence of and of algebraic degrees at most , such that , is equivalent to the existence of of algebraic degree at most such that (i.e., is an annihilator of ) or (i.e., is an annihilator of , this case happens when in the original equality ).

Other algebraic-like attacks exist.

The Rønjom--Helleseth attack (Rønjom and Helleseth, 2007) is efficient if the algebraic degree of the function is not close to and has complexity .

- The Fast Algebraic Attack (FAA), introduced by Courtois (2003) shortly after the standard algebraic attack is efficient if one can find of low algebraic degree and of algebraic degree significantly less than (but maybe larger than ) such that . Its complexity is roughly of the order (see Hawkes and Rose, 2004):

*Fast Algebraic Complexity*: whose value is invariant by changing into , is bounded above by , and is bounded below by the Fast Algebraic Immunity (Carlet, 2013; Carlet and Tang, 2015).

Let an -variable Boolean function , with (optimal) algebraic immunity . We assume it is used either as a combiner or as a filter in a stream cipher, and that the linear part of the pseudorandom generator (from which the Boolean function takes its input bits to produce a bit of the keystream at each clock-cycle) has size . Then the complexity of an algebraic attack using one annihilator of degree is roughly (Courtois and Meier, 2003). Let us choose (which is usual), then the complexity of the algebraic attack is at least (which is considered nowadays a sufficient complexity) for .

#### 4.2.2 Constraints and Bounds on Boolean Functions for Stream Ciphers

The function must be

*balanced*since otherwise the attacker would be able to distinguish a pair “plaintext-ciphertext” (or even a part of the plaintext and a part of the ciphertext, both at the same positions) from a randomly chosen pair by calculating the Hamming distance between the two texts of the pair.The function must have large

*algebraic degree*to allow resistance to the Berlekamp--Massey attack (Massey, 1969) and to the already mentioned Rønjom--Helleseth attack (Rønjom and Helleseth, 2007). In practice, we want an algebraic degree close to (the maximum for a balanced function).The function needs also to allow resistance to the fast correlation attack (Meier and Staffelbach, 1988). This attack uses the existence of an affine Boolean function (i.e., a function of algebraic degree at most 1) whose Hamming distance to is small; that is, to withstand this attack, the

*nonlinearity*should be high. The nonlinearity of an -variable Boolean function is bounded above by . The functions whose nonlinearity equals this maximal value are*bent*. A function used in a stream cipher should not be bent because bent functions are never balanced, but it should have nonlinearity near this maximum, since the fast correlation attack has an online complexity proportional to . Note that the average nonlinearity of random -variable functions lies asymptotically near (see Rodier, 2006), which is not bad. A large nonlinearity is analogous to the small bias of nonlinearity . The smallest possible bias is and is reached for bent functions. Since bent functions are unbalanced, we must accept larger bias, for instance approximately equal to or even ).The function must have

*algebraic immunity*close to , and fast algebraic complexity near ; for this, having a*fast algebraic immunity*near is sufficient. A difference by 1 in the algebraic immunity can have a large impact on the complexity of the algebraic attack (see Canteaut, 2006; Carlet, 2010).In some applications it is important to

*minimize the number of terms in the ANF*(see, e.g., Méaux et al., 2016).- A particular model of PRG is the combiner model, in which the input bits to the Boolean function are the outputs to independent Linear Feedback Shift Registers (LFSRs). A divide-and-conquer attack called a
*correlation attack*obliges the function to be -resilient for a sufficiently large value of (see Camion et al., 1992; Siegenthaler, 1985); that is, satisfy for every vector of Hamming weight at most . No class of sufficiently high order resilient functions having optimal algebraic immunity has been found yet. Moreover, the Siegenthaler bound shows that resilient functions cannot resist the Rønjom--Helleseth attack or the FAA. For a resilient Boolean function where and , the following Siegenthaler bound holds (Siegenthaler, 2006): Therefore, the correlation immunity and the algebraic degree properties are conflicting and it is not possible to obtain a Boolean function with both properties optimal. Hence, the combiner model is not widely used nowadays and the filter model (where the linear part is in general made of a single LFSR generating a sequence of maximal period) is preferred.

The goal in our *second* experimental set is to evolve Boolean functions that possess the aforesaid properties and have dimensions of practical interest, that is, 13 inputs and higher.

## 5 Experimental Setup

In this section, we briefly present the algorithms we use as well as the experimental setup and fitness functions for our research scenarios. For all the experiments, we use Evolutionary Computation Framework (ECF) (Jakobovic, 2014). We note that ECF has all the algorithms and fitness functions we define here.

### 5.1 Genetic Algorithm

The GA represents the individuals as strings of bits which present truth tables of Boolean functions. We use a 3-tournament selection, where the worst from the 3 randomly selected individuals is eliminated (Eiben and Smith, 2003). A new individual is created by applying crossover to the remaining two and then a mutation with given probability (see Algorithm 1).

Mutation is selected uniformly at random between a simple mutation, where a single bit is inverted, and a mixed mutation, which randomly shuffles the bits in a randomly selected subset. When the balancedness property is included in the fitness function, we include a balanced simple mutation; in this operator two bits are randomly inverted if the solution is balanced, and only one bit otherwise. The crossover operators are one-point and uniform crossover, performed uniformly at random for each new offspring. For each of the objectives we experiment with population sizes of 50, 100, 500, and 1 000 and individual mutation probabilities of 0.1, 0.3, 0.5, 0.7, and 0.9. The mutation probability is used to select whether an individual would be mutated or not, and the mutation operator is executed only once on a given individual. For example, if the mutation probability is 0.7, then on average 7 out of every 10 new individuals will be mutated and one mutation will be performed on each of those individuals.

### 5.2 Genetic Programming

Genetic programming (GP) uses a representation where individuals are trees of Boolean primitives which are then evaluated according to the truth table they produce. The function set for GP in all experiments is OR, XOR, AND, XNOR, and AND with one input inverted. Terminals correspond to Boolean variables. GP uses the same selection as in Algorithm 1 with tournament size 3. The crossover is performed with five different tree-based crossover operators selected at random: a simple tree crossover with 90% bias for functional nodes, uniform crossover, size fair, one-point, and context preserving crossover (Poli et al., 2008). We use a single mutation type, a subtree mutation applied with 30% probability, and experiment with maximum tree depth sizes of 5, 7, 8, and 9 and population sizes of 200, 500, 1 000, and 2 000.

### 5.3 Cartesian Genetic Programming

Cartesian genetic programming (CGP) (Miller, 1999) solutions are directed graphs with Boolean primitives as nodes, that are also evaluated using the truth table they produce. The function set for the CGP is the same as for the GP. The number of input connections for each node is two and the number of program output connections is one. Setting the number of rows to be 1 and levels-back parameter to be equal to the number of columns is regarded as the best and most general choice (Miller, 2011). We experiment with genotype sizes of 500, 1 000, 2 000, and 3 000 nodes and mutation rates of 1%, 4%, 7%, 10%, and 13% per node. For CGP individual selection we use a (1 + 4)-ES in which the offspring are favored over parents when they have a fitness better than or equal to the fitness of the parent. The population size for CGP equals five in all our experiments. The mutation operator is one-point mutation where the mutation point is chosen with a fixed mutation rate. The number of genes mutated is defined as a fixed percentage of the total number of genes. The single output gene is never mutated and is taken from the last node in the genotype. A gene chosen for mutation might be a node representing an input connection or a function.

### 5.4 Common Parameters

The number of independent runs for each experiment is 30. For the stopping condition we use the number of evaluations, which we set to 1 000 000.

### 5.5 Fitness Functions

#### 5.5.1 The First Design Problem

*maximization*:

Here, represents the Hamming weight of a Boolean function that has all ones in its truth table (i.e., , where represents the number of inputs of a Boolean function), is the order of the correlation immunity we want to obtain, and finally, represents the cardinality of the support of a Boolean function. The function consists of two parts: the first part rewards Boolean functions with smaller support, while the second part acts as a penalty for solutions with different target correlation immunity. The penalty part is multiplied with maximum value of the reward part, so that any solution with the right CI is always better than any other solution with CI different that the target value. This way, the distance to the target CI is regarded as the primary, and the support as the secondary objective.

Furthermore, we use a condition that eliminates trivial solutions, that is, Boolean functions with Hamming weight equal to either zero or . If such a solution is encountered, it is given the worst fitness value, which is .

#### Multiobjective Optimization

*minimized*, while the second criterion, , is

*maximized*.

In our experiments we apply the well known NSGA-II algorithm for multiobjective optimization (Deb et al., 2002). Note that NSGA-II can be paired with any of the Boolean representations (i.e., truth table in GA, tree in GP, and graph in CGP), but based on the performance in the initial round of experiments, we present only the results of the tree representation (GP) with the multiobjective evolution.

#### 5.5.2 The Second Design Problem

The second design problem is concerned with evolving Boolean functions with properties as described in Section 4.2.2. In this scenario we aim to optimize multiple properties that can oppose each other; therefore, we decide to apply a bottom up approach where we begin with a simple fitness function including only a single property, and then optimize more complex fitness functions that include additional properties of interest. We follow this approach since our previous results indicate that concentrating on a fitness function which includes all the properties can actually result in a solution with poor property values (Picek et al., 2013).

Before presenting fitness functions for the second design problem, we give a short discussion on two properties of interest: algebraic immunity and fast algebraic immunity. To evaluate a Boolean function with regards to those two properties is exponentially dependent on the number of inputs in computation time. Each of these properties for input dimension greater than 12 require several orders more time to calculate them (for instance, to calculate fast algebraic immunity of a Boolean function with 14 inputs we need more than 5 minutes). Therefore, it is not practical to include these two properties in our fitness functions since we work with 1 000 000 evaluations and 30 independent runs. Because of that, we calculate those two properties a posteriori and present the results. We note that this also presents one serious drawback of heuristics when evolving such Boolean functions, since it is unrealistic to run complex calculations for a great number of evaluations. Still, we show that evolving Boolean functions with regards to other properties can result in functions that are also relatively fit in regards to these two properties.

*maximization*. We start our experiments with the simplest fitness function that consists of one constraint (the function needs to be balanced) and the nonlinearity property. Here we use a two-stage fitness in which a fitness bonus equal to the nonlinearity is awarded only to a genotype that is perfectly balanced (this occurs when ); otherwise, the fitness is only the balancedness penalty. The balancedness penalty is defined as the negative difference up to the balancedness (i.e., the number of bits that need to be changed to reach balancedness). The delta function takes the value one when and is zero otherwise.

We note that to calculate the nonlinearity property, it is also possible to use the spectrum based cost function (Clark et al., 2004). However, we opted not to follow that line of work for the following reasons: the first one is that good values for additional parameters and for Boolean functions of sizes 13 up to 16 inputs are not well researched. The second reason lies in the fact that our experiments show that it is possible to evolve highly nonlinear Boolean functions even with the simplest calculation of the nonlinearity part of the fitness function.

#### Multiobjective Optimization

## 6 Results and Discussion

In this section, we present the results for each of the algorithms used as well as all fitness functions. The first part deals with the results for the evolution of minimal Hamming weight Boolean functions with different orders of correlation immunity. The second part presents Boolean functions for usage in combiner and filter generators. We report the performance of the selected algorithms, and additionally present the best obtained values in order to compare EAs with the existing results. The evolutionary algorithms are first compared with each other using basic statistical indicators to assess their performance. After that, we select only the single best results obtained by any of the algorithms and compare them with the values found in the related literature.

### 6.1 Tuning Phase

There is a large number of experiments; more precisely, for the first design problem, we have 81 combinations for each of the target CI values and function sizes (11–16 inputs) and 16 combinations for the second design problem (13–16 inputs, 4 fitness functions). In accordance with that, we conduct a parameter tuning phase based on the first problem for a medium-sized Boolean function of 14 inputs and the target correlation immunity order of four. Parameter tuning phase has a stopping condition of 1 000 000 evaluations. Afterward, we use the best obtained set of parameters for all test scenarios.

### 6.2 The First Design Problem

Here, we give results obtained for all algorithms and fitness functions as in Eq. (13). We display for each algorithm the best value it obtained for a certain input size and the correlation immunity level as well as the number of times it reaches that value. We do not give further statistics in the form of average values or similar since we do not consider it descriptive. As an example, consider a set of very good solutions and one extremely bad solution. Although that algorithm performs very well, the penalty given to that bad solution will skew the complete results. Furthermore, in this application domain the user is typically interested only in the best obtained solution.

#### 6.2.1 Genetic Algorithm

The results for GA in this application were very poor; in the tuning phase we were unable to obtain a single solution with the desired correlation immunity for any of the parameter settings. This suggests that bitstring GA with the common genetic operators we used cannot traverse the search space towards better solutions. The GA does succeed in finding the desired values, but only for very small problem sizes (e.g., for up to six variables), where the size of the solution is not large. However, since these cases are not representative of the problem, we do not experiment with GA in the rest of the article when considering the first design problem.

#### 6.2.2 Genetic Programming Results

In the tuning phase the results indicate the best maximum depth is 5; for different population sizes for GP there are practically no statistical differences, so we opted for the size of 1 000. Depending on the number of bits and the target CI, the problem at hand may be so easy that every combination of parameters reaches the same solution every time, or so hard that there is a very small probability of reaching the best (optimal) fitness value regardless of the parameters.

The results for the GP for all combinations of sizes and target CI are given in Table 1. Because of the previous remark, the results are given in the form of the best obtained value (i.e., the minimal Hamming weight) and the number of runs (out of 30) in which that value was reached.

t\n . | 11 . | 12 . | 13 . | 14 . | 15 . | 16 . |
---|---|---|---|---|---|---|

2 | 16/25 | 16/19 | 32/14 | 64/20 | 128/16 | 256/23 |

3 | 32/29 | 32/17 | 32/2 | 64/15 | 128/12 | 256/11 |

4 | 128/11 | 256/30 | 256/12 | 256/1 | 2 048/14 | 4 096/12 |

5 | 256/13 | 256/5 | 512/7 | 1 024/16 | 2 048/23 | 4 096/23 |

6 | 512/13 | 1 024/27 | 1 024/5 | 2 048/5 | 4 096/3 | 4 096/4 |

7 | 1 024/16 | 1 024/8 | 2 048/4 | 4 096/8 | 8 192/18 | 8 192/1 |

8 | 1 024/30 | 2 048/30 | 4 096/30 | 8 192/30 | 8 192/1 | 16 384/4 |

9 | 1 024/30 | 2 048/30 | 4 096/30 | 8 192/30 | 16 384/30 | 16 384/2 |

10 | 1 024/30 | 2 048/30 | 4 096/30 | 8 192/30 | 16 384/30 | 32 768/30 |

11 | 2 048/30 | 4 096/30 | 8 192/30 | 16 384/30 | 32 768/30 | |

12 | 4 096/30 | 8 192/30 | 16 384/30 | 32 768/30 | ||

13 | 8 192/30 | 16 384/30 | 32 768/30 | |||

14 | 16 384/30 | 32 768/30 | ||||

15 | 32 768/30 |

t\n . | 11 . | 12 . | 13 . | 14 . | 15 . | 16 . |
---|---|---|---|---|---|---|

2 | 16/25 | 16/19 | 32/14 | 64/20 | 128/16 | 256/23 |

3 | 32/29 | 32/17 | 32/2 | 64/15 | 128/12 | 256/11 |

4 | 128/11 | 256/30 | 256/12 | 256/1 | 2 048/14 | 4 096/12 |

5 | 256/13 | 256/5 | 512/7 | 1 024/16 | 2 048/23 | 4 096/23 |

6 | 512/13 | 1 024/27 | 1 024/5 | 2 048/5 | 4 096/3 | 4 096/4 |

7 | 1 024/16 | 1 024/8 | 2 048/4 | 4 096/8 | 8 192/18 | 8 192/1 |

8 | 1 024/30 | 2 048/30 | 4 096/30 | 8 192/30 | 8 192/1 | 16 384/4 |

9 | 1 024/30 | 2 048/30 | 4 096/30 | 8 192/30 | 16 384/30 | 16 384/2 |

10 | 1 024/30 | 2 048/30 | 4 096/30 | 8 192/30 | 16 384/30 | 32 768/30 |

11 | 2 048/30 | 4 096/30 | 8 192/30 | 16 384/30 | 32 768/30 | |

12 | 4 096/30 | 8 192/30 | 16 384/30 | 32 768/30 | ||

13 | 8 192/30 | 16 384/30 | 32 768/30 | |||

14 | 16 384/30 | 32 768/30 | ||||

15 | 32 768/30 |

#### 6.2.3 Cartesian Genetic Programming Results

When using CGP, we observed that the parameter values play a much more significant role than in the GP case. For Boolean functions up to a size of 13 inputs, the genotype size of 3 000 and mutation probability of 13% performs the best, but for larger sizes the best results are obtained for genotype size of 1 000 and mutation rate of 10%. In Table 2, we give the statistics for CGP in the form of the best obtained solutions and the number of times those solutions were reached. When the algorithm is not able to find a correct solution (by correct, we mean with a desired order of CI), we denote it with “-”.

t\n . | 11 . | 12 . | 13 . | 14 . | 15 . | 16 . |
---|---|---|---|---|---|---|

2 | 16/4 | 32/7 | 32/8 | 64/3 | 128/4 | 256/2 |

3 | 32/1 | 64/2 | 128/2 | 256/2 | 2 048/8 | 1 024/1 |

4 | 256/2 | 512/3 | 512/1 | 2 048/3 | 4 096/2 | 8 192/1 |

5 | 512/2 | 1 024/4 | 2 048/7 | 4 096/3 | 8 192/1 | 8 192/1 |

6 | 1 024/25 | 2 048/26 | 4 096/28 | 8 192/27 | 16 384/17 | 32 768/8 |

7 | 1 024/23 | 2 048/24 | 4 096/27 | 8 192/23 | 16 384/8 | 32 768/3 |

8 | 1 024/22 | 2 048/20 | 4 096/21 | 8 192/20 | 16 384/4 | – |

9 | 1 024/14 | 2 048/12 | 4 096/15 | 8 192/8 | 16 384/1 | – |

10 | 1 024/5 | 2 048/8 | 4 096/9 | 8 192/2 | – | – |

11 | 2 048/1 | 4 096/5 | – | – | – | |

12 | 4 096/1 | – | – | – | ||

13 | – | – | – | |||

14 | – | – | ||||

15 | – |

t\n . | 11 . | 12 . | 13 . | 14 . | 15 . | 16 . |
---|---|---|---|---|---|---|

2 | 16/4 | 32/7 | 32/8 | 64/3 | 128/4 | 256/2 |

3 | 32/1 | 64/2 | 128/2 | 256/2 | 2 048/8 | 1 024/1 |

4 | 256/2 | 512/3 | 512/1 | 2 048/3 | 4 096/2 | 8 192/1 |

5 | 512/2 | 1 024/4 | 2 048/7 | 4 096/3 | 8 192/1 | 8 192/1 |

6 | 1 024/25 | 2 048/26 | 4 096/28 | 8 192/27 | 16 384/17 | 32 768/8 |

7 | 1 024/23 | 2 048/24 | 4 096/27 | 8 192/23 | 16 384/8 | 32 768/3 |

8 | 1 024/22 | 2 048/20 | 4 096/21 | 8 192/20 | 16 384/4 | – |

9 | 1 024/14 | 2 048/12 | 4 096/15 | 8 192/8 | 16 384/1 | – |

10 | 1 024/5 | 2 048/8 | 4 096/9 | 8 192/2 | – | – |

11 | 2 048/1 | 4 096/5 | – | – | – | |

12 | 4 096/1 | – | – | – | ||

13 | – | – | – | |||

14 | – | – | ||||

15 | – |

#### 6.2.4 Multiobjective Optimization Results

In this problem, the multiobjective (MO) approach did not prove competitive; the results for a medium-sized problem with 13 bits are given in Table 3. Table entries with a dash denote that the algorithm did not reach the target CI value in this case. Since the MO results were significantly worse, we do not include them in further comparisons.

#### 6.2.5 Best Obtained Solutions

In this section, we compile the list of the best obtained solutions over all algorithms; the results are given in Table 4. The values given in bold present the ones that were previously unknown. Regarding the individual algorithms, if a solution is found only with GP, we leave the background color of a cell white. The cells with gray background, on the other hand, denote solutions found both with GP and CGP.

t\n . | 11 . | 12 . | 13 . | 14 . | 15 . | 16 . |
---|---|---|---|---|---|---|

2 | 16 | 16 | 32 | 64 | 128 | 256 |

3 | 32 | 32 | 32 | 64 | 128 | 256 |

4 | 128 | 256 | 256 | 256 | 2 048 | 4 096 |

5 | 256 | 256 | 512 | 1 024 | 2 048 | 4 096 |

6 | 512 | 1 024 | 1 024 | 2 048 | 4 096 | 4 096 |

7 | 1 024 | 1 024 | 2 048 | 4 096 | 8 192 | 8 192 |

8 | 1 024 | 2 048 | 4 096 | 8 192 | 8 192 | 16 384 |

9 | 1 024 | 2 048 | 4 096 | 8 192 | 16 384 | 16 384 |

10 | 1 024 | 2 048 | 4 096 | 8 192 | 16 384 | 32 768 |

11 | 2 048 | 4 096 | 8 192 | 16 384 | 32 768 | |

12 | 4 096 | 8 192 | 16 384 | 32 768 | ||

13 | 8 192 | 16 384 | 32 768 | |||

14 | 16 384 | 32 768 | ||||

15 | 32 768 |

t\n . | 11 . | 12 . | 13 . | 14 . | 15 . | 16 . |
---|---|---|---|---|---|---|

2 | 16 | 16 | 32 | 64 | 128 | 256 |

3 | 32 | 32 | 32 | 64 | 128 | 256 |

4 | 128 | 256 | 256 | 256 | 2 048 | 4 096 |

5 | 256 | 256 | 512 | 1 024 | 2 048 | 4 096 |

6 | 512 | 1 024 | 1 024 | 2 048 | 4 096 | 4 096 |

7 | 1 024 | 1 024 | 2 048 | 4 096 | 8 192 | 8 192 |

8 | 1 024 | 2 048 | 4 096 | 8 192 | 8 192 | 16 384 |

9 | 1 024 | 2 048 | 4 096 | 8 192 | 16 384 | 16 384 |

10 | 1 024 | 2 048 | 4 096 | 8 192 | 16 384 | 32 768 |

11 | 2 048 | 4 096 | 8 192 | 16 384 | 32 768 | |

12 | 4 096 | 8 192 | 16 384 | 32 768 | ||

13 | 8 192 | 16 384 | 32 768 | |||

14 | 16 384 | 32 768 | ||||

15 | 32 768 |

### 6.3 The Second Design Problem

In this section, we first present the results only as the obtained values of the related fitness functions, to compare different algorithms. For each combination of input size and fitness function, we give the best obtained value and the average value. Afterward, we single out only the best solutions for each fitness and input size and show the values of the relevant cryptographic properties. Table 5 shows the theoretical optimal values (maximal) of the observed properties. We omit the correlation immunity property from the table since its value is conflicting with the algebraic degree property and equals 0 when algebraic degree reaches the optimal value. Furthermore, we do not give data for the minimal number of terms in the ANF representation since those data are currently unknown.

#### 6.3.1 Genetic Algorithm Results

For this problem we initially applied the GA only for Boolean functions with 13 inputs. The results are shown in Table 6 where for each fitness function we give the best and the mean obtained result from 30 runs. The GA is significantly worse than the other methods; this is not surprising, since already for 13 bits the size of the truth table is 8 192. Since the difference in the performance increases only for a larger number of inputs, we omit the other GA results in this section.

#### 6.3.2 Genetic Programming Results

Table 7 shows the best and mean fitness values for all the combinations of the number of bits and the defined fitness functions. It can be shown that all the results are significantly better than the GA results (the analysis is not included). One possible downside is that the evaluation takes considerably longer than the GA, since for each GP individual a complete truth table must be constructed to evaluate the relevant properties.

n\ fitness . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|

13 | 4 032/4 032 | 4 042/4 039.1 | 4 040/4 039.9 | 4 043/4 039.9 |

14 | 8 064/8 064 | 8 074/8 072.7 | 8 078/8 073.9 | 8 078/8 074.5 |

15 | 16 256/16 256 | 16 264/16 263.8 | 16 265/16 264.9 | 16 265/16 264.5 |

16 | 32 512/32 512 | 32 615/32 523.875 | 32 616/32 524.1 | 32 526/32 523.2 |

n\ fitness . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|

13 | 4 032/4 032 | 4 042/4 039.1 | 4 040/4 039.9 | 4 043/4 039.9 |

14 | 8 064/8 064 | 8 074/8 072.7 | 8 078/8 073.9 | 8 078/8 074.5 |

15 | 16 256/16 256 | 16 264/16 263.8 | 16 265/16 264.9 | 16 265/16 264.5 |

16 | 32 512/32 512 | 32 615/32 523.875 | 32 616/32 524.1 | 32 526/32 523.2 |

#### 6.3.3 Cartesian Genetic Programming Results

We give the best and mean values for the CGP algorithm in Table 8. It can be seen there are virtually no differences between fitness functions 2, 3, and 4. This indicates that the significance of the ANFMinimize and CI properties in the fitness function is not enough to drive the search. To rectify this, one can add the weight factor on the ANFMinimize and the correlation immunity parts or restrict algebraic degree only to a certain level (and thus allowing the CI property to assume higher values).

n\ fitness . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|

13 | 4 032/3 981.6 | 4 042/3 988.8 | 4 042/3 988.8 | 4 042/3 991.9 |

14 | 8 064/7 413.2 | 8 075/7 431.1 | 8 075/7 432.7 | 8 075/7 420.5 |

15 | 16 128/15 346.4 | 16 140/15 351.9 | 16 140/15 351.9 | 16 138/15 364.1 |

16 | 32 256/31 636.6 | 32 268/31 587.45 | 32 268/31 587.5 | 32 264/31 602.4 |

n\ fitness . | 1 . | 2 . | 3 . | 4 . |
---|---|---|---|---|

13 | 4 032/3 981.6 | 4 042/3 988.8 | 4 042/3 988.8 | 4 042/3 991.9 |

14 | 8 064/7 413.2 | 8 075/7 431.1 | 8 075/7 432.7 | 8 075/7 420.5 |

15 | 16 128/15 346.4 | 16 140/15 351.9 | 16 140/15 351.9 | 16 138/15 364.1 |

16 | 32 256/31 636.6 | 32 268/31 587.45 | 32 268/31 587.5 | 32 264/31 602.4 |

#### 6.3.4 Best Obtained Solutions

Finally, we combine all the single-objective optimization results and select only the best obtained solutions for every combination of input size and fitness function. For each selected solution, we present the observed cryptographic properties in Table 9. Apart from the number of terms in ANF representation (#ANF), all the properties are maximized. We do not explicitly write which algorithm obtained a certain solution since it is not always easy to select a single best solution. For instance, when and fitness function is as in Eq. (18), the best GP result is given in Table 9. However, CGP finds a solution that has the nonlinearity equal to 4 030, algebraic degree equal to 12 and equal to 61. They both have the same fitness value, but the choice of which Boolean function is better is dependent on the setting one considers.

n . | fitness . | . | . | . | . | . | . |
---|---|---|---|---|---|---|---|

13 | 1 | 4 032 | 5 | 28 | 0 | 4 | 5 |

2 | 4 032 | 10 | 55 | 0 | 3 | 4 | |

3 | 4 032 | 7 | 8 | 0 | 3 | 4 | |

14 | 1 | 8 064 | 10 | 50 | 0 | 3 | 4 |

2 | 8 064 | 13 | 11 | 0 | 3 | 4 | |

3 | 8 064 | 12 | 14 | 1 | 3 | 4 | |

15 | 1 | 16 256 | 5 | 33 | 0 | 3 | 4 |

2 | 16 256 | 8 | 25 | 0 | 3 | 4 | |

3 | 16 256 | 8 | 9 | 0 | 3 | 4 | |

16 | 1 | 32 512 | 12 | 99 | 0 | 4 | 6 |

2 | 32 512 | 14 | 12 | 0 | 3 | 4 | |

3 | 32 512 | 12 | 21 | 1 | 4 | 6 |

n . | fitness . | . | . | . | . | . | . |
---|---|---|---|---|---|---|---|

13 | 1 | 4 032 | 5 | 28 | 0 | 4 | 5 |

2 | 4 032 | 10 | 55 | 0 | 3 | 4 | |

3 | 4 032 | 7 | 8 | 0 | 3 | 4 | |

14 | 1 | 8 064 | 10 | 50 | 0 | 3 | 4 |

2 | 8 064 | 13 | 11 | 0 | 3 | 4 | |

3 | 8 064 | 12 | 14 | 1 | 3 | 4 | |

15 | 1 | 16 256 | 5 | 33 | 0 | 3 | 4 |

2 | 16 256 | 8 | 25 | 0 | 3 | 4 | |

3 | 16 256 | 8 | 9 | 0 | 3 | 4 | |

16 | 1 | 32 512 | 12 | 99 | 0 | 4 | 6 |

2 | 32 512 | 14 | 12 | 0 | 3 | 4 | |

3 | 32 512 | 12 | 21 | 1 | 4 | 6 |

Clearly, the single-objective fitness gives insufficient weight to the CI property to be able to reach higher values. This can be remedied with the multiobjective optimization, as described next.

#### 6.3.5 Multiobjective Optimization Results

In order to compare the attained property values, we present the multiobjective results after the best results from single-objective optimization. Table 10 shows the best nondominated results from multiobjective optimization. The bold values in each row represent the property whose value was the best obtained, with other properties selected as large as possible.

n | |||

13 | 4 032 | 7 | 0 |

13 | 3 966 | 12 | 0 |

13 | 0 | 1 | 11 |

14 | 8 064 | 7 | 1 |

14 | 7 806 | 13 | 0 |

14 | 4 096 | 2 | 11 |

15 | 16 128 | 4 | 1 |

15 | 15 728 | 13 | 0 |

15 | 0 | 1 | 11 |

16 | 32 512 | 7 | 0 |

16 | 30 660 | 14 | 1 |

16 | 0 | 1 | 12 |

n | |||

13 | 4 032 | 7 | 0 |

13 | 3 966 | 12 | 0 |

13 | 0 | 1 | 11 |

14 | 8 064 | 7 | 1 |

14 | 7 806 | 13 | 0 |

14 | 4 096 | 2 | 11 |

15 | 16 128 | 4 | 1 |

15 | 15 728 | 13 | 0 |

15 | 0 | 1 | 11 |

16 | 32 512 | 7 | 0 |

16 | 30 660 | 14 | 1 |

16 | 0 | 1 | 12 |

In this problem the MO approach did succeed in finding larger values for certain properties, such as correlation immunity. The best values for nonlinearity and algebraic degree are about the same as in the single-objective case. If we are aiming for a more general all-round function with majority of properties being “good enough,” the single-objective is still the best option. However, if we need to maximize a single property, then the MO approach gives a viable alternative.

### 6.4 Discussion and Future Work

When considering the first design problem and the performance of CGP, we observe that the results are somewhat worse than those from related work (Picek, Carlet et al., 2015; Picek, Jakobovic et al., 2015). Actually, CGP often outperforms the GP on Boolean problems and one can ask why is this not the case here. Naturally, it is also hard to expect that CGP will always outperform GP because there is a high influence of the stopping condition and fitness function. We observed in our experiments that CGP would benefit from a higher number of evaluations in this problem, but to keep the comparison as fair as possible we retained the same number of evaluations as GP.

One can see that this function uses additional information and that is the target level of support (). Naturally, if one knows what level one wants to achieve, this does not pose a problem, but in the case where the best support value (and that is the value one is usually interested in) is unknown, this represents a serious drawback. Of course, it is possible to iteratively decrease the level of support until the algorithm cannot find such a value anymore, but that substantially prolongs the evolution process. However, we note that when using this fitness function, CGP easily outperforms GP. A plausible scenario for this fitness function could be to first use fitness as in Eq. (13) with GP to find out what is the best support value, and then to use fitness as in Eq. (24) and CGP to find as many as possible Boolean functions that have the desired support value.

When discussing the optimality of the results for the first design problem, we follow the conjectures from Bhasin et al. (2013). Here, represents the lowest weight of nonzero function of variables. The conjecture was made (Bhasin et al., 2013, Sec. C.2) that the values in each column of Table 4 are nondecreasing. The values for in Table 4 are interesting from this viewpoint: if the conjecture is true, then they are optimal since they cannot be smaller than for , but the conjecture may be false; further investigations are needed to clarify this point. If represent the minimal possible values, then since it is known from Bhasin et al. (2013, Sec. C.1) that , then the solution for and has also a minimal Hamming weight. Finally, for , by following the same reasoning, Hamming weights for are again the optimal values, if a solution for and has optimal values. Actually, the value was already known (see Hedayat et al., 1999, Table 12.1, page 319); hence it does not appear as a boldface entry in Table 4.

When examining the results for the second design problem, we can again see that the GP outperforms all other algorithms. As in the first problem, the GA is not able to produce any competitive results already for the smallest size that equals 13 inputs. The results for the CGP are somewhat worse than expected if considering some of the results from related work (Picek, Carlet et al., 2015; Picek, Jakobovic et al., 2015). However, those works consider only 8 inputs for Boolean functions so any comparison is hard to make. As in the first design problem, it seems that the CGP needs many more evaluations to compete with the GP on such large Boolean functions. Naturally, this should not pose a problem since CGP is much faster than the GP, but we refrain here from adding more evaluations in an effort to give as fair as possible a comparison. When evaluating the best obtained solutions, it is not always easy to say which of the solutions are the best ones. The same fitness values are produced by various combinations of properties so that depends on the setting one considers. Furthermore, we see that even the nonlinearity and the algebraic degree values are often far from the theoretical best values. This is interesting since for smaller sizes (e.g., 8 inputs) previous results suggest that EC is highly competitive and is able to reach theoretical bounds for those properties.

We note that the most of the evolved functions have highly suboptimal values for the AI and FAI properties, which suggests that these properties must be parts of the fitness function. However, as already stated, their computational complexity is too high to be used in evolutionary experiments when considering the number of necessary evaluations. Naturally, this does not mean that some of the evolved solutions do not have good values of AI and FAI properties, but only that it is unrealistic to expect that all evolved solutions are highly fit with regards to those properties. Furthermore, if the evolutionary process creates a large number of individuals that must be evaluated a posteriori for those two properties, then one is again facing the same problem since the high computational complexity could prevent such an analysis. When considering the minimization of the number of terms in the ANF representation, we see that EAs are able to minimize it significantly. Therefore, the main drawback when using the EAs to evolve balanced Boolean functions with large number of inputs stems from the complexity of certain parts of the evaluation and not from some inability of EAs to work with such sizes of Boolean functions.

## 7 Conclusion

In this article, we conducted an extensive analysis on the performance of EAs when evolving Boolean functions for cryptography. Furthermore, we consider only the Boolean function dimensions that are of practical relevance. In the first design problem, we evolve Boolean functions with minimal Hamming weight and different orders of CI properties. By doing so, we are able to find the previously unknown values for sizes 11, 12, and 13. In our opinion, 16 inputs is approaching the upper limit of what is manageable for EAs and larger sizes would pose a problem. However, EAs should definitely present a viable option when generating Boolean functions, especially since most of the algebraic constructions are not suitable for the problem at hand.

When considering the second design problem, the situation becomes more complicated. Since our fitness functions are now more complex, the evolution process lasts longer and it is a question whether one can consider 16 inputs as a manageable size. Furthermore, two of the properties that are relevant (namely, AI and FAI) are computationally too complex to be parts of the fitness function.

Naturally, one can evolve Boolean functions and disregard these properties in the fitness, but later choose the best overall solutions by a posteriori evaluating those criteria. We investigated that approach by iteratively making our fitness function more complicated (i.e., consisting of more terms), but without significant success. If the quality of the AI and FAI properties is essential, then we believe that algebraic constructions pose a better option since they require a smaller number of evaluations.

We also emphasize that this is the first work, as far as we are aware, that considers the minimization of ANF representation as a design criterion. Since this property has a big impact on the implementation cost, we believe it is of high importance and should be considered in the future. When looking at the big picture, we conclude that EAs are powerful enough to evolve Boolean functions that are of practical dimensions for usages in cryptography. However, at least by following the representations and fitness functions we consider, 16 inputs also seems to be close to the upper limit for EAs.

## Acknowledgments

This work has been supported in part by the Croatian Science Foundation under the project IP-2014-09-4882. In addition, this work was supported in part by the Research Council KU Leuven (C16/15/058), (CREA/14/005), and IOF project EDA-DSE (HB/13/020). Claude Carlet and Sylvain Guilley thank funding from ERA-NET CHIST-ERA “SECODE” (https://secode.enst.fr/).

## References

*Proceedings of the 2nd ACM Conference on Computer and Communications Security*