Jean-Pierre Lozi
email: rlozi@unice.fr

Oleg Garasym

René Lozi

J.-P Lozi

The Challenging Problem of Industrial Applications of Multicore-Generated Iterates of Nonlinear Mappings

Keywords: Chaos, Cryptography, Mappings, Chaotic pseudorandom numbers Attractors AMS Subject Classification 37N30, 37D45, 65C10, 94A60 p. 6

The study of nonlinear dynamics is relatively recent with respect to the long historical development of early mathematics since the Egyptian and the Greek civilization, even if one includes in this field of research the pioneer works of Gaston Julia and Pierre Fatou related to one-dimensional maps with a complex variable, nearly a century ago. In France, Igor Gumosky and Christian Mira began their mathematical researches in 1958; in Japan, the Hayashi' School (with disciples such as Yoshisuke Ueda and Hiroshi Kawakami), a few years later, was motivated by applications to electric and electronic circuits. In Ukraine, Alexander Sharkovsky found the intriguing Sharkovsky's order, giving the periods of periodic orbits of such nonlinear maps in 1962, although these results were only published in 1964. In 1983, Leon O. Chua invented a famous electronic circuit that generates chaos, built with only two capacitors, one inductor and one nonlinear negative resistance. Since then, thousands of papers have been published on the general topic of chaos. However, the pace of mathematics is slow, because any progress is based on strictly rigorous proof. Therefore, numerous problems still remain unsolved. For example, the long-term dynamics of the Hénon map, the first example of a strange attractor for mappings, remain unknown close to the classical parameter values from a strictly mathematical point of view, 40 years after its original publication. In spite of this lack of rigorous mathematical proofs, nowadays, engineers are actively working on applications of chaos for several purposes: global optimization, genetic algorithms, CPRNG (Chaotic Pseudorandom Number Generators), cryptography, and so on. They use nonlinear maps for practical applications without the need of sophisticated theorems. In this chapter, after giving some prototypical examples of the industrial

Introduction

The last few decades have seen the tremendous development of new IT technologies that incessantly increase the need for new and more secure cryptosystems.

For instance, the recently invented Bitcoin cryptocurrency is based on the secure Blockchain system that involves hash functions [START_REF] Delahaye | Cryptocurrencies and blockchains[END_REF]. This technology, used for information encryption, is pushing forward the demand for more efficient and secure pseudorandom number generators [START_REF] Menezes | Handbook of Applied Cryptography[END_REF] which, in the scope of chaos-based cryptography, were first introduced by Matthews in the 1990s [START_REF] Matthews | On the derivation of chaotic encryption algorithm[END_REF]. Contrarily to most algorithms that are used nowadays and based on a limited number of arithmetic or algebraic methods (like elliptic curves), networks of coupled chaotic maps offer quasi-infinite possibilities to generate parallel streams of pseudorandom numbers (PRN) at a rapid pace when they are executed on modern multicore processors. Chaotic maps are able to generate independent and secure pseudorandom sequences (used as information carriers or directly involved in the process of encryption/decryption [START_REF] Lozi | Noise-resisting ciphering based on a chaotic multi-stream pseudorandom number generator[END_REF]). However, the majority of well-known chaotic maps are not naturally suitable for encryption [START_REF] Li | Period extension and randomness enhancement using high-throughput reseeding-mixing PRNG[END_REF] and most of them do not exhibit even satisfactory properties for such a purpose.

In this chapter, we explore the novel idea of coupling a symmetric tent map with a logistic map, following several network topologies. We add a specific injection mechanism to capture the escaping orbits. In the goal of extending our results to industrial mathematics, we implement these networks on multicore machines and we test up to 100 trillion iterates of such mappings, in order to make sure that the obtained results are firmly grounded and able to be used in industrial contexts such as e-banking, e-purchasing, or the Internet of Things (IoT).

The chaotic maps, when used in the sterling way, could generate not only chaotic numbers, but also pseudorandom numbers as shown in [START_REF] Noura | Design of a fast and robust chaos-based cryptosystem for image encryption[END_REF] and as we show in this chapter with more sophisticated numerical experiments.

Various choices of PNR Generators (PRNGs) and crypto-algorithms are currently necessary to implement continuous, reliable security systems. We use a software approach because it is easy to change a cryptosystem to support protection, whereas p. 2

Chapter 4

replacing hardware used for True Random Number Generators would be costly and time-consuming. For instance, after the secure software protocol Wi-Fi Protected Access (WPA) was broken, it was simply updated and no expensive hardware had to be replaced.

It is a very challenging task to design CPRNGs (Chaotic Pseudo Random Number Generators) that are applicable to cryptography: numerous numerical tests must ensure that their properties are satisfactory. We mainly focus on two-to fivedimension maps, although upper dimensions can be very easily explored with modern multicore machines. Nevertheless, in four and five dimensions, the studied CRPNGs are efficient enough for cryptography.

In Sect. 4.2, we briefly recall the dawn and the maturity of researches on chaos. In Sect. 4.3, we explore two-dimensional topologies of networks of coupled chaotic maps. In Sect. [START_REF] Lozi | Noise-resisting ciphering based on a chaotic multi-stream pseudorandom number generator[END_REF].4, we study more thoroughly a mapping in higher dimensions (up to 5) far beyond the NIST tests which are limited to a few millions of iterates and which seem not robust enough for industrial applications, although they are routinely used worldwide. In order to check the portability of the computations on multicore architectures, we have implemented all our numerical experiments on several different multicore machines. We conclude this chapter in Sect. 4.5.

The Dawn and the Maturity of Researches on Chaos

The study of nonlinear dynamics is relatively recent with respect to the long historical development of early mathematics since the Egyptian and the Greek civilizations (and even before). The first alleged artifact of mankind's mathematical thinking goes back to the Upper Paleolithic era. Dating as far back as 22,000 years ago, the Ishango bone is a dark brown bone which happens to be the fibula of a baboon, with a sharp piece of quartz affixed to one end for engraving. It was first thought to be a tally stick, as it has a series of what has been interpreted as tally marks carved in three columns running the length of the tool [START_REF] Bogoshi | The oldest mathematical artifact[END_REF].

Twenty thousand years later, the Rhind Mathematical Papyrus is the best example of Egyptian mathematics. It dates back to around 1650 BC. Its author is the scribe Ahmes who indicated that he copied it from an earlier document dating from the 12th dynasty, around 1800 BC. It is a practical handbook, whose the first part consists of reference tables and a collection of 20 arithmetic and 20 algebraic problems and linear equations. Problem 32 for instance corresponds (in modern notation) to solving x + x 3 + x 4 = 2 for x [START_REF] Smith | History of Mathematics[END_REF]. Since those early times, mathematics have known great improvements, flourishing in many different fields such as geometry, algebra (both linked, thanks to the invention of Cartesian coordinates by René Descartes [START_REF] Descartes | Discours de la méthode[END_REF]), analysis, probability, number and set theory, and so on.

However, nonlinear problems are very difficult to handle, because, as shown by Galois' theory of algebraic equations which provides a connection between field theory and group theory, it is impossible to solve any polynomial equation p. 3 of degree equal or greater than 5 using only the usual algebraic operations (addition, subtraction, multiplication, division) and the application of radicals (square roots, cube roots, etc.) [START_REF] Galois | Mémoire sur les conditions de résolubilité des équations par radicaux (mémoire manuscrit de 1830)[END_REF].

The beginning of the study of nonlinear equation systems goes back to the original works of Gaston Julia and Pierre Fatou regarding to one-dimensional maps with a complex variable, nearly a century ago [START_REF] Julia | Mémoire sur l'itération des fonctions rationnelles[END_REF][START_REF] Fatou | Sur l'itération des fonctions transcendantes entières[END_REF]. Compared to thousands of years of mathematical development, a century is a very short period. In France, 30 years later, Igor Gumosky and Christian Mira began their mathematical researches with the help of a computer in 1958 [START_REF] Gumowski | Recurrence and Discrete Dynamics systems[END_REF]. They developed very elaborate studies of iterations. One of the best-known formulas they published is

x n+1 = f (x n ) + by n y n+1 = f (x n+1 ) -x n , with f (x) = ax + 2(1 -a) x 2 1 + x 2 (4.1)
which can be considered as a non-autonomous mapping from the plane R 2 onto itself that exhibits esthetic chaos. Surprisingly, slight variations of the parameter value lead to very different shapes of the attractor (Fig. 4.1).

In Ukraine, Alexander Sharkovsky found the intriguing Sharkovsky's order, giving the periods of periodic orbits of such nonlinear maps in 1962, although these results were only published in 1964 [START_REF] Sharkovskiȋ | Coexistence of cycles of a continuous map of the line into itself[END_REF]. In Japan the Hayashi' School (with disciples like Yoshisuke Ueda and Hiroshi Kawakami), a few years later, was motivated by applications to electric and electronic circuits. Ikeda proposed the Ikeda attractor [START_REF] Ikeda | Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system[END_REF][START_REF] Ikeda | Optical turbulence: chaotic behavior of transmitted light from a ring cavity[END_REF] which is a chaotic attractor for u ≥ 0.6 (Fig. 4.2).

x n+1 = 1 + u(x n cos t n -y n sin t n ) y n+1 = u(x n sin t n + y n cos t n ), with t n = 0.4 - 6 1 + x 2 n + y 2 n (4.2)
In 1983, Leon O. Chua invented a famous electronic circuit that generates chaos built with only two capacitors, one inductor and one nonlinear negative resistance [START_REF] Chua | The double scroll family[END_REF]. Since then, thousands of papers have been published on the general topic of chaos. However the pace of mathematics is slow, because any progress is based on strictly rigorous proof. Therefore numerous problems still remain unsolved. For example, the long-term dynamics of the Hénon map [START_REF] Hénon | Two-dimensional mapping with a strange attractor[END_REF], the first example of a strange attractor for mappings, remains unknown close to the classical parameter values from a strictly mathematical point of view, 40 years after its original publication.

Nevertheless, in spite of this lack of rigorous mathematical results, nowadays, engineers are actively working on applications of chaos for several purposes: global optimization, genetic algorithms, CPRNG, cryptography, and so on. They use nonlinear maps for practical applications without the need of sophisticated theorems. During the last 20 years, several chaotic image encryption methods have been proposed in the literature.

Dynamical systems which present a mixing behavior and that are highly sensitive to initial conditions are called chaotic. Small differences in initial conditions (such as those due to rounding errors in numerical computation) yield widely diverging outcomes for chaotic systems. This effect, popularly known as the butterfly effect, renders long-term predictions impossible in general [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF]. This happens even though these systems are deterministic, meaning that their future behavior is fully determined by their initial conditions, with no random elements involved. In other words, the deterministic nature of these systems does not make them predictable. Mastering the global properties of those dynamical systems is a challenging issue nowadays that we try to fix by exploring several network topologies of coupled maps.

In this chapter, after giving some prototypical examples of industrial applications of iterations of nonlinear maps, we focus on the exploration of topologies of coupled nonlinear maps that have a very rich potential of complex behavior. Very long computations on multicore machines are used, generating up to one hundred trillion iterates, in order to assess such topologies. We show the emergence of randomness from chaos and discuss the promising future of chaos theory for cryptographic security.

p. 5

Miscellaneous Network Topologies of Coupled Chaotic Maps

Tent-Logistic Entangled Map

In this section we consider only two 1-D maps: the logistic map

f µ (x) ≡ L µ (x) = 1 -µx 2 (4.3)
and the symmetric tent map

f µ (x) ≡ T µ (x) = 1 -µ|x| (4.4)
both associated to the dynamical system

x n+1 = f µ (x n ), ( 4.5)
where µ is a control parameter which impacts the chaotic degree. Both mappings are sending the one-dimensional interval [-1,1] onto itself.

Since the first study by R. May [START_REF] May | Stability and Complexity of Models Ecosystems[END_REF][START_REF] May | Biological populations with nonoverlapping generations: stable points, stable cycles, and chaos[END_REF] of the logistic map in the frame of nonlinear dynamical systems, both the logistic (4.3) and the symmetric tent map (4.4) have been fully explored with the aim to easily generate pseudorandom numbers [START_REF] Lozi | Giga-periodic orbits for weakly coupled tent and logistic discretized maps[END_REF].

However, the collapse of iterates of dynamical systems [START_REF] Yuan | Collapsing of chaos in one dimensional maps[END_REF] or at least the existence of very short periodic orbits, their non-constant invariant measure, and the easilyrecognized shape of the function in the phase space, could lead to avoid the use of such one-dimensional maps (logistic, baker, tent, etc.) or two-dimensional maps (Hénon, Standard, Belykh, etc.) as PRNGs (see [START_REF] Lozi | Can we trust in numerical computations of chaotic solutions of dynamical systems?[END_REF] for a survey). Yet, the very simple implementation as computer programs of chaotic dynamical systems led some authors to use them as a base for cryptosystems [25,[START_REF] Ariffin | Modified baptista type chaotic cryptosystem via matrix secret key[END_REF]. Even if the logistic and tent maps are topologically conjugates (i.e., they have similar topological properties: distribution, chaoticity, etc.), their numerical behavior differs drastically due to the structure of numbers in computer realization [START_REF] Lanford | Informal remarks on the orbit structure of discrete approximations to chaotic maps[END_REF].

As said above, both logistic and tent maps are never used in serious cryptography articles because they have weak security properties (collapsing effect) if applied alone. Thus, these maps are often used in modified form to construct CPRNGs [START_REF] Wong | A modified chaotic cryptographic method[END_REF][START_REF] Nejati | A realizable modified tent map for true random number generation[END_REF][START_REF] Lozi | Mathematical chaotic circuits: an efficient tool for shaping numerous architectures of mixed chaotic/pseudo random number generator[END_REF].

Recently, Lozi et al. proposed innovative methods in order to increase randomness properties of the tent and logistic maps over their coupling and sub-sampling [START_REF] Lozi | Emergence of randomness from chaos[END_REF][START_REF] Rojas | New alternate ring-coupled map for multirandom number generation[END_REF][START_REF] Garasym | Robust PRNG based on homogeneously distributed chaotic dynamics[END_REF]. Nowadays, hundreds of publications on industrial applications of chaosbased cryptography are available [START_REF] Jallaouli | Design and analyses of two stream ciphers based on chaotic coupling and multiplexing techniques[END_REF][START_REF] Garasym | Application of observer-based chaotic synchronization and identifiability to the original CSK model for secure information transmission[END_REF][START_REF] Farajallah | Fast and secure chaos-based cryptosystem for images[END_REF][START_REF] Taralova | Chaotic generator synthesis: dynamical and statistical analysis[END_REF].

In this chapter, we explore more thoroughly the original idea of combining features of tent (T µ ) and logistic (L µ ) maps to produce a new map with improved properties, through combination in several network topologies. This idea was recently introduced [START_REF] Garasym | How useful randomness for cryptography can emerge from multicore-implemented complex networks of chaotic maps?[END_REF]39] in order to improve previous CPRNGs. Looking at both Eqs. (4.3) and (4.4), it is possible to reverse the shape of the graph of the tent map T and to entangle it with the graph of the logistic map L. We obtain the combined map

f µ (x) ≡ T L µ (x) = µ|x| -µx 2 = µ(|x| -x 2 ) (4.6)
When used in more than one dimension, the T L µ map can be considered as a twovariable map

T L µ (x (i) , x ( j) ) = µ(|x (i) | -(x ( j) ) 2 ), i = j (4.7)
Moreover, we can combine again the T L µ map with T µ in various ways. If with choose, for instance, a network with a ring shape (Fig. 4.3). It is possible to define a mapping M µ, p : J p → J p where

J p = [-1, 1] p ⊂ R p : M µ, p ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ x (1) n x (2) n . . . x ( p) n ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ x (1) n+1 x (2) n+1 . . . x ( p) n+1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ T µ (x (1) n ) + T L µ (x (1) n , x (2) n ) T µ (x (2) n ) + T L µ (x (2) n , x (3) n ) . . . T µ (x ( p) n ) + T L µ (x ( p) n , x (1) n ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (4.8)
However, if used in this form, system (4.8) has unstable dynamics and iterated points x (1) n , x (2) n , . . . ,

x ( p)
n quickly spread out. Therefore, to solve the problem of keeping dynamics in the torus J p = [-1, 1] p ⊂ R p , the following injection mechanism has to be used in conjunction with (4.8) if (x (i) n+1 < -1) then add 2 if (x (i) n+1 > 1) then subtract 2

, i = 1, 2, . . . , p. (4.9)
p. 7 The T L µ function is a powerful tool to change dynamics. Used in conjunction with T µ , the map T L µ makes it possible to establish mutual influence between system components x (i) n in M µ, p . This multidimensional coupled mapping is interesting because it performs contraction and distance stretching between components, improving chaotic distribution.

The coupling of components has an excellent effect in achieving chaos, because they interact with global system dynamics, being a part of them. Component interaction has a global effect. In order to study this new mapping, we use a graphical approach, however other theoretical assessing functions are also involved.

Note that system (4.8) can be made more generic by introducing constants k i which generalize considered topologies. Let k = (k 1 , k 2 , . . . , k p ), we define

M k µ, p ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ x (1) n x (2) n . . . x ( p) n ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ x (1) n+1 x (2) n+1 . . . x ( p) n+1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ T µ (x (1) n ) + k 1 × T L µ (x (i) n , x ( j) n ), i, j = (1, 2) or (2, 1) T µ (x (2) n ) + k 2 × T L µ (x (i) n , x ( j) n ) i, j = (2, 3) or (3, 2) . . . T µ (x ( p) n ) + k p × T L µ (x (i) n , x ( j) n ) i, j = ( p, 1) or (1, p) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (4.10) System (4.10) is called alternate if k i = (-1) i or k i = (-1) i+1 , 1 ≤ i ≤ p, or non-alternate if k i = +1, or k i = -1. It
i j i' j' #1 +1 +1 1 2 1 2 #2 +1 -1 1 2 1 2 #3 -1 +1 1 2 1 2 #4 -1 -1 1 2 1 2 #5 +1 +1 2 1 2 1 #6 +1 -1 2 1 2 1 #7 -1 +1 2 1 2 1 #8 -1 -1 2 1 2 1 #9 +1 +1 1 2 2 1 #10 +1 -1 1 2 2 1 #11 -1 +1 1 2 2 1 #12 -1 -1 1 2 2 1 #13 +1 +1 2 1 1 2 #14 +1 -1 2 1 1 2 #15 -1 +1 2 1 1 2 #16 -1 -1 2 1 1 2

Two-Dimensional Network Topologies

We first consider the simplest coupling case, in which only two equations are coupled.

The first condition needed to obtain a multidimensional mapping, in the aim of building a new CPRNG, is to obtain excellent uniform distribution of the iterated points. The second condition is that the CPRNG must be assessed positively by the NIST tests [START_REF] Rukhin | Statistical test suite for random and pseudorandom number generators for cryptographic applications[END_REF]. In [START_REF] Garasym | How useful randomness for cryptography can emerge from multicore-implemented complex networks of chaotic maps?[END_REF]39] this two-dimensional case is studied in detail. Using a bifurcation diagram and computation of Lyapunov exponents, it is shown that the best value for the parameter is µ = 2. Therefore, in the rest of this chapter we use this parameter value and we only briefly recall the results found with this value in both of those articles. The general form of

M k 2,2 is then M k 2,2 x (1) n x (2) n = x (1) n+1 x (2) n+1 = T 2 (x (1) n ) + k 1 × T L 2 (x (i) n , x ( j) n ) T 2 (x (2) n ) + k 2 × T L 2 (x (i ) n , x ( j ) n ) (4.11)
with i, j, i , j = 1 or 2, i = j, and i = j . Considering this general form, it is possible to define 16 different maps (Table 4.1). Among this set of maps, we study case #3 and case #13. The map of case #3 is called Single-Coupled alternate due to the shape of the corresponding network and denoted T T L SC 2 ,

p. 9

T T L SC 2 = ⎧ ⎨ ⎩ x (1) n+1 = 1 -2|x (1) n | -2(|x (1) n | -(x (2) n ) 2 ) = T 2 (x (1) n ) -T L 2 ((x (1) n ), (x (2) n )) x (2) n+1 = 1 -2|x (2) n | + 2(|x (1) n | -(x (2) n ) 2 ) = T 2 (x (2) n ) + T L 2 ((x (1) n ), (x (2) n ))
(4.12) and case #13 is called Ring-Coupled non-alternate and denoted T T L RC 2 ,

T T L RC 2 = ⎧ ⎨ ⎩ x (1) n+1 = 1 -2|x (1) n | + 2(|x (2) n | -(x (1) n ) 2 ) = T 2 (x (1) n ) + T L 2 ((x (2) n ), (x (1) n )) x (2) n+1 = 1 -2|x (2) n | + 2(|x (1) n | -(x (2) n ) 2 ) = T 2 (x (2) n ) + T L 2 ((x (1) n ), (x (2) n )) (4.13)
Both systems were selected because they have balanced contraction and stretching processes between components. They allow achieving uniform distribution of the chaotic dynamics. Equations (4.12) and (4.13) are used, of course, in conjunction with injection mechanism (4.9). The largest torus where points mapped by (4.12) and (4.13) are sent is [-2, 2] 2 . The confinement from torus [-2, 2] 2 to torus [-1, 1] 2 of the dynamics obtained by this mechanism is shown in Figs. 4.5 and 4.6: dynamics cross from the negative region (in blue) to the positive one, and conversely to the negative region, if the points stand in the positive regions (in red). Through this operation, the system's dynamics are trapped inside [-1, 1] 2 . In addition, after this operation is done, the resulting system exhibits more complex dynamics with additional nonlinearity, which is advantageous for chaotic encryption (since it improves security).

A careful distribution analysis of both T T L SC 2 and T T L RC 2 has been performed using approximated invariant measures. x (1) n ≡ x (1) n -2; if x (1) n < -1 then x (1) n ≡ x (1) n + 2 (from [START_REF] Garasym | How useful randomness for cryptography can emerge from multicore-implemented complex networks of chaotic maps?[END_REF])

p. 10

n ≡ x (2) n -2; if x (2) n < -1 then x (2) n ≡ x (2)
n + 2 (from [START_REF] Garasym | How useful randomness for cryptography can emerge from multicore-implemented complex networks of chaotic maps?[END_REF])

Approximated Invariant Measures

We recall in this section the definition of approximated invariant measures which are important tools for assessing the uniform distribution of iterates. We have previously introduced them for the first studies of the weakly coupled symmetric tent map [START_REF] Lozi | Giga-periodic orbits for weakly coupled tent and logistic discretized maps[END_REF]. We first define an approximation P M,N (x) of the invariant measure, also called the probability distribution function linked to the one-dimensional map f (Eq. (4.5)) when computed with floating numbers (or numbers in double precision). To this goal, we consider a regular partition of M small intervals (boxes) r i of J = [-1, 1] defined by

s i = -1 + 2i M , i = 0, M, (4.14)
r i = [s i , s i+1 [, i = 0, M -2, (4.15) r M-1 = [s M-1 , 1], (4.16
)

J = M-1 0 r i . (4.17)
The length of each box r i is equal to

s i+1 -s i = 2 M (4.18)
All iterates f (n) (x) belonging to these boxes are collected (after a transient regime of Q iterations decided a priori, i.e., the first Q iterates are discarded). Once the p. 11 computation of N + Q iterates is completed, the relative number of iterates with respect to N /M in each box r i represents the value P N (s i ). The approximated P N (x) defined in this article is therefore a step function, with M steps. Since M may vary, we define

P M,N (s i ) = 1 2 M N (#r i ) (4.19)
where #r i is the number of iterates belonging to the interval r i and the constant 1/2 allows the normalisation of P M,N (x) on the interval J .

P M,N (x) = P M,N (s i ), ∀ x ∈ r i (4.20)
In the case of p-coupled maps, we are more interested by the distribution of each component x (1) , x (2) , . . . , x ( p) of the vector

X = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
x (1) x (2) . . .

x ( p) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
rather than by the distribution of the variable X itself in J p . We then consider the approximated probability distribution function P M,N (x ( j) ) associated to one component of X . In this chapter, we use either N disc for M or N iter for N , depending on which is more explicit. The discrepancies E 1 (in norm L 1 ), E 2 (in norm L 2 ), and E ∞ (in norm L ∞ ) between P N disc ,N iter (x) and the Lebesgue measure, which is the invariant measure associated to the symmetric tent map, are defined by

E 1,N disc ,Niter (x) = P N disc ,Niter (x) -0.5 L 1 (4.21) E 2,N disc ,Niter (x) = P N disc ,Niter (x) -0.5 L 2 (4.22) E ∞,N disc ,Niter (x) = P N disc ,Niter (x) -0.5 L ∞ (4.23)
In the same way, an approximation of the correlation distribution function C M,N (x, y) is obtained by numerically building a regular partition of M 2 small squares (boxes) of J 2 , embedded in the phase subspace (x l , x m )

s i = -1 + 2i M , t j = -1 + 2 j M , i, j = 0, M (4.24) r i, j = [s i , s i+1 [×[t j , t j+1 [, i, j = 0, M -2 (4.25) r M-1, j = [s M-1 , 1] × [t j , t j+1 ], j = 0, M -2 (4.26) r i,M-1 = [s i , s i+1 [×[t M-1 , 1], j = 0, M -2 (4.27) r M-1,M-1 = [s M-1 , 1] × [t M-1 , 1] (4.28)
p. 12

The measure of the area of each box is

(s i+1 -s i ).(t i+1 -t i ) = 2 M 2 (4.29)
Once N + Q iterated points (x l n , x m n ) belonging to these boxes are collected, the relative number of iterates with respect to N /M 2 in each box r i, j represents the value C N (s i , t j ). The approximated probability distribution function C N (x, y) defined here is then a two-dimensional step function, with M 2 steps. Since M can take several values in the next sections, we define

C M,N (s i , t j ) = 1 4 M 2 N (#r i, j ) (4.30)
where #r i, j is the number of iterates belonging to the square r i, j and the constant 1/4 allows the normalisation of C M,N (x, y) on the square J 2 .

C M,N (x, y) = C M,N (s i , t j ) ∀(x, y) ∈ r i, j (4.31)
The discrepancies y) and the uniform distribution on the square are defined by

E C 1 (in norm L 1 ), E C 2 (in norm L 2 ) and E C ∞ (in norm L ∞ ) between C N disc ,N iter (x,
E C 1 ,N disc ,N iter (x, y) = C N disc ,N iter (x, y) -0.25 L 1 (4.32) E C 2 ,N disc ,N iter (x, y) = C N disc ,N iter (x, y) -0.25 L 2 (4.33) E C ∞ ,N disc ,N iter (x, y) = C N disc ,N iter (x, y) -0.25 L ∞ (4.34)
Finally, let AC N disc ,N iter be the autocorrelation distribution function which is the correlation function C N disc ,N iter of (4.31), defined in the delay space (x (i) n , x (i) n+1 ) instead of the phase (x l , x m ) space. We define in the same manner than (4.32), (4.33), and (4.34)

E C 1 ,N disc ,N iter (x, y), E C 2 ,N disc ,N iter (x, y), and E C ∞ ,N disc ,N iter (x, y).

Study of Randomness of TTL SC

2 and TTL RC

, and Other Topologies

Using numerical computations, we assess the randomness properties of the two-dimensional maps T T L SC 2 and T T L RC 2 . If all requirements 1-8 of Fig. 4.7 are verified, the dynamical systems associated to those maps can be considered as pseudorandom and their application to cryptosystems is possible.

p. 13 Whenever one among the eight criteria is not satisfied for a given map, one cannot consider that the associated dynamical system is a good CPRNG candidate. As said above, when µ = 2, the Lyapunov exponents of both considered maps are positive.

In the phase space, we plot the iterates in the system of coordinates x (1) n versus x (2) n in order to analyze the density of the points' distribution. Based on such an analysis, it is possible to assess the complexity of the behavior of dynamics, noticing any weakness or inferring on the nature of randomness. We also use the approximate invariant measures to assess more precisely the distribution of iterates.

The graphs of the attractor in phase space for the T T L RC 2 non-alternate (Fig. 4.8) and T T L SC 2 alternate (Fig. 4.9) maps are different. The T T L SC 2 map has wellscattered points in the whole pattern, but there are some more "concentrated" regions forming curves on the graph. Instead, the map T T L RC 2 has good repartition. Some other numerical results we do not report in this chapter show that even if those maps have good random properties, it is possible to improve mapping randomness by modifying slightly network topologies.

p. 14

T T L SC 2 (x (1) n , x (2) n ) = x (1) n+1 = 1 + 2(x (2) n ) 2 -4|x (1) n | x (2) n+1 = 1 -2(x (2) n ) 2 + 2(|x (1) n | -|x (2) n |) (4.35)
In [START_REF] Garasym | How useful randomness for cryptography can emerge from multicore-implemented complex networks of chaotic maps?[END_REF], it is shown that if the impact of component x (1) n is reduced, randomness is improved. Hence, the following MT T L SC 2 map is introduced

MT T L SC 2 (x (1) n , x (2) n ) = x (1) n+1 = 1 + 2(x (2) n ) 2 -2|x (1) n | x (2) n+1 = 1 -2(x (2) n ) 2 + 2(|x (1) n | -|x (2) n |) (4.36)
and the injection mechanism (4.9) is used as well, but it is restricted to three phases: (1) n+1 > 1) then subtract 2 if (x (2) n+1 < -1) then add 2 if (x (2) n+1 > 1) then subtract 2 (4.37)

⎧ ⎪ ⎨ ⎪ ⎩ if (x
This injection mechanism allows the regions containing iterates to match excellently (Fig. 4.10).

The change of topology leading to MT T L SC 2 greatly improves the density of iterates in the phase space (Fig. 4.11) where 10 9 points are plotted. The point distribution of iterates in phase delay for the variable x (2) is quite good as well (Fig. 4.12). On both pictures, a grid of 200 × 200 boxes is generated to use the box counting method defined in Sect. 4.3.3. Moreover, the largest Lyapunov exponent is equal to 0.5905, indicating a strong chaotic behavior.

p. 15 2 alternative map, on the (x (1) , x (2) ) plane (from [START_REF] Garasym | How useful randomness for cryptography can emerge from multicore-implemented complex networks of chaotic maps?[END_REF]) However, regarding the phase delay for the variable x (1) , results are not satisfactory. We have plotted in Fig. 4.13 10 9 iterates of MT T L SC 2 in the delay plane, and in Fig. 4.14 the same iterates using the counting box method.

When such a great number of iterates is computed, one has to be cautious with raw graphical methods because irregularities of the density repartition are masked due to the huge number of plotted points. Therefore, these figures highlight the necessity of using the tools we have defined in Sect. 4.3.3.

Nevertheless, NIST tests were used to check randomness properties of MT T L SC 2 . Since they only require binary sequences, we generated 4 × 10 6 iterates whose 5 × 10 5 first ones were cut off. The rest of the sequence was converted to binary form according to the IEEE-754 standard (32-bit single-precision floating point).

p. 16

Fig. 4.12 Approximate density function of the MT T L SC

2 alternative map, on the (x (1) n , x (1) n+1 ) plane (from [START_REF] Garasym | How useful randomness for cryptography can emerge from multicore-implemented complex networks of chaotic maps?[END_REF]) As said in the introduction, networks of coupled chaotic maps offer quasi-infinite possibilities to generate parallel streams of pseudorandom numbers. For example, in [39], the following modification of MT T L SC 2 is also studied and shows good randomness properties

N T T L SC 2 (x (1) n , x (2) n ) = ⎧ ⎪ ⎨ ⎪ ⎩ x (1) n+1 = 1 -2|x (2) n | = T 2 (x (2) n ) x (2) n+1 = 1 -(2x (2) n ) 2 -2(|x (2) n | -|x (1) n |) = L 2 (x (2) n ) + T 2 (x (2) n ) -T 2 (x (1) n ) (4.38)
p. 17

Mapping in Higher Dimension

Higher dimensional systems make it possible to achieve better randomness and uniform point distribution, because more perturbations and nonlinear mixing are involved. In this section, we focus on a particular realization of the M k µ, p map (4.10) from dimension two to dimension five. (2) (from [START_REF] Garasym | How useful randomness for cryptography can emerge from multicore-implemented complex networks of chaotic maps?[END_REF])

p. 18

Usually, three or four dimensions are complex enough to create robust random sequences as we show here. Thus, it is advantageous if the system can increase its dimension. Since the MT T L SC 2 alternative map cannot be nested in higher dimensions, we describe how to improve randomness and to obtain the best distribution of points, and how to produce more complex dynamics than the T T L SC 2 (x (2) , x (1) ) alternative map in dimension greater than 2. Let

T T L RC, pD 2 = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ x (1) n+1 = 1 -2|x (1) n | + 2(|x (2) n | -(x (1) n ) 2 ) x (2) n+1 = 1 -2|x (2) n | + 2(|x (3) n | -(x (2) n ) 2 ) . . . x ( p) n+1 = 1 -2|x ( p) n | + 2(|x (1) n | -(x ( p) n ) 2 ) (4.39)
be this realization. We show in Figs. 4.17 and 4.18 successful NIST tests for T T L RC, pD 2 in 3-D and 4-D, for the variable x (1) .

Fig. 4.17 N I ST test for T T L RC,3D

2 for x (1) (from [START_REF] Garasym | How useful randomness for cryptography can emerge from multicore-implemented complex networks of chaotic maps?[END_REF]) p. 19 2 for x (1) (from [START_REF] Garasym | How useful randomness for cryptography can emerge from multicore-implemented complex networks of chaotic maps?[END_REF])

Numerical Experiments

All NIST tests for dimensions three to five for every variable are successful, showing that these realizations in 3-D up to 5-D are good CPRNGs. In addition to those tests, we study the mapping more thoroughly, far beyond the NIST tests which are limited to a few million iterates and which seem not robust enough for industrial mathematics, although they are routinely used worldwide.

In order to check the portability of the computations on multicore architectures, we have implemented all our numerical experiments on several different multicore machines.

Checking the Uniform Repartition of Iterated Points

We first compute the discrepancies E 1 (in norm L 1 ), E 2 (in norm L 2 )n and E ∞ (in norm E ∞ ) between P N disc ,N iter (x) and the Lebesgue measure which is the uniform measure on the interval J = [-1, 1]. We set M = N iter = 200, and vary the number N iter of iterated points in the range 10 4 to 10 14 . From our knowledge, this article is the first one that checks such a huge number of iterates (in conjunction with [39]). We compare E 1,200,N iter (x (1) ) for T T L RC, pD 2 with p = 2 to 5 (Table 4.2, Fig. 4

.19).

As shown in Fig. 4.19, E 1,200,N iter (x (1) ) decreases steadily when N iter increases. However, the decreasing process is promptly (with respect to N iter ) bounded below for p = 2. This is also the case for other values of p, however, the boundary decreases with p, therefore showing better randomness properties for higher dimensional mappings.

Table 4.3 compares x (1) , x (2) , …,x ( p) for T T L RC,5D

2

, for different values of N iter . It is obvious that the same quality of randomness is obtained for each one of them, contrarily to the results obtained for MT T L SC 2 .

p. 20 with p = 2 to 5, with respect to N iter (horizontal axis, logarithmic value) Table 4.3 E 1,200,N iter (x (i) ) for T T L RC,5D 2 for i = 1 to 5

N iter

x (1) x (2) x (3) x (4) x (5) 0.000160547 0.000159192 0.000160014 0.000159213 0.000159159 10 13 5.04473e-05 5.03574e-05 5.05868e-05 5.04694e-05 5.01681e-05 10 14 1.59929e-05 1.60291e-05 1.59282e-05 1.59832e-05 1.60775e-05 p. 21 (1) ), E 2,200,N iter (x (1) ), and E ∞,200,N iter (x (1) ) for T T L RC,5D 0.000160547 0.000201102 0.0008602 10 13 5.04473e-05 6.32233e-05 0.00026894 10 14 1.59929e-05 2.00533e-05 9.89792e-05 Fig. 4. [START_REF] May | Stability and Complexity of Models Ecosystems[END_REF] Comparison between E 1,200,N iter (x (1) ), E 2,200,N iter (x (1) ), and E ∞,200,N iter (x (1) ) (vertical axis) for T T L RC,5D 2 with respect to N iter (horizontal axis, logarithmic value)

The comparisons between E 1,200,N iter (x (1) ), E 2,200,N iter (x (1) ), and E ∞,N iter (x (1) ) for T T L RC,5D 2 in Table 4.4 and Fig. 4.20 show that E 1,200,N iter (x (1) ) < E 2,200,N iter (x (1) ) < E ∞,N iter (x (1) ) (4.40)

for every value of N iter .

Autocorrelation Study in the Delay Space

In this section, we assess autocorrelation errors

E AC 1 ,N disc ,N iter (x, y), E AC 2 ,N disc ,N iter (x,
n+2 ), and E AC1,200,N iter (x (1) n , x (1) n+3 (1) n , x (1) n+1 ), E AC1,200,N iter (x (1) n , x (1) n+2 ), and ments for M = 20 to 20, 000, however, in this chapter, we only present the results for M = 200. We first compare E AC 1 ,200,N iter (x (1) n , x (1) n+1 ) with E AC 1 ,200,N iter (x (1) n , x (1) n+2 ) and E AC 1 ,200,N iter (x (1) n , x (1) n+3 ) for T T L RC, pD 2 when the dimension of the system is within the range p = 2 to 5 (Tables 4.5, 4.6, 4.7 and 4.8). It is possible to see that better randomness properties are obtained for higher dimensional mappings.

) for T T L RC,2D 2 N iter (x (1) n , x (1) n+1 ) (x (1) n , x (1) n+2 ) (x (1) n , x (1)
E AC1,200,N iter (x (1) n , x (1) n+3 ) for T T L RC,3D 2 N iter (x (1) n , x (1) n+1 ) (x (1) n , x (1) n+2 ) (x (1) n , x (1) n+3 ) 10 4 1
The comparison between E AC 1 ,200,N iter (x (1) n , x (1) n+1 ), E AC 2 ,200,N iter (x (1) n , x (1) n+1 ), and E AC ∞ ,200,N iter (x (1) n , x (1) n+1 ) for T T L RC,5D 2 in Table 4.9 shows that numerically p. 23

n+2 ), and E AC1,200,N iter (x (1) n , x (1) n+3 ) for T T L RC,4D

2 0.000160547 0.000159144 0.000159246 10 13 5.0394e-05 10 14 1.59929e-05

N iter (x (1) n , x (1) n+1 ) (x (1) n , x (1) n+2 ) (x (1) n , x (1)
E AC 1 ,200,N iter (x (1) n , x (1) n+1 ) < E AC 2 ,200,N iter (x (1) n , x (1) n+1 ) < E AC ∞ ,200,N iter (x (1) n , x (1) n+1 ) (4.41) Equation (4.41) is not only valid for M = 200, but also for other values of M and every component of X .

In order to illustrate the numerical results displayed in these tables, we plot in Fig. 4.21 the repartition of iterates of T T L RC,5D 2 in the delay plane (x (1) n , x (1) n+1 ), using the box counting method. On a grid of 200 × 200 boxes (N iter = M = 200), p. 24 (1) n , x (1) n+1 ) for T T L RC,5D 0.000160547 0.000201102 0.0008602 10 13 5.0394e-05 6.31756e-05 0.000280168 10 14 1.59929e-05 2.00533e-05 9.89792e-05 we have generated 10 6 points. The horizontal axis is x (1) n , and the vertical axis is x (1) n+1 . In order to check very carefully the repartition of the iterates of T T L RC,5D 2 , we have also plotted the repartition in the delay planes (x (1) n , x (1) n+2 ), (x (1) n , x (1) n+3 ), and (x (1) n , x (1) n+4 ) (Figs. 4.22, 4.23, and 4.24). This repartition is uniform everywhere as shown also in Table 4.8.

We find the same regularity for every component x (2) , x (3) , x (4) , and x (5) , as shown in Figs. 4.25, 4.26, 4.27, 4.28, and in

Autocorrelation Study in the Phase Space

Finally, in this section, we assess the autocorrelation errors E C 1 ,N disc ,N iter (x, y), E C 2 ,N disc ,N iter (x, y), and E C ∞ ,N disc ,N iter (x, y), defined by Eqs. (4.32), (4.33), and (4.34), in the phase space. We checked all combinations of the components. Due to space limitations, we only provide part of the numerical computations we have performed to carefully check the randomness of T T L RC, pD 2 for p = 2, 5 and i = 1, p. Like in the previous section, we only provide the results for M = 200. We first compare E C 1 ,200,N iter (x (1) n , x (2) n ), E C 2 ,200,N iter (x (1) n , x (2) n ), and E C ∞ ,200,N iter (x (1) n , x (2) n ) (Table 4.11), and our other results verified that p. 26 (2) n , x (2) n+1 ) of T T L RC,5D 2 ; box counting method, 10 6 points are generated on a grid of 200 × 200 boxes, the horizontal axis is x (2) n , and the vertical axis is x (1) n ,

(2) n+1 E C 1 ,N disc ,N iter (x
x (2) n ) < E C 2 ,N disc ,N iter (x (1) n , x (2) n ) < E C ∞ ,N disc ,N iter (x (1) n , x (2) n ) (4.42)
We have also assessed the autocorrelation errors

E C 1 ,N disc ,N iter (x (i) n , x ( j)
n ) for i, j = 1, 5, i = j, and various values of the number of iterates for T T L RC,5D 2 (Table 4.12). We have performed the same experiments for E C 1 ,N disc ,N iter (x (1) n , x (2) n ) for p = 1, 5 (Table 4.13).

p. 27 Our numerical experiments all show a similar trend: T T L RC, pD 2 is a good candidate for a CPRNG, and the randomness performance of such mappings increases in higher dimensions.

Checking the Influence of Discretization in Computation of Approximated Invariant Measures

In order to verify that the computations we have performed using the discretization M = N disc = 200 of the phase space and the delay space in the numerical experip. 28

(i) n , x (i) n+1 ), E AC1,200,N iter (x (i) n , x (i) n+2 ), and E AC1,200,N iter (x (i) n , x (i) n+3 ) for T T L RC,5D 2 for i = 1 to 5 N iter i (x (i) n , x (i) n+1 ) (x (i) n , x (i) n+2 ) (x (i) n , x (i) n+3 ) 10 4 1

Computation Time of PRNs

The numerical experiments performed in this section have involved several multicore machines. We show in Table 4.15 different computation times (in seconds) for the generation of N iter PRNs for T T L RC, pD 2 with p = 2 to 5, and various values of the number of iterates (N iter ). The machine used is a laptop computer with a Core i7 4980HQ processor with eight logical cores. Table 4.16 shows the computation time of only one PRN in the same experiment. Time is expressed in 10 -10 s.

p. 30 These results show that the pace of computation is very high. When T T L RC,5D 2 is the mapping tested, and the machine used is a laptop computer with a Core i7 4980HQ processor with 8 logical cores, computing 10 11 iterates with five parallel streams of PRNs leads to around 2 billion PRNs being produced per second. Since these PRNs are computed in the standard double precision format, it is possible to extract from each 50 random bits (the size of the mantissa being 52 bits for a double precision floating-point number in standard IEEE-754). Therefore, T T L RC,5D

Conclusion

In this chapter, we thoroughly explored the novel idea of combining features of a tent map (T µ ) and a logistic map (L µ ) to produce a new map with improved properties, through combination in several network topologies. This idea was recently introduced [START_REF] Garasym | How useful randomness for cryptography can emerge from multicore-implemented complex networks of chaotic maps?[END_REF]39] in order to improve previous CPRNGs. We have summarized the previously explored topologies in dimension two. We have presented new results of numerical experiments in higher dimensions (up to five) for the mapping T T L RC, pD 2 on multicore machines and shown that T T L RC,5D 2 is a very good CPRNG which is fit for industrial applications. The pace of generation of random bits can be incredibly high (up to 200 billion random bits per second).

Fig. 4 . 1 Fig. 4 . 2

 4142 Fig. 4.1 Gumowski-Mira attractor for parameter values a = 0.92768 and a = 0.93333

Fig. 4 . 3

 43 Fig. 4.3 Auto and ring-coupling of the T L µ and T µ maps (from [38])

Fig. 4 . 4

 44 Fig. 4.4 Return mechanism from the [-2, 2] p torus to [-1, 1] p (from [38])

Fig. 4 . 5

 45 Fig. 4.5 Injection mechanism of the iterates from torus [-2, 2] 2 to torus [-1, 1] 2 . If x (1) n > 1 thenx(1) n ≡ x(1) n -2; if x(1)

Fig. 4 . 6

 46 Fig. 4.6 If x (2) n > 1 then x (2)

Fig. 4 . 7 Fig. 4 . 8

 4748 Fig. 4.7 The main criteria for assessing CPRNG (from [34])

Fig. 4 . 9

 49 Fig. 4.9 Phase space behavior of T T L SC 2 alternative (4.18), plot of 20, 000 points

Fig. 4 . 10

 410 Fig. 4.10 Injection mechanism (4.21) of the MT T L SC 2 alternative map (From [38])

Fig. 4 .

 4 Fig. 4.13 Plot of one billion iterates of MT T L SC 2 in the delay plane

Fig. 4 .Fig. 4 .

 44 Fig. 4.14 Plot of one billion iterates of MT T L SC 2 using the counting box method

Fig. 4 .

 4 Fig. 4.16 N I ST tests for the variable x(2) (from[START_REF] Garasym | How useful randomness for cryptography can emerge from multicore-implemented complex networks of chaotic maps?[END_REF])

Fig. 4 .

 4 Fig. 4.18 N I ST test for T T L RC,4D

Fig. 4 .

 4 Fig. 4.21 Repartition of iterates in the delay plane (x (1) n , x (1) n+1 ) of T T L RC,5D2

Fig. 4 . 2 , 21 Fig. 4 .

 42214 Fig. 4.[START_REF] Lozi | Giga-periodic orbits for weakly coupled tent and logistic discretized maps[END_REF] Repartition of iterates in the delay plane (x(1) n , x(1) n+2 ) of T T L RC,5D 2 , as in Fig.4.21

Fig. 4 . 2 , 21 Fig. 4 .

 42214 Fig.4.[START_REF] Lozi | Can we trust in numerical computations of chaotic solutions of dynamical systems?[END_REF] Repartition of iterates in the delay plane (x(1) n , x(1) n+4 ) of T T L RC,5D 2 , as in Fig.4.21

Fig. 4 . 2 ,Fig. 4 .

 424 Fig. 4.26 Repartition of iterates in the delay plane (x (3) n , x (3) n+1 ) of T T L RC,5D 2 , as in Fig. 4.25

Fig. 4 . 2 ,

 42 Fig. 4.[START_REF] Wong | A modified chaotic cryptographic method[END_REF] Repartition of iterates in the delay plane (x(5) n , x(5) n+1 ) of T T L RC,5D 2 , as in Fig.4.25

Fig. 4 . 29 2 ,

 4292 Fig. 4.29 Comparison between E C1,N disc ,N iter (x (1) n , y (2) n ), for T T L RC,4D 2 , M = N disc = 20, 200, 2000, 20, 000, and various values of the number of iterates

Table 4 . 1

 41 can be a mix of alternate and non-alternate if k i = +1 or -1 randomly. The sixteen maps defined by Eq.(4.11)

	p. 8

Table 4 .

 4 2 E 1,200,N iter (x(1) ) for T T L

			RC, pD 2	with p = 2 to 5
	N iter	p = 2	p = 3	p = 4	p = 5
	10 4	1.5631	1.5553	1.5587	1.5574
	10 5	0.55475	0.5166	0.51315	0.5154
	10 6	0.269016	0.159306	0.158548	0.158436
	10 7	0.224189	0.050509	0.0501934	0.0505558
	10 8	0.219427	0.0164173	0.0159175	0.0160018
	10 9	0.218957	0.00640196	0.00505021	0.00509754
	10 10	0.218912	0.00420266	0.00160505	0.00160396
	10 11	0.218913	0.00392507	0.000513833	0.000505591
	10 12	0.218913	0.00389001	0.000189371	0.000160547
	10 13	0.218914	0.00388778	0.000112764	5.04473e-05
	10 14	0.218914	0.003887	0.000101139	1.59929e-05

Fig. 4.19 Graph of E 1,200,N iter (x (1) ) for T T L RC, pD 2

Table 4 . 4

 44 Comparison between E 1,200,N iter (x

Table 4 . 5

 45 Comparison between E AC1,200,N iter (x

	(1) n , x	(1) n+1 ), E AC1,200,N iter (x	(1)

n , x

Table 4 . 6

 46 Comparison between E AC1,200,N iter (x

	n+3 )

Table 4 . 7

 47 Comparison between E AC1,200,N iter (x

	(1) n , x	(1) n+1 ), E AC1,200,N iter (x	(1)

n , x

Table 4 . 9

 49 Comparison between E AC1,200,N iter (x

	(1) n , x	(1) n+1 ), E AC2,200,N iter (x	(1) n , x	(1) n+1 ), and

Table 4 .

 4 10.

	p. 25

Table 4 .

 4 10 Comparison between E AC1,200,N iter (x

Table 4 .

 4 11 Comparison between E AC1,200,N iter (x

	(1) n , x	(2) n ), E AC2,200,N iter (x	(1) n , x	(2) n ), and

Table 4 .

 4 [START_REF] Fatou | Sur l'itération des fonctions transcendantes entières[END_REF] Comparison between E C1,200,N iter (x (i) n , x

				( j) n ), for i, j = 1 to 5, i = j, and for various
	values of number of iterates for T T L RC,5D 2
	N iter	10 6	10 8	10 10	10 12	10 14
	x(1), x(2)	0.158058	0.0160114	0.0015927	0.000158795 1.60489e-05
	x(1), x(3)	0.158956	0.0159261	0.00159456	0.000159326 1.73852e-05
	x(1), x(4)	0.15943	0.0160321	0.00160091	0.000160038 1.74599e-05
	x(1), x(5)	0.159074	0.0158962	0.00160204	0.000159048 1.59133e-05
	x(2), x(3)	0.15825	0.0159754	0.00159442	0.000160659 1.60419e-05
	x(2), x(4)	0.159248	0.0159668	0.00159961	0.000160313 1.73507e-05
	x(2), x(5)	0.15889	0.0160116	0.0015934	0.000160462 1.73496e-05
	x(3), x(4)	0.159136	0.0158826	0.00158123	0.000158758 1.59451e-05
	x(3), x(5)	0.159216	0.0159341	0.00161268	0.000159079 1.75013e-05
	x(4), x(5)	0.158918	0.0160516	0.0016008	0.000159907 1.59445e-05

Table 4 .

 4 [START_REF] Gumowski | Recurrence and Discrete Dynamics systems[END_REF] Comparison between E C1,200,N iter (x

			(i) n , x	( j) n ), for T T L	RC, pD 2	for p = 2, . . . , 5, and
	various values of the number of iterates
	N iter	p = 2	p = 3	p = 4		p = 5
	10 4	1.5624	1.5568	1.55725		1.55915
	10 5	0.57955	0.5163	0.51083		0.514
	10 6	0.330084	0.160282	0.158256	0.158058
	10 7	0.294918	0.0509584	0.0504002	0.0505508
	10 8	0.291428	0.0176344	0.0157924	0.0160114
	10 9	0.291012	0.00911485	0.00506758	0.00507915
	10 10	0.291025	0.00783204	0.00159046	0.0015927
	10 11	0.291033	0.00771201	0.000521561	0.000506086
	10 12	0.291036	0.00769998	0.000209109	0.000158795
	10 13		0.00769867	0.000150031	5.03666e-05
	10 14		0.00769874	0.000144162	1.60489e-05

Table 4 .

 4 14 Comparison E C1,N disc ,N iter (x

	(1) n , x	(2) n ),for T T L RC,4D 2	M = N disc =

Table 4 .

 4 15 Comparison of computation times (in second) for the generation of N iter PRNs for T T L RC, pD 2 with p = 2 to 5, and various values of N iter iterates

	N iter	p = 2	p = 3	p = 4	p = 5
	10 4	0.000146	0.000216	0.000161	0.000142
	10 5	0.000216	0.000277	0.000262	0.000339
	10 6	0.001176	0.002403	0.001681	0.002467
	10 7	0.011006	0.016195	0.018968	0.022351
	10 8	0.113093	0.161776	0.166701	0.227638
	10 9	1.09998	1.58949	1.60441	2.29003
	10 10	11.4901	18.0142	18.537	26.1946
	10 11	123.765	183.563	185.449	257.244

Table 4 .

 4 [START_REF] Ikeda | Optical turbulence: chaotic behavior of transmitted light from a ring cavity[END_REF] Comparison of computation times (in 10 -10 s) for the generation of only one PRN for T T L

	RC, pD 2	with p = 2 to 5, and various values of the number of iterates
	N iter	p	p = 3	p 4	p = 5
	10 4	73.0	72.0	40.25	28.4
	10 5	10.8	9.233	6.55	6.78
	10 6	5.88	8.01	4.2025	4.934
	10 7	5.503	5.39833	4.742	4.702
	10 8	5.65465	4.0444	4.16753	4.55276
	10 9	5.4999	5.2983	4.01103	4.58006
	10 10	5.74505	4.50335	4.63425	5.23892
	10 11	6.18825	6.11877	4.63622	5.14488

can produce 100 billion random bits per second, an incredible pace! With a machine with 4 Intel Xeon E7-4870 processors having a total of 80 logical cores, the computation is twice as fast, producing 2 × 10 11 random bits per second.p.


Nicolas Drouhin

Theoretical considerations on the retirement consumption puzzle and the optimal age of retirement

Keywords: C61 D91 J26 life cycle theory of consumption and saving; optimal retirement, retirement consumption puzzle, discontinuous optimal control

 principle of optimality, it provides a very general and parsimonious formula for determining the optimal age of retirement taking into account the possible discontinuity of the optimal consumption profile at the age of retirement.

Introduction

In this article, I build a model that address at the same time the retirement consumption puzzle and and the optimal age of retirement.

Since Hamermesh (1984a) many empirical studies document a drop in consumption at retirement, the retirement consumption puzzle [START_REF] Banks | Is there a retirement-savings puzzle?[END_REF][START_REF] Bernheim | What accounts for the variation in retirement wealth among us households?[END_REF]Battistin et al., 2009, among others). This phenomena is seen as puzzling and "paradoxical" because it seems in contradiction with the idea that, within the intertemporal choice model, which is the backbone of modern economics, when preferences are convex, consumption smoothing is the rule. Then, explanation of this paradox has been searched in relaxing some assumptions of the model of a fully rational forward looking agent. For example the agent may systematically underestimate the drop in earnings associated with retirement (Hamermesh, 1984a). Or, the agent may not be fully time consistent as in the hyperbolic discounting model [START_REF] Angeletos | The hyperbolic consumption model: Calibration, simulation, and empirical evaluation[END_REF].

Without denying that those phenomena may be important traits of "real" agents behavior, building on an insight of [START_REF] Banks | Is there a retirement-savings puzzle?[END_REF] in their conclusion, this paper will emphasize the point that a closer look at the intertemporal choice model of consumption and savings in continuous time allows to understand that what is smooth in the model is not necessarily consumption, but marginal utility of consumption. Of course, if consumption is the only variable of the utility function the two properties are equivalent. But if utility is multi-variate, any discontinuity in a dimension, may imply an optimal discontinuity response in the others. I will illustrate that insight into a very general model of inter-temporal choice that can be considered as a realistic generalisation of the basic one. Two ingredients will be required. First, I will assume a bi-variate, additively intertemporaly separable utility function that depends on consumption and leisure. Second I will assume, realistically, that retirement is not a smooth process with a per period duration of labor that tend progressively to zero, but a discontinuous process.

I will show that, as long as the per period utility function is not additively separable in consumption and leisure, discontinuity of the consumption function is the rule in this general model. However, as insightful is the preceding statement, it is not so easy to prove it formally with all generality because the assumptions imply a discontinuous payoff function, a case that is not standard with usual intertemporal optimization techniques in continuous time. I will provide a general and simple lemma that will make the problem tractable and it's resolution at the same time rigourous and insightful.

So, if we want to solve the paradox within a quite standard model of intertemporal choice, we have to drop additive separability of utility of consumption and leisure. And if we want to extend the problem to the choice of the optimal retirement age, we have to carry on with this non-separability. However, as pointed by d 'Albis et al. (2012) most of the study addressing the question has been made precisely under the assumption of additive separability in consumption and leisure (see d [START_REF] Albis | Endogenous retirement and monetary cycles[END_REF][START_REF] Bloom | Optimal retirement with increasing longevity[END_REF][START_REF] Boucekkine | Vintage human capital, demographic trends, and endogenous growth[END_REF][START_REF] Hazan | Longevity and lifetime labor supply: Evidence and implications[END_REF][START_REF] Heijdra | The individual life-cycle, annuity market imperfections and economic growth[END_REF][START_REF] Heijdra | Retirement, pensions, and ageing[END_REF][START_REF] Kalemli-Ozcan | Mortality change, the uncertainty effect, and retirement[END_REF][START_REF] Prettner | Increasing life expectancy and optimal retirement in general equilibrium[END_REF]Sheshinski, 1978, among others). And if there are some important papers that study a general life cycle model of consumption and savings, without additive separability of consumption and leisure [START_REF] Heckman | Life cycle consumption and labor supply: An explanation of the relationship between income and consumption over the life cycle[END_REF][START_REF] Heckman | A life-cycle model of earnings, learning, and consumption[END_REF][START_REF] Bütler | Neoclassical life-cycle consumption: a textbook example[END_REF] they are mostly focused on the the explanation of co-movement of earnings and consumption all over the life-cycle. Hamermesh (1984b) and [START_REF] Chang | Uncertain lifetimes, retirement and economic welfare[END_REF] study the retirement decision with non separability of consumption and leisure, but they fully endogenize the work decision, without any granularity concerning the per-period duration of worktime, and thus without any discontinuity of per period labor supply, implying model that are unable to explain at the same time retirement consumption paradox and the retirement decision. The model I propose can easily be expanded to endogenize the retirement decision and provide very general condition that fulfills the optimal age of retirement.

I will show that when optimal consumption is discontinuous at the age of retirement, this condition is qualitatively very different than in the traditional case.

consumption paradox

Let's assume that we are in a very standard continuous time life-cycle model of consumption and savings with preference for leisure and retirement.

P            max c T t e -θ(s-t) u (c (s) , l(s)) ds s.t.∀s ∈ [t, T ], ȧ(s) = ra(s) + w(s)(1 -l(s)) + b(s) -c(s) a(t)
given and a(T ) ≥ 0 t is the decision date and T is life duration, u is the per-period bi-variate utility function that depends on consumption and leisure. c, is the intertemporal consumption profile, the control variable of the program, l, is the intertemporal leisure profile, that I will assume, in a first stance, to be exogenous. a, is a life-cycle asset, the state variable of the program, that brings interest at the rate r. w, is labor income per period when the individual spend all this time working. b is social security income profile, interpreted as social security benefit when positive (typically after retirement) and social security contribution when negative (typically before retirement). l, c, w and b are assumed to be piecewise continuous and a is assumed to be piecewise smooth, assumptions that are fully compatible with the use of standard optimal control theory. I assume that the utility function includes standard minimum requirements of the microeconomic theory of consumption/leisure trade-off: u 1 > 0, u 2 > 0, u 11 < 0, u 22 < 0 and quasi-concavity (i.e. the indifference curves are convex). It implies that -

u 11 u 2 2 + 2u 12 u 1 u 2 -u 22 u 2 1 > 0.
It is important to notice that, without further assumptions, the sign of the second order crossed derivative is undetermined.

I will assume that there exists a retirement age t R such that:

     ∀s ∈ [t, t R ), l(s) = κ < 1 ∀s ∈ [t R , T ], l(s) = 1
Of course this assumption is a simplification, but it allows to characterize directly the central idea of the paper: retirement is fundamentally a discontinuity in the labor/leisure profile. This assumption seems much more realistic than usual idea that retirement is the smooth process with per period work duration tending to zero at the age of retirement.1

I denote c * the optimal consumption profile, solution of the program P and a * the associated value of the state variable. Of course those optimal functions are parameterized by all the given of the problem (t, t R , T, a(t), r, w, b, l).

I denote V * the optimal value of the problem i. e.

V * (t, t R , T, a(t), r, w) = T t e -θ(s-t) u (c * (t, t R , T, a(t), r, w, b, l, s), l(s)) ds
Because of the discontinuity of the instantaneous payoff function in t R , the problem is non standard. Therefore it is useful to decompose the problem in two separate ones:

P 0 P 1                  max c t R t e -θ(s-t) u (c (s) , κ) ds s.t. ȧ(s) = ra(s) + (1 -κ)w(s) + b(s) -c(s) a(t), a(t R ) given                  max c T t R e -θ(s-t) u (c (s) , 1)) ds s.t. ȧ(s) = ra(s) + b(s) -c(s) a(t R )
given and a(T ) ≥ 0 P 0 and P 1 . As c * they are also implicit functions of the parameter of their respective program and I can define the value function of P 0 and P 1 .

V 0 (t, t R , a(t), a(t R ), r, w, b, κ) = t R t e -θ(s-t) u c 0 (t, t R , a(t), a(t R ), r, w, b, κ), κ ds V 1 (t R , T, a(t R ), r, b) = T t R e -θ(s-t) u c 1 (t R , T, a(t R ), r, b, s), 1 ds
The two programs are linked by the asset level at the age of retirement. By application of the optimality principle, I can deduce:

Lemma 1 (A Principle of Optimality).

If (c * , a * ) is an admissible pair solution of program P then we have: [START_REF] Bellman | Dynamic programming[END_REF] principle of optimality.

1. V * (t, t R , T, a(t), r, w) = V 0 * (t, t R , a(t), a * (t R ), r, w) + V 1 * (t R , T, a * (t R ), r, w) 2. a * (t R ) = argmax a(t R ) {V 0 (t, t R , a(t), a(t R ), r, w) + V 1 (t R , T, a(t R ), r, w)} Proof: It is a direct application of
I have now all the material to solve the program P.

Proposition 1 (Discontinuity of the consumption profile).

If I denote c 0 (t R ) def = lim s→t R c 0 (s)
, and restrict my analysis to per period utility with a second order cross derivative that is either, everywhere strictly positive, everywhere strictly negative or everywhere equal to zero:

1. The optimal consumption profile solution of program P is unique.

The optimal consumption profile solution of program

P is continuous for every age s in [t, t R ) (t R , T ]. 3. In t R , u 1 (c 0 (t R ), κ) = u 1 (c 1 (t R ), 1
) and the continuity of the optimal consumption profile is determined solely the cross derivative of the per period utility function.

(a) c 0 (t R ) > c 1 (t R ) ⇔ u 12 (c, l) < 0 (b) c 0 (t R ) = c 1 (t R ) ⇔ u 12 (c, l) = 0 (c) c 0 (t R ) < c 1 (t R ) ⇔ u 12 (c, l) > 0
Proof: Relying on Lemma 1, I start by solving the program P 0 and P 1 for a given a(t R ). Denoting µ 0 the costate variable, the Hamiltonian of the Program P 0 is:

H 0 (c(s), a(s), µ 0 (s), s) = e -θ(s-t) u (c (s) , κ) + µ 0 (s) [r a(s) + (1 -κ)w(s) + b(s) -c(s)] (1)
According to Pontryagin maximum principle the necessary condition for optimality is:

∀s ∈ [t, t R ), ∂H 0 (•) ∂c(s) = 0 ⇒ µ 0 (s) = e -θ(s-t) u 1 (c (s) , κ) (2) ∀s ∈ [t, t R ), ∂H 0 (•) ∂a(s) = -μ0 (s) ⇒ μ0 (s) = -r µ 0 (s) (3) ∀s ∈ [t, t R ), ȧ(s) = ra(s) + (1 -κ)w(s) + b(s) -c(s) (4)
Moreover by construction of the Hamiltonian and Pontryagin maximum principle it is well known that:

∂V 0 (t, t R , a(t), a(t R ), r, w, b, κ) ∂a(t R ) = -µ 0 (t R ) (5)
Similarly for program P 1 , we have:

H 1 (c(s), a(s), µ 1 (s), s) = e -θ(s-t) u (c (s) , 1) + µ 1 (s) [r a(s) + b(s) -c(s)] (6) ∀s ∈ (t R , T ], ∂H 1 (•) ∂c(s) = 0 ⇒ µ 1 (s) = e -θ(s-t) u 1 (c (s) , 1) (7) ∀s ∈ (t R , T ], ∂H 1 (•) ∂a(s) = -μ1 (s) ⇒ μ1 (s) = -r µ 1 (s) (8) ∀s ∈ (t R , T ], ȧ(s) = ra(s) + b(s) -c(s) (9) ∂V 1 (t, t R , a(t), a(t R ), r, b) ∂a(t R ) = µ 1 (t R ) (10)
Moreover, P 1 being a constrained endpoint problem, we have to fulfill the transversality condition:

µ 1 (T )a(T ) = 0 ⇒ a(T ) = 0 (11) P 0 and P 1 verifying the standard strict concavity condition of their respective Hamiltonian, they both admit continuous and unique solution on their respective domain.

Let us now turn to the solution problem of the optimal value of the asset at retirement date a * (t R ). Relaying on the principle of optimality (Lemma 1), a necessary condition for a * (t R ) to be a maximum of (V

0 (•) + V 1 (•)) is: ∂V 0 (•) ∂a(t R ) + V 1 (•) ∂a(t R ) = -µ 0 (t R ) + µ 1 (t R ) = 0 ⇔ u 1 (c 0 (t R ), κ) = u 1 (c 1 (t R ), 1) (12)
It is easy to check that the left hand term of the last equality is increasing in a(t R )

while the right hand one is decreasing, assuring the uniqueness of a

* (t R ). If for all c, l in R + × [0, 1], u 12 < 0, then u 1 (c 0 (t R ), κ) < u 1 (c 0 (t R ), 1). Because u 11 < 0, we can only have u 1 (c 0 (t R ), κ) = u 1 (c 1 (t R ), 1), if and only if c 0 (t R ) > c 1 (t R ).
The reasoning is the same for the two other cases.

In this setting, a negative cross derivative of the per period utility of consumption and leisure is necessary to obtain a discontinuous drop in consumption at the age of retirement, i.e. to resolve the retirement consumption puzzle. It means that, if we believe that the model is a proper simplification of the intertemporal choice of agent in the real world, the observation of that kind of drop, informs us on the negative sign of the cross derivative. It may seems strange because many workhorse utility function in labor economics such as the cobb-Douglas or the CES utility function are characterized by a positive cross derivative.

However, it is important to notice that relying on a different model of intertemporal choice with full endogeneity of labor, [START_REF] Heckman | Life cycle consumption and labor supply: An explanation of the relationship between income and consumption over the life cycle[END_REF] also conclude that a negative cross derivative of the per period utility of consumption and leisure was required to explain the hump shape of the intertemporal consumption profile.

In this part, I have given a complete theoretical treatment of an idea that was alluded in [START_REF] Banks | Is there a retirement-savings puzzle?[END_REF] and in the "back-of-the-envelope calculation" in [START_REF] Battistin | The retirement consumption puzzle: Evidence from a regression discontinuity approach[END_REF]. This calculation was grounded on the following parametrical form:

u(c, l) = (c α l 1-α ) 1-γ 1 -γ
with γ > 0 interpreted as the reciprocal of the intertemporal elasticity of substitution. They rightfully conclude that to solve the retirement consumption puzzle in this model, γ > 1 is required, but they miss the right insight for explaining that. Because, in this model, γ fully capture the intensity of the response of consumption to a variation of the rate of interest only when leisure is fully endogenous, but in this case there will be no discontinuity in the consumption function. As we have shown, explaining such a discontinuity, requires leisure to be exogenous at the age of retirement 2 , then it is -c u 11 /u 1 = α(γ -1) + 1 that will capture the intensity of response of consumption to a change of the rate of interest. Moreover, if the model is based on a Cobb-Douglass utility function, it is in fact a power transformation of a Cobb-Douglass, a transformation that can alter the sign of the second order cross derivative.

We have u 12 = α(1 -α)(1 -γ)c α(1-γ)-1 l (1-α)(1-γ)-1 . With this special parametrical form, the sign of the cross derivative of utility is fully given by the position of γ with respect to unity. When γ is higher than one, this cross derivative is negative explaining the downward discontinuity in consumption, as confirmed by the general statement of Proposition 1.3. The effect has nothing to do with the the intertemporal elasticity of substitution per se.

2 Or at least a constraint for a minimum per-period work duration that is binding.

3 Optimal age of retirement I have solved the program P with the age of retirement, t R , being a parameter. I have all the material to characterize the optimal age of retirement, the one that maximizes the value of the program. In particular, the decomposition of the general Program in two sub-programs delimited by the age of retirement, allows to derive this optimal age of retirement in a parsimonious and elegant manner.

Proposition 2 (The optimal age of retirement).

When an interior solution exists, and denoting b

0 (t R ) def = lim s→t - R b (
s) < 0, the optimal age of retirement tR is such that:

u(c 1 ( tR ), 1) -u(c 0 ( tR ), κ) = u 1 (c 0 (t R ), κ) ((1 -κ)w( tR ) + b 0 ( tR )) -b( tR )) + (c 1 ( tR ) -c 0 ( tR )) (13) Proof: tR is a solution of max t R
V * (t, t R , T, a(t), r, w). Because V * is continuous and differentiable in t R , a necessary condition for having an interior solution is:

∂V * (t, t R , T, a(t), r, w) ∂t R = 0 (14)
Relying on Lemma 1 and noting that by construction of the Hamiltonian and Pontryagin maximum principle:

∂V 0 (t, t R , a(t), r, w) ∂t R = H 0 (c 0 (t R ), a 0 (t R ), µ 0 (t R ), t R ) and ∂V 1 (t R , T, a(t), r, w) ∂t R = -H 1 (c 1 (t R ), a 1 (t R ), µ 1 (t R ), t R )
we can easily conclude that tR is such that:

H 0 (c 0 ( tR ), a * ( tR ), µ 0 ( tR ), tR ) = H 1 (c 1 ( tR ), a * ( tR ), µ 1 ( tR ), tR ) (15)
Using the definitions of the Hamiltonian and first order conditions of program P 0 and P 1 and remembering that, in any case, a is continuous in t R , we get the right hand side.

This is a standard marginal condition for optimality. The left hand side of Equation ( 13) is the direct cost in utility of a marginal increase in the retirement age, while the right hand side is the indirect gain in utility due to supplementary resources generated by a longer work duration. The important and innovative point is that when taking into account the retirement consumption puzzle, the endogenous drop of consumption implies that less resources are required to maintain a same level of utility. Thus the earnings differential can be higher when the agents decide to retire.

Proposition 2 provides a very general characterisation of the optimal retirement age.

Moreover, when expanding consumption before and after retirement as implicit function of the parameters of the problem, and when endogenizing the budgetary constraint of the social security system, it allows to derive comparative static results on the optimal age of retirement.

Conclusion

This short paper provides a general methodology to resolve the retirement consumption puzzle and the choice of the optimal age of retirement. The principle is illustrated in a simple model of intertemporal choice in which utility depend on consumption and leisure with certain horizon. To solve the puzzle we need only two assumptions: 1.

retirement implies a discontinuity in the leisure intertemporal profile and, 2. the crossderivative of the utility function is negative. The method is general and can easily be extended in more realistic models with uncertain lifetime3 .

This idea could be generalized by endogenizing per period work duration taking into account a granularity assumption. In general for organizational reason work duration can be zero or something significantly different from zero.

In a companion paper, I am actually working on a calibrated version of the model taking into account a realistic modeling of uncertain lifetime in the spirit of[START_REF] Drouhin | A rank-dependent utility model of uncertain lifetime[END_REF] and the possibility of a non stationary intertemporal utility, allowing for per period utility to change with age in the spirit of[START_REF] Drouhin | Non stationary additive utility and time consistency[END_REF]