1 \documentclass{article}
2 %\documentclass[10pt,journal,letterpaper,compsoc]{IEEEtran}
3 %\documentclass[preprint,12pt]{elsarticle}
4 \usepackage[utf8]{inputenc}
5 \usepackage[T1]{fontenc}
12 \usepackage[ruled,vlined]{algorithm2e}
14 \usepackage[standard]{ntheorem}
15 \usepackage{algorithmic}
23 % Pour mathds : les ensembles IR, IN, etc.
26 % Pour avoir des intervalles d'entiers
30 % Pour faire des sous-figures dans les figures
31 \usepackage{subfigure}
34 \newtheorem{notation}{Notation}
36 \newcommand{\X}{\mathcal{X}}
37 \newcommand{\Go}{G_{f_0}}
38 \newcommand{\B}{\mathds{B}}
39 \newcommand{\N}{\mathds{N}}
40 \newcommand{\BN}{\mathds{B}^\mathsf{N}}
43 \newcommand{\alert}[1]{\begin{color}{blue}\textit{#1}\end{color}}
47 \title{Efficient and Cryptographically Secure Generation of Chaotic Pseudorandom Numbers on GPU}
50 %% \author{Jacques M. Bahi}
51 %% \ead{jacques.bahi@univ-fcomte.fr}
52 %% \author{ Rapha\"{e}l Couturier \corref{cor1}}
53 %% \ead{raphael.couturier@univ-fcomte.fr}
54 %% \cortext[cor1]{Corresponding author}
55 %% \author{ Christophe Guyeux}
56 %% \ead{christophe.guyeux@univ-fcomte.fr}
57 %% \author{ Pierre-Cyrille Héam }
58 %% \ead{pierre-cyrille.heam@univ-fcomte.fr}
60 \author{Christophe Guyeux \and Rapha\"{e}l Couturier \and Pierre-Cyrille Héam \and Jacques M. Bahi\\
61 FEMTO-ST Institute, UMR 6174 CNRS,\\ University of Franche Comte, Belfort, France}
67 %\IEEEcompsoctitleabstractindextext{
69 In this paper we present a new pseudorandom number generator (PRNG) on
70 graphics processing units (GPU). This PRNG is based on the so-called chaotic iterations. It
71 is firstly proven to be chaotic according to the Devaney's formulation. We thus propose an efficient
72 implementation for GPU that successfully passes the {\it BigCrush} tests, deemed to be the hardest
73 battery of tests in TestU01. Experiments show that this PRNG can generate
74 about 20 billion of random numbers per second on Tesla C1060 and NVidia GTX280
76 It is then established that, under reasonable assumptions, the proposed PRNG can be cryptographically
78 A chaotic version of the Blum-Goldwasser asymmetric key encryption scheme is finally proposed.
84 % pseudo random number\sep parallelization\sep GPU\sep cryptography\sep chaos
89 %\IEEEdisplaynotcompsoctitleabstractindextext
90 %\IEEEpeerreviewmaketitle
93 \section{Introduction}
95 Randomness is of importance in many fields such as scientific simulations or cryptography.
96 ``Random numbers'' can mainly be generated either by a deterministic and reproducible algorithm
97 called a pseudorandom number generator (PRNG), or by a physical non-deterministic
98 process having all the characteristics of a random noise, called a truly random number
100 In this paper, we focus on reproducible generators, useful for instance in
101 Monte-Carlo based simulators or in several cryptographic schemes.
102 These domains need PRNGs that are statistically irreproachable.
103 In some fields such as in numerical simulations, speed is a strong requirement
104 that is usually attained by using parallel architectures. In that case,
105 a recurrent problem is that a deflation of the statistical qualities is often
106 reported, when the parallelization of a good PRNG is realized.
107 This is why ad-hoc PRNGs for each possible architecture must be found to
108 achieve both speed and randomness.
109 On the other side, speed is not the main requirement in cryptography: the great
110 need is to define \emph{secure} generators able to withstand malicious
111 attacks. Roughly speaking, an attacker should not be able in practice to make
112 the distinction between numbers obtained with the secure generator and a true random
113 sequence. Or, in an equivalent formulation, he or she should not be
114 able (in practice) to predict the next bit of the generator, having the knowledge of all the
115 binary digits that have been already released. ``Being able in practice'' refers here
116 to the possibility to achieve this attack in polynomial time, and to the exponential growth
117 of the difficulty of this challenge when the size of the parameters of the PRNG increases.
120 Finally, a small part of the community working in this domain focuses on a
121 third requirement, that is to define chaotic generators.
122 The main idea is to take benefits from a chaotic dynamical system to obtain a
123 generator that is unpredictable, disordered, sensible to its seed, or in other word chaotic.
124 Their desire is to map a given chaotic dynamics into a sequence that seems random
125 and unassailable due to chaos.
126 However, the chaotic maps used as a pattern are defined in the real line
127 whereas computers deal with finite precision numbers.
128 This distortion leads to a deflation of both chaotic properties and speed.
129 Furthermore, authors of such chaotic generators often claim their PRNG
130 as secure due to their chaos properties, but there is no obvious relation
131 between chaos and security as it is understood in cryptography.
132 This is why the use of chaos for PRNG still remains marginal and disputable.
134 The authors' opinion is that topological properties of disorder, as they are
135 properly defined in the mathematical theory of chaos, can reinforce the quality
136 of a PRNG. But they are not substitutable for security or statistical perfection.
137 Indeed, to the authors' mind, such properties can be useful in the two following situations. On the
138 one hand, a post-treatment based on a chaotic dynamical system can be applied
139 to a PRNG statistically deflective, in order to improve its statistical
140 properties. Such an improvement can be found, for instance, in~\cite{bgw09:ip,bcgr11:ip}.
141 On the other hand, chaos can be added to a fast, statistically perfect PRNG and/or a
142 cryptographically secure one, in case where chaos can be of interest,
143 \emph{only if these last properties are not lost during
144 the proposed post-treatment}. Such an assumption is behind this research work.
145 It leads to the attempts to define a
146 family of PRNGs that are chaotic while being fast and statistically perfect,
147 or cryptographically secure.
148 Let us finish this paragraph by noticing that, in this paper,
149 statistical perfection refers to the ability to pass the whole
150 {\it BigCrush} battery of tests, which is widely considered as the most
151 stringent statistical evaluation of a sequence claimed as random.
152 This battery can be found in the well-known TestU01 package~\cite{LEcuyerS07}.
153 More precisely, each time we performed a test on a PRNG, we ran it
154 twice in order to observe if all $p-$values are inside [0.01, 0.99]. In
155 fact, we observed that few $p-$values (less than ten) are sometimes
156 outside this interval but inside [0.001, 0.999], so that is why a
157 second run allows us to confirm that the values outside are not for
158 the same test. With this approach all our PRNGs pass the {\it
159 BigCrush} successfully and all $p-$values are at least once inside
161 Chaos, for its part, refers to the well-established definition of a
162 chaotic dynamical system defined by Devaney~\cite{Devaney}.
164 In a previous work~\cite{bgw09:ip,guyeux10} we have proposed a post-treatment on PRNGs making them behave
165 as a chaotic dynamical system. Such a post-treatment leads to a new category of
166 PRNGs. We have shown that proofs of Devaney's chaos can be established for this
167 family, and that the sequence obtained after this post-treatment can pass the
168 NIST~\cite{Nist10}, DieHARD~\cite{Marsaglia1996}, and TestU01~\cite{LEcuyerS07} batteries of tests, even if the inputted generators
170 The proposition of this paper is to improve widely the speed of the formerly
171 proposed generator, without any lack of chaos or statistical properties.
172 In particular, a version of this PRNG on graphics processing units (GPU)
174 Although GPU was initially designed to accelerate
175 the manipulation of images, they are nowadays commonly used in many scientific
176 applications. Therefore, it is important to be able to generate pseudorandom
177 numbers inside a GPU when a scientific application runs in it. This remark
178 motivates our proposal of a chaotic and statistically perfect PRNG for GPU.
180 allows us to generate almost 20 billion of pseudorandom numbers per second.
181 Furthermore, we show that the proposed post-treatment preserves the
182 cryptographical security of the inputted PRNG, when this last has such a
184 Last, but not least, we propose a rewriting of the Blum-Goldwasser asymmetric
185 key encryption protocol by using the proposed method.
188 {\bf Main contributions.} In this paper a new PRNG using chaotic iteration
189 is defined. From a theoretical point of view, it is proven that it has fine
190 topological chaotic properties and that it is cryptographically secured (when
191 the initial PRNG is also cryptographically secured). From a practical point of
192 view, experiments point out a very good statistical behavior. An optimized
193 original implementation of this PRNG is also proposed and experimented.
194 Pseudorandom numbers are generated at a rate of 20GSamples/s, which is faster
195 than in~\cite{conf/fpga/ThomasHL09,Marsaglia2003} (and with a better
196 statistical behavior). Experiments are also provided using
197 \begin{color}{red} the well-known Blum-Blum-Shub
201 random generator. The generation speed is significantly weaker.
202 %Note also that an original qualitative comparison between topological chaotic
203 %properties and statistical tests is also proposed.
208 The remainder of this paper is organized as follows. In Section~\ref{section:related
209 works} we review some GPU implementations of PRNGs. Section~\ref{section:BASIC
210 RECALLS} gives some basic recalls on the well-known Devaney's formulation of chaos,
211 and on an iteration process called ``chaotic
212 iterations'' on which the post-treatment is based.
213 The proposed PRNG and its proof of chaos are given in Section~\ref{sec:pseudorandom}.
214 %Section~\ref{The generation of pseudorandom sequence} illustrates the statistical
215 %improvement related to the chaotic iteration based post-treatment, for
216 %our previously released PRNGs and a new efficient
217 %implementation on CPU.
218 Section~\ref{sec:efficient PRNG} %{sec:efficient PRNG
220 describes and evaluates theoretically new effective versions of
221 our pseudorandom generators, in particular with a GPU implementation.
222 Such generators are experimented in
223 Section~\ref{sec:experiments}.
224 We show in Section~\ref{sec:security analysis} that, if the inputted
225 generator is cryptographically secure, then it is the case too for the
226 generator provided by the post-treatment.
228 security evaluation is also outlined in Section~\ref{sec:Practicak evaluation}.
229 Such a proof leads to the proposition of a cryptographically secure and
230 chaotic generator on GPU based on the famous Blum Blum Shub
231 in Section~\ref{sec:CSGPU} and to an improvement of the
232 Blum-Goldwasser protocol in Sect.~\ref{Blum-Goldwasser}.
233 This research work ends by a conclusion section, in which the contribution is
234 summarized and intended future work is presented.
239 \section{Related work on GPU based PRNGs}
240 \label{section:related works}
242 Numerous research works on defining GPU based PRNGs have already been proposed in the
243 literature, so that exhaustivity is impossible.
244 This is why authors of this document only give reference to the most significant attempts
245 in this domain, from their subjective point of view.
246 The quantity of pseudorandom numbers generated per second is mentioned here
247 only when the information is given in the related work.
248 A million numbers per second will be simply written as
249 1MSample/s whereas a billion numbers per second is 1GSample/s.
251 In \cite{Pang:2008:cec} a PRNG based on cellular automata is defined
252 with no requirement to an high precision integer arithmetic or to any bitwise
253 operations. Authors can generate about
254 3.2MSamples/s on a GeForce 7800 GTX GPU, which is quite an old card now.
255 However, there is neither a mention of statistical tests nor any proof of
256 chaos or cryptography in this document.
258 In \cite{ZRKB10}, the authors propose different versions of efficient GPU PRNGs
259 based on Lagged Fibonacci or Hybrid Taus. They have used these
260 PRNGs for Langevin simulations of biomolecules fully implemented on
261 GPU. Performances of the GPU versions are far better than those obtained with a
262 CPU, and these PRNGs succeed to pass the {\it BigCrush} battery of TestU01.
263 However the evaluations of the proposed PRNGs are only statistical ones.
266 Authors of~\cite{conf/fpga/ThomasHL09} have studied the implementation of some
267 PRNGs on different computing architectures: CPU, field-programmable gate array
268 (FPGA), massively parallel processors, and GPU. This study is of interest, because
269 the performance of the same PRNGs on different architectures are compared.
270 FPGA appears as the fastest and the most
271 efficient architecture, providing the fastest number of generated pseudorandom numbers
273 However, we notice that authors can ``only'' generate between 11 and 16GSamples/s
274 with a GTX 280 GPU, which should be compared with
275 the results presented in this document.
276 We can remark too that the PRNGs proposed in~\cite{conf/fpga/ThomasHL09} are only
277 able to pass the {\it Crush} battery, which is far easier than the {\it Big Crush} one.
279 Lastly, Cuda has developed a library for the generation of pseudorandom numbers called
280 Curand~\cite{curand11}. Several PRNGs are implemented, among
282 Xorwow~\cite{Marsaglia2003} and some variants of Sobol. The tests reported show that
283 their fastest version provides 15GSamples/s on the new Fermi C2050 card.
284 But their PRNGs cannot pass the whole TestU01 battery (only one test is failed).
287 We can finally remark that, to the best of our knowledge, no GPU implementation has been proven to be chaotic, and the cryptographically secure property has surprisingly never been considered.
289 \section{Basic Recalls}
290 \label{section:BASIC RECALLS}
292 This section is devoted to basic definitions and terminologies in the fields of
293 topological chaos and chaotic iterations. We assume the reader is familiar
294 with basic notions on topology (see for instance~\cite{Devaney}).
297 \subsection{Devaney's Chaotic Dynamical Systems}
298 \label{subsec:Devaney}
299 In the sequel $S^{n}$ denotes the $n^{th}$ term of a sequence $S$ and $V_{i}$
300 denotes the $i^{th}$ component of a vector $V$. $f^{k}=f\circ ...\circ f$
301 is for the $k^{th}$ composition of a function $f$. Finally, the following
302 notation is used: $\llbracket1;N\rrbracket=\{1,2,\hdots,N\}$.
305 Consider a topological space $(\mathcal{X},\tau)$ and a continuous function $f :
306 \mathcal{X} \rightarrow \mathcal{X}$.
309 The function $f$ is said to be \emph{topologically transitive} if, for any pair of open sets
310 $U,V \subset \mathcal{X}$, there exists $k>0$ such that $f^k(U) \cap V \neq
315 An element $x$ is a \emph{periodic point} for $f$ of period $n\in \mathds{N}^*$
316 if $f^{n}(x)=x$.% The set of periodic points of $f$ is denoted $Per(f).$
320 $f$ is said to be \emph{regular} on $(\mathcal{X}, \tau)$ if the set of periodic
321 points for $f$ is dense in $\mathcal{X}$: for any point $x$ in $\mathcal{X}$,
322 any neighborhood of $x$ contains at least one periodic point (without
323 necessarily the same period).
327 \begin{definition}[Devaney's formulation of chaos~\cite{Devaney}]
328 The function $f$ is said to be \emph{chaotic} on $(\mathcal{X},\tau)$ if $f$ is regular and
329 topologically transitive.
332 The chaos property is strongly linked to the notion of ``sensitivity'', defined
333 on a metric space $(\mathcal{X},d)$ by:
336 \label{sensitivity} The function $f$ has \emph{sensitive dependence on initial conditions}
337 if there exists $\delta >0$ such that, for any $x\in \mathcal{X}$ and any
338 neighborhood $V$ of $x$, there exist $y\in V$ and $n > 0$ such that
339 $d\left(f^{n}(x), f^{n}(y)\right) >\delta $.
341 The constant $\delta$ is called the \emph{constant of sensitivity} of $f$.
344 Indeed, Banks \emph{et al.} have proven in~\cite{Banks92} that when $f$ is
345 chaotic and $(\mathcal{X}, d)$ is a metric space, then $f$ has the property of
346 sensitive dependence on initial conditions (this property was formerly an
347 element of the definition of chaos). To sum up, quoting Devaney
348 in~\cite{Devaney}, a chaotic dynamical system ``is unpredictable because of the
349 sensitive dependence on initial conditions. It cannot be broken down or
350 simplified into two subsystems which do not interact because of topological
351 transitivity. And in the midst of this random behavior, we nevertheless have an
352 element of regularity''. Fundamentally different behaviors are consequently
353 possible and occur in an unpredictable way.
357 \subsection{Chaotic Iterations}
358 \label{sec:chaotic iterations}
361 Let us consider a \emph{system} with a finite number $\mathsf{N} \in
362 \mathds{N}^*$ of elements (or \emph{cells}), so that each cell has a
363 Boolean \emph{state}. Having $\mathsf{N}$ Boolean values for these
364 cells leads to the definition of a particular \emph{state of the
365 system}. A sequence which elements belong to $\llbracket 1;\mathsf{N}
366 \rrbracket $ is called a \emph{strategy}. The set of all strategies is
367 denoted by $\llbracket 1, \mathsf{N} \rrbracket^\mathds{N}.$
370 \label{Def:chaotic iterations}
371 The set $\mathds{B}$ denoting $\{0,1\}$, let
372 $f:\mathds{B}^{\mathsf{N}}\longrightarrow \mathds{B}^{\mathsf{N}}$ be
373 a function and $S\in \llbracket 1, \mathsf{N} \rrbracket^\mathds{N}$ be a ``strategy''. The so-called
374 \emph{chaotic iterations} are defined by $x^0\in
375 \mathds{B}^{\mathsf{N}}$ and
377 \forall n\in \mathds{N}^{\ast }, \forall i\in
378 \llbracket1;\mathsf{N}\rrbracket ,x_i^n=\left\{
380 x_i^{n-1} & \text{ if }S^n\neq i \\
381 \left(f(x^{n-1})\right)_{S^n} & \text{ if }S^n=i.
386 In other words, at the $n^{th}$ iteration, only the $S^{n}-$th cell is
387 \textquotedblleft iterated\textquotedblright . Note that in a more
388 general formulation, $S^n$ can be a subset of components and
389 $\left(f(x^{n-1})\right)_{S^{n}}$ can be replaced by
390 $\left(f(x^{k})\right)_{S^{n}}$, where $k<n$, describing for example,
391 delays transmission~\cite{Robert1986,guyeux10}. Finally, let us remark that
392 the term ``chaotic'', in the name of these iterations, has \emph{a
393 priori} no link with the mathematical theory of chaos, presented above.
396 Let us now recall how to define a suitable metric space where chaotic iterations
397 are continuous. For further explanations, see, e.g., \cite{guyeux10}.
399 Let $\delta $ be the \emph{discrete Boolean metric}, $\delta
400 (x,y)=0\Leftrightarrow x=y.$ Given a function $f$, define the function
401 $F_{f}: \llbracket1;\mathsf{N}\rrbracket\times \mathds{B}^{\mathsf{N}}
402 \longrightarrow \mathds{B}^{\mathsf{N}}$
405 & (k,E) & \longmapsto & \left( E_{j}.\delta (k,j)+ f(E)_{k}.\overline{\delta
406 (k,j)}\right) _{j\in \llbracket1;\mathsf{N}\rrbracket}%
409 \noindent where + and . are the Boolean addition and product operations.
410 Consider the phase space:
412 \mathcal{X} = \llbracket 1 ; \mathsf{N} \rrbracket^\mathds{N} \times
413 \mathds{B}^\mathsf{N},
415 \noindent and the map defined on $\mathcal{X}$:
417 G_f\left(S,E\right) = \left(\sigma(S), F_f(i(S),E)\right), \label{Gf}
419 \noindent where $\sigma$ is the \emph{shift} function defined by $\sigma
420 (S^{n})_{n\in \mathds{N}}\in \llbracket 1, \mathsf{N} \rrbracket^\mathds{N}\longrightarrow (S^{n+1})_{n\in
421 \mathds{N}}\in \llbracket 1, \mathsf{N} \rrbracket^\mathds{N}$ and $i$ is the \emph{initial function}
422 $i:(S^{n})_{n\in \mathds{N}} \in \llbracket 1, \mathsf{N} \rrbracket^\mathds{N}\longrightarrow S^{0}\in \llbracket
423 1;\mathsf{N}\rrbracket$. Then the chaotic iterations proposed in
424 Definition \ref{Def:chaotic iterations} can be described by the following iterations:
428 X^0 \in \mathcal{X} \\
434 With this formulation, a shift function appears as a component of chaotic
435 iterations. The shift function is a famous example of a chaotic
436 map~\cite{Devaney} but its presence is not sufficient enough to claim $G_f$ as
438 To study this claim, a new distance between two points $X = (S,E), Y =
439 (\check{S},\check{E})\in
440 \mathcal{X}$ has been introduced in \cite{guyeux10} as follows:
442 d(X,Y)=d_{e}(E,\check{E})+d_{s}(S,\check{S}),
448 \displaystyle{d_{e}(E,\check{E})} & = & \displaystyle{\sum_{k=1}^{\mathsf{N}%
449 }\delta (E_{k},\check{E}_{k})}, \\
450 \displaystyle{d_{s}(S,\check{S})} & = & \displaystyle{\dfrac{9}{\mathsf{N}}%
451 \sum_{k=1}^{\infty }\dfrac{|S^k-\check{S}^k|}{10^{k}}}.%
457 This new distance has been introduced to satisfy the following requirements.
459 \item When the number of different cells between two systems is increasing, then
460 their distance should increase too.
461 \item In addition, if two systems present the same cells and their respective
462 strategies start with the same terms, then the distance between these two points
463 must be small because the evolution of the two systems will be the same for a
464 while. Indeed, both dynamical systems start with the same initial condition,
465 use the same update function, and as strategies are the same for a while, furthermore
466 updated components are the same as well.
468 The distance presented above follows these recommendations. Indeed, if the floor
469 value $\lfloor d(X,Y)\rfloor $ is equal to $n$, then the systems $E, \check{E}$
470 differ in $n$ cells ($d_e$ is indeed the Hamming distance). In addition, $d(X,Y) - \lfloor d(X,Y) \rfloor $ is a
471 measure of the differences between strategies $S$ and $\check{S}$. More
472 precisely, this floating part is less than $10^{-k}$ if and only if the first
473 $k$ terms of the two strategies are equal. Moreover, if the $k^{th}$ digit is
474 nonzero, then the $k^{th}$ terms of the two strategies are different.
475 The impact of this choice for a distance will be investigated at the end of the document.
477 Finally, it has been established in \cite{guyeux10} that,
480 Let $f$ be a map from $\mathds{B}^\mathsf{N}$ to itself. Then $G_{f}$ is continuous in
481 the metric space $(\mathcal{X},d)$.
484 The chaotic property of $G_f$ has been firstly established for the vectorial
485 Boolean negation $f_0(x_1,\hdots, x_\mathsf{N}) = (\overline{x_1},\hdots, \overline{x_\mathsf{N}})$ \cite{guyeux10}. To obtain a characterization, we have secondly
486 introduced the notion of asynchronous iteration graph recalled bellow.
488 Let $f$ be a map from $\mathds{B}^\mathsf{N}$ to itself. The
489 {\emph{asynchronous iteration graph}} associated with $f$ is the
490 directed graph $\Gamma(f)$ defined by: the set of vertices is
491 $\mathds{B}^\mathsf{N}$; for all $x\in\mathds{B}^\mathsf{N}$ and
492 $i\in \llbracket1;\mathsf{N}\rrbracket$,
493 the graph $\Gamma(f)$ contains an arc from $x$ to $F_f(i,x)$.
494 The relation between $\Gamma(f)$ and $G_f$ is clear: there exists a
495 path from $x$ to $x'$ in $\Gamma(f)$ if and only if there exists a
496 strategy $s$ such that the parallel iteration of $G_f$ from the
497 initial point $(s,x)$ reaches the point $x'$.
498 We have then proven in \cite{bcgr11:ip} that,
502 \label{Th:Caractérisation des IC chaotiques}
503 Let $f:\mathds{B}^\mathsf{N}\to\mathds{B}^\mathsf{N}$. $G_f$ is chaotic (according to Devaney)
504 if and only if $\Gamma(f)$ is strongly connected.
507 Finally, we have established in \cite{bcgr11:ip} that,
509 Let $f: \mathds{B}^{n} \rightarrow \mathds{B}^{n}$, $\Gamma(f)$ its
510 iteration graph, $\check{M}$ its adjacency
512 a $n\times n$ matrix defined by
514 M_{ij} = \frac{1}{n}\check{M}_{ij}$ %\textrm{
516 $M_{ii} = 1 - \frac{1}{n} \sum\limits_{j=1, j\neq i}^n \check{M}_{ij}$ otherwise.
518 If $\Gamma(f)$ is strongly connected, then
519 the output of the PRNG detailed in Algorithm~\ref{CI Algorithm} follows
520 a law that tends to the uniform distribution
521 if and only if $M$ is a double stochastic matrix.
525 These results of chaos and uniform distribution have led us to study the possibility of building a
526 pseudorandom number generator (PRNG) based on the chaotic iterations.
527 As $G_f$, defined on the domain $\llbracket 1 ; \mathsf{N} \rrbracket^{\mathds{N}}
528 \times \mathds{B}^\mathsf{N}$, is built from Boolean networks $f : \mathds{B}^\mathsf{N}
529 \rightarrow \mathds{B}^\mathsf{N}$, we can preserve the theoretical properties on $G_f$
530 during implementations (due to the discrete nature of $f$). Indeed, it is as if
531 $\mathds{B}^\mathsf{N}$ represents the memory of the computer whereas $\llbracket 1 ; \mathsf{N}
532 \rrbracket^{\mathds{N}}$ is its input stream (the seeds, for instance, in PRNG, or a physical noise in TRNG).
533 Let us finally remark that the vectorial negation satisfies the hypotheses of both theorems above.
535 \section{Application to Pseudorandomness}
536 \label{sec:pseudorandom}
538 \subsection{A First Pseudorandom Number Generator}
540 We have proposed in~\cite{bgw09:ip} a new family of generators that receives
541 two PRNGs as inputs. These two generators are mixed with chaotic iterations,
542 leading thus to a new PRNG that
543 should improve the statistical properties of each
544 generator taken alone.
545 Furthermore, the generator obtained in this way possesses various chaos properties that none of the generators used as input present.
549 \begin{algorithm}[h!]
551 \KwIn{a function $f$, an iteration number $b$, an initial configuration $x^0$
553 \KwOut{a configuration $x$ ($n$ bits)}
555 $k\leftarrow b + PRNG_1(b)$\;
558 $s\leftarrow{PRNG_2(n)}$\;
559 $x\leftarrow{F_f(s,x)}$\;
563 \caption{An arbitrary round of $Old~ CI~ PRNG_f(PRNG_1,PRNG_2)$}
570 This generator is synthesized in Algorithm~\ref{CI Algorithm}.
571 It takes as input: a Boolean function $f$ satisfying Theorem~\ref{Th:Caractérisation des IC chaotiques};
572 an integer $b$, ensuring that the number of executed iterations
573 between two outputs is at least $b$
574 and at most $2b+1$; and an initial configuration $x^0$.
575 It returns the new generated configuration $x$. Internally, it embeds two
576 inputted generators $PRNG_i(k), i=1,2$,
577 which must return integers
578 uniformly distributed
579 into $\llbracket 1 ; k \rrbracket$.
580 For instance, these PRNGs can be the \textit{XORshift}~\cite{Marsaglia2003},
581 being a category of very fast PRNGs designed by George Marsaglia
582 that repeatedly uses the transform of exclusive or (XOR, $\oplus$) on a number
583 with a bit shifted version of it. Such a PRNG, which has a period of
584 $2^{32}-1=4.29\times10^9$, is summed up in Algorithm~\ref{XORshift}.
585 This XORshift, or any other reasonable PRNG, is used
586 in our own generator to compute both the number of iterations between two
587 outputs (provided by $PRNG_1$) and the strategy elements ($PRNG_2$).
589 %This former generator has successively passed various batteries of statistical tests, as the NIST~\cite{bcgr11:ip}, DieHARD~\cite{Marsaglia1996}, and TestU01~\cite{LEcuyerS07} ones.
592 \begin{algorithm}[h!]
594 \KwIn{the internal configuration $z$ (a 32-bit word)}
595 \KwOut{$y$ (a 32-bit word)}
596 $z\leftarrow{z\oplus{(z\ll13)}}$\;
597 $z\leftarrow{z\oplus{(z\gg17)}}$\;
598 $z\leftarrow{z\oplus{(z\ll5)}}$\;
602 \caption{An arbitrary round of \textit{XORshift} algorithm}
607 \subsection{A ``New CI PRNG''}
609 In order to make the Old CI PRNG usable in practice, we have proposed
610 an adapted version of the chaotic iteration based generator in~\cite{bg10:ip}.
611 In this ``New CI PRNG'', we prevent a given bit from changing twice between two outputs.
612 This new generator is designed by the following process.
614 First of all, some chaotic iterations have to be done to generate a sequence
615 $\left(x^n\right)_{n\in\mathds{N}} \in \left(\mathds{B}^{32}\right)^\mathds{N}$
616 of Boolean vectors, which are the successive states of the iterated system.
617 Some of these vectors will be randomly extracted and our pseudorandom bit
618 flow will be constituted by their components. Such chaotic iterations are
619 realized as follows. Initial state $x^0 \in \mathds{B}^{32}$ is a Boolean
620 vector taken as a seed and chaotic strategy $\left(S^n\right)_{n\in\mathds{N}}\in
621 \llbracket 1, 32 \rrbracket^\mathds{N}$ is
622 an \emph{irregular decimation} of $PRNG_2$ sequence, as described in
623 Algorithm~\ref{Chaotic iteration1}.
625 Then, at each iteration, only the $S^n$-th component of state $x^n$ is
626 updated, as follows: $x_i^n = x_i^{n-1}$ if $i \neq S^n$, else $x_i^n = \overline{x_i^{n-1}}$.
627 Such a procedure is equivalent to achieving chaotic iterations with
628 the Boolean vectorial negation $f_0$ and some well-chosen strategies.
629 Finally, some $x^n$ are selected
630 by a sequence $m^n$ as the pseudorandom bit sequence of our generator.
631 $(m^n)_{n \in \mathds{N}} \in \mathcal{M}^\mathds{N}$ is computed from $PRNG_1$, where $\mathcal{M}\subset \mathds{N}^*$ is a finite nonempty set of integers.
633 The basic design procedure of the New CI generator is summarized in Algorithm~\ref{Chaotic iteration1}.
634 The internal state is $x$, the output state is $r$. $a$ and $b$ are those computed by the two input
635 PRNGs. Lastly, the value $g(a)$ is an integer defined as in Eq.~\ref{Formula}.
636 This function must be chosen such that the outputs of the resulted PRNG are uniform in $\llbracket 0, 2^\mathsf{N}-1 \rrbracket$. Function of \eqref{Formula} achieves this
637 goal (other candidates and more information can be found in ~\cite{bg10:ip}).
644 0 \text{ if }0 \leqslant{y^n}<{C^0_{32}},\\
645 1 \text{ if }{C^0_{32}} \leqslant{y^n}<\sum_{i=0}^1{C^i_{32}},\\
646 2 \text{ if }\sum_{i=0}^1{C^i_{32}} \leqslant{y^n}<\sum_{i=0}^2{C^i_{32}},\\
647 \vdots~~~~~ ~~\vdots~~~ ~~~~\\
648 N \text{ if }\sum_{i=0}^{N-1}{C^i_{32}}\leqslant{y^n}<1.\\
654 \textbf{Input:} the internal state $x$ (32 bits)\\
655 \textbf{Output:} a state $r$ of 32 bits
656 \begin{algorithmic}[1]
659 \STATE$d_i\leftarrow{0}$\;
662 \STATE$a\leftarrow{PRNG_1()}$\;
663 \STATE$k\leftarrow{g(a)}$\;
664 \WHILE{$i=0,\dots,k$}
666 \STATE$b\leftarrow{PRNG_2()~mod~\mathsf{N}}$\;
667 \STATE$S\leftarrow{b}$\;
670 \STATE $x_S\leftarrow{ \overline{x_S}}$\;
671 \STATE $d_S\leftarrow{1}$\;
676 \STATE $k\leftarrow{ k+1}$\;
679 \STATE $r\leftarrow{x}$\;
682 \caption{An arbitrary round of the new CI generator}
683 \label{Chaotic iteration1}
688 We have shown in~\cite{bfg12a:ip} that the use of chaotic iterations
689 implies an improvement of the statistical properties for all the
690 inputted defective generators we have investigated.
691 For instance, when considering the TestU01 battery with its 588 tests, we obtained 261
692 failures for a PRNG based on the logistic map alone, and
693 this number of failures falls below 138 in the Old CI(Logistic,Logistic) generator.
694 In the XORshift case (146 failures when considering it alone), the results are more amazing,
695 as the chaotic iterations post-treatment makes it fails only 8 tests.
696 Further investigations have been systematically realized in \cite{bfg12a:ip}
697 using a large set of inputted defective PRNGs, the three most used batteries of
698 tests (DieHARD, NIST, and TestU01), and for all the versions of generators we have proposed.
699 In all situations, an obvious improvement of the statistical behavior has
700 been obtained, reinforcing the impression that chaos leads to statistical
701 enhancement~\cite{bfg12a:ip}.
703 \subsection{Improving the Speed of the Former Generator}
705 Instead of updating only one cell at each iteration, we now propose to choose a
706 subset of components and to update them together, for speed improvement. Such a proposition leads
707 to a kind of merger of the two sequences used in Algorithms
708 \ref{CI Algorithm} and \ref{Chaotic iteration1}. When the updating function is the vectorial negation,
709 this algorithm can be rewritten as follows:
714 x^0 \in \llbracket 0, 2^\mathsf{N}-1 \rrbracket, S \in \llbracket 0, 2^\mathsf{N}-1 \rrbracket^\mathds{N} \\
715 \forall n \in \mathds{N}^*, x^n = x^{n-1} \oplus S^n,
718 \label{equation Oplus}
720 where $\oplus$ is for the bitwise exclusive or between two integers.
721 This rewriting can be understood as follows. The $n-$th term $S^n$ of the
722 sequence $S$, which is an integer of $\mathsf{N}$ binary digits, presents
723 the list of cells to update in the state $x^n$ of the system (represented
724 as an integer having $\mathsf{N}$ bits too). More precisely, the $k-$th
725 component of this state (a binary digit) changes if and only if the $k-$th
726 digit in the binary decomposition of $S^n$ is 1.
728 The single basic component presented in Eq.~\ref{equation Oplus} is of
729 ordinary use as a good elementary brick in various PRNGs. It corresponds
730 to the following discrete dynamical system in chaotic iterations:
733 \forall n\in \mathds{N}^{\ast }, \forall i\in
734 \llbracket1;\mathsf{N}\rrbracket ,x_i^n=\left\{
736 x_i^{n-1} & \text{ if } i \notin \mathcal{S}^n \\
737 \left(f(x^{n-1})\right)_{S^n} & \text{ if }i \in \mathcal{S}^n.
741 where $f$ is the vectorial negation and $\forall n \in \mathds{N}$,
742 $\mathcal{S}^n \subset \llbracket 1, \mathsf{N} \rrbracket$ is such that
743 $k \in \mathcal{S}^n$ if and only if the $k-$th digit in the binary
744 decomposition of $S^n$ is 1. Such chaotic iterations are more general
745 than the ones presented in Definition \ref{Def:chaotic iterations} because, instead of updating only one term at each iteration,
746 we select a subset of components to change.
749 Obviously, replacing the previous CI PRNG Algorithms by
750 Equation~\ref{equation Oplus}, which is possible when the iteration function is
751 the vectorial negation, leads to a speed improvement
752 (the resulting generator will be referred as ``Xor CI PRNG''
755 of chaos obtained in~\cite{bg10:ij} have been established
756 only for chaotic iterations of the form presented in Definition
757 \ref{Def:chaotic iterations}. The question is now to determine whether the
758 use of more general chaotic iterations to generate pseudorandom numbers
759 faster, does not deflate their topological chaos properties.
761 \subsection{Proofs of Chaos of the General Formulation of the Chaotic Iterations}
763 Let us consider the discrete dynamical systems in chaotic iterations having
764 the general form: $\forall n\in \mathds{N}^{\ast }$, $ \forall i\in
765 \llbracket1;\mathsf{N}\rrbracket $,
770 x_i^{n-1} & \text{ if } i \notin \mathcal{S}^n \\
771 \left(f(x^{n-1})\right)_{S^n} & \text{ if }i \in \mathcal{S}^n.
776 In other words, at the $n^{th}$ iteration, only the cells whose id is
777 contained into the set $S^{n}$ are iterated.
779 Let us now rewrite these general chaotic iterations as usual discrete dynamical
780 system of the form $X^{n+1}=f(X^n)$ on an ad hoc metric space. Such a formulation
781 is required in order to study the topological behavior of the system.
783 Let us introduce the following function:
786 \chi: & \llbracket 1; \mathsf{N} \rrbracket \times \mathcal{P}\left(\llbracket 1; \mathsf{N} \rrbracket\right) & \longrightarrow & \mathds{B}\\
787 & (i,X) & \longmapsto & \left\{ \begin{array}{ll} 0 & \textrm{if }i \notin X, \\ 1 & \textrm{if }i \in X, \end{array}\right.
790 where $\mathcal{P}\left(X\right)$ is for the powerset of the set $X$, that is, $Y \in \mathcal{P}\left(X\right) \Longleftrightarrow Y \subset X$.
792 Given a function $f:\mathds{B}^\mathsf{N} \longrightarrow \mathds{B}^\mathsf{N} $, define the function:
793 $F_{f}: \mathcal{P}\left(\llbracket1;\mathsf{N}\rrbracket \right) \times \mathds{B}^{\mathsf{N}}
794 \longrightarrow \mathds{B}^{\mathsf{N}}$
797 (P,E) & \longmapsto & \left( E_{j}.\chi (j,P)+f(E)_{j}.\overline{\chi(j,P)}\right) _{j\in \llbracket1;\mathsf{N}\rrbracket}%
800 where + and . are the Boolean addition and product operations, and $\overline{x}$
801 is the negation of the Boolean $x$.
802 Consider the phase space:
804 \mathcal{X} = \mathcal{P}\left(\llbracket 1 ; \mathsf{N} \rrbracket\right)^\mathds{N} \times
805 \mathds{B}^\mathsf{N},
807 \noindent and the map defined on $\mathcal{X}$:
809 G_f\left(S,E\right) = \left(\sigma(S), F_f(i(S),E)\right), %\label{Gf} %%RAPH, j'ai viré ce label qui existe déjà avant...
811 \noindent where $\sigma$ is the \emph{shift} function defined by $\sigma
812 (S^{n})_{n\in \mathds{N}}\in \mathcal{P}\left(\llbracket 1 ; \mathsf{N} \rrbracket\right)^\mathds{N}\longrightarrow (S^{n+1})_{n\in
813 \mathds{N}}\in \mathcal{P}\left(\llbracket 1 ; \mathsf{N} \rrbracket\right)^\mathds{N}$ and $i$ is the \emph{initial function}
814 $i:(S^{n})_{n\in \mathds{N}} \in \mathcal{P}\left(\llbracket 1 ; \mathsf{N} \rrbracket\right)^\mathds{N}\longrightarrow S^{0}\in \mathcal{P}\left(\llbracket 1 ; \mathsf{N} \rrbracket\right)$.
815 Then the general chaotic iterations defined in Equation \ref{general CIs} can
816 be described by the following discrete dynamical system:
820 X^0 \in \mathcal{X} \\
826 Once more, a shift function appears as a component of these general chaotic
829 To study the Devaney's chaos property, a distance between two points
830 $X = (S,E), Y = (\check{S},\check{E})$ of $\mathcal{X}$ must be defined.
833 d(X,Y)=d_{e}(E,\check{E})+d_{s}(S,\check{S}),
836 \noindent where $ \displaystyle{d_{e}(E,\check{E})} = \displaystyle{\sum_{k=1}^{\mathsf{N}%
837 }\delta (E_{k},\check{E}_{k})}$ is once more the Hamming distance, and
838 $ \displaystyle{d_{s}(S,\check{S})} = \displaystyle{\dfrac{9}{\mathsf{N}}%
839 \sum_{k=1}^{\infty }\dfrac{|S^k\Delta {S}^k|}{10^{k}}}$,
840 %%RAPH : ici, j'ai supprimé tous les sauts à la ligne
843 %% \begin{array}{lll}
844 %% \displaystyle{d_{e}(E,\check{E})} & = & \displaystyle{\sum_{k=1}^{\mathsf{N}%
845 %% }\delta (E_{k},\check{E}_{k})} \textrm{ is once more the Hamming distance}, \\
846 %% \displaystyle{d_{s}(S,\check{S})} & = & \displaystyle{\dfrac{9}{\mathsf{N}}%
847 %% \sum_{k=1}^{\infty }\dfrac{|S^k\Delta {S}^k|}{10^{k}}}.%
851 where $|X|$ is the cardinality of a set $X$ and $A\Delta B$ is for the symmetric difference, defined for sets A, B as
852 $A\,\Delta\,B = (A \setminus B) \cup (B \setminus A)$.
856 The function $d$ defined in Eq.~\ref{nouveau d} is a metric on $\mathcal{X}$.
860 $d_e$ is the Hamming distance. We will prove that $d_s$ is a distance
861 too, thus $d$, as being the sum of two distances, will also be a distance.
863 \item Obviously, $d_s(S,\check{S})\geqslant 0$, and if $S=\check{S}$, then
864 $d_s(S,\check{S})=0$. Conversely, if $d_s(S,\check{S})=0$, then
865 $\forall k \in \mathds{N}, |S^k\Delta {S}^k|=0$, and so $\forall k, S^k=\check{S}^k$.
866 \item $d_s$ is symmetric
867 ($d_s(S,\check{S})=d_s(\check{S},S)$) due to the commutative property
868 of the symmetric difference.
869 \item Finally, $|S \Delta S''| = |(S \Delta \varnothing) \Delta S''|= |S \Delta (S'\Delta S') \Delta S''|= |(S \Delta S') \Delta (S' \Delta S'')|\leqslant |S \Delta S'| + |S' \Delta S''|$,
870 and so for all subsets $S,S',$ and $S''$ of $\llbracket 1, \mathsf{N} \rrbracket$,
871 we have $d_s(S,S'') \leqslant d_e(S,S')+d_s(S',S'')$, and the triangle
872 inequality is obtained.
877 Before being able to study the topological behavior of the general
878 chaotic iterations, we must first establish that:
881 For all $f:\mathds{B}^\mathsf{N} \longrightarrow \mathds{B}^\mathsf{N} $, the function $G_f$ is continuous on
882 $\left( \mathcal{X},d\right)$.
887 We use the sequential continuity.
888 Let $(S^n,E^n)_{n\in \mathds{N}}$ be a sequence of the phase space $%
889 \mathcal{X}$, which converges to $(S,E)$. We will prove that $\left(
890 G_{f}(S^n,E^n)\right) _{n\in \mathds{N}}$ converges to $\left(
891 G_{f}(S,E)\right) $. Let us remark that for all $n$, $S^n$ is a strategy,
892 thus, we consider a sequence of strategies (\emph{i.e.}, a sequence of
894 As $d((S^n,E^n);(S,E))$ converges to 0, each distance $d_{e}(E^n,E)$ and $d_{s}(S^n,S)$ converges
895 to 0. But $d_{e}(E^n,E)$ is an integer, so $\exists n_{0}\in \mathds{N},$ $%
896 d_{e}(E^n,E)=0$ for any $n\geqslant n_{0}$.\newline
897 In other words, there exists a threshold $n_{0}\in \mathds{N}$ after which no
898 cell will change its state:
899 $\exists n_{0}\in \mathds{N},n\geqslant n_{0}\Rightarrow E^n = E.$
901 In addition, $d_{s}(S^n,S)\longrightarrow 0,$ so $\exists n_{1}\in %
902 \mathds{N},d_{s}(S^n,S)<10^{-1}$ for all indexes greater than or equal to $%
903 n_{1}$. This means that for $n\geqslant n_{1}$, all the $S^n$ have the same
904 first term, which is $S^0$: $\forall n\geqslant n_{1},S_0^n=S_0.$
906 Thus, after the $max(n_{0},n_{1})^{th}$ term, states of $E^n$ and $E$ are
907 identical and strategies $S^n$ and $S$ start with the same first term.\newline
908 Consequently, states of $G_{f}(S^n,E^n)$ and $G_{f}(S,E)$ are equal,
909 so, after the $max(n_0, n_1)^{th}$ term, the distance $d$ between these two points is strictly less than 1.\newline
910 \noindent We now prove that the distance between $\left(
911 G_{f}(S^n,E^n)\right) $ and $\left( G_{f}(S,E)\right) $ is convergent to
912 0. Let $\varepsilon >0$. \medskip
914 \item If $\varepsilon \geqslant 1$, we see that the distance
915 between $\left( G_{f}(S^n,E^n)\right) $ and $\left( G_{f}(S,E)\right) $ is
916 strictly less than 1 after the $max(n_{0},n_{1})^{th}$ term (same state).
918 \item If $\varepsilon <1$, then $\exists k\in \mathds{N},10^{-k}\geqslant
919 \varepsilon > 10^{-(k+1)}$. But $d_{s}(S^n,S)$ converges to 0, so
921 \exists n_{2}\in \mathds{N},\forall n\geqslant
922 n_{2},d_{s}(S^n,S)<10^{-(k+2)},
924 thus after $n_{2}$, the $k+2$ first terms of $S^n$ and $S$ are equal.
926 \noindent As a consequence, the $k+1$ first entries of the strategies of $%
927 G_{f}(S^n,E^n)$ and $G_{f}(S,E)$ are the same ($G_{f}$ is a shift of strategies) and due to the definition of $d_{s}$, the floating part of
928 the distance between $(S^n,E^n)$ and $(S,E)$ is strictly less than $%
929 10^{-(k+1)}\leqslant \varepsilon $.
932 %%RAPH : ici j'ai rajouté une ligne
933 %%TOF : ici j'ai rajouté un commentaire
936 \forall \varepsilon >0,$ $\exists N_{0}=max(n_{0},n_{1},n_{2})\in \mathds{N}
937 ,$ $\forall n\geqslant N_{0},$
938 $ d\left( G_{f}(S^n,E^n);G_{f}(S,E)\right)
939 \leqslant \varepsilon .
941 $G_{f}$ is consequently continuous.
945 It is now possible to study the topological behavior of the general chaotic
946 iterations. We will prove that,
949 \label{t:chaos des general}
950 The general chaotic iterations defined on Equation~\ref{general CIs} satisfy
951 the Devaney's property of chaos.
954 Let us firstly prove the following lemma.
956 \begin{lemma}[Strong transitivity]
958 For all couples $X,Y \in \mathcal{X}$ and any neighborhood $V$ of $X$, we can
959 find $n \in \mathds{N}^*$ and $X' \in V$ such that $G^n(X')=Y$.
963 Let $X=(S,E)$, $\varepsilon>0$, and $k_0 = \lfloor log_{10}(\varepsilon)+1 \rfloor$.
964 Any point $X'=(S',E')$ such that $E'=E$ and $\forall k \leqslant k_0, S'^k=S^k$,
965 are in the open ball $\mathcal{B}\left(X,\varepsilon\right)$. Let us define
966 $\check{X} = \left(\check{S},\check{E}\right)$, where $\check{X}= G^{k_0}(X)$.
967 We denote by $s\subset \llbracket 1; \mathsf{N} \rrbracket$ the set of coordinates
968 that are different between $\check{E}$ and the state of $Y$. Thus each point $X'$ of
969 the form $(S',E')$ where $E'=E$ and $S'$ starts with
970 $(S^0, S^1, \hdots, S^{k_0},s,\hdots)$, verifies the following properties:
972 \item $X'$ is in $\mathcal{B}\left(X,\varepsilon\right)$,
973 \item the state of $G_f^{k_0+1}(X')$ is the state of $Y$.
975 Finally the point $\left(\left(S^0, S^1, \hdots, S^{k_0},s,s^0, s^1, \hdots\right); E\right)$,
976 where $(s^0,s^1, \hdots)$ is the strategy of $Y$, satisfies the properties
977 claimed in the lemma.
980 We can now prove the Theorem~\ref{t:chaos des general}.
982 \begin{proof}[Theorem~\ref{t:chaos des general}]
983 Firstly, strong transitivity implies transitivity.
985 Let $(S,E) \in\mathcal{X}$ and $\varepsilon >0$. To
986 prove that $G_f$ is regular, it is sufficient to prove that
987 there exists a strategy $\tilde S$ such that the distance between
988 $(\tilde S,E)$ and $(S,E)$ is less than $\varepsilon$, and such that
989 $(\tilde S,E)$ is a periodic point.
991 Let $t_1=\lfloor-\log_{10}(\varepsilon)\rfloor$, and let $E'$ be the
992 configuration that we obtain from $(S,E)$ after $t_1$ iterations of
993 $G_f$. As $G_f$ is strongly transitive, there exists a strategy $S'$
994 and $t_2\in\mathds{N}$ such
995 that $E$ is reached from $(S',E')$ after $t_2$ iterations of $G_f$.
997 Consider the strategy $\tilde S$ that alternates the first $t_1$ terms
998 of $S$ and the first $t_2$ terms of $S'$:
999 %%RAPH : j'ai coupé la ligne en 2
1001 S=(S_0,\dots,S_{t_1-1},S'_0,\dots,S'_{t_2-1},S_0,$$$$\dots,S_{t_1-1},S'_0,\dots,S'_{t_2-1},S_0,\dots).$$ It
1002 is clear that $(\tilde S,E)$ is obtained from $(\tilde S,E)$ after
1003 $t_1+t_2$ iterations of $G_f$. So $(\tilde S,E)$ is a periodic
1004 point. Since $\tilde S_t=S_t$ for $t<t_1$, by the choice of $t_1$, we
1005 have $d((S,E),(\tilde S,E))<\epsilon$.
1009 %\section{Statistical Improvements Using Chaotic Iterations}
1011 %\label{The generation of pseudorandom sequence}
1014 %Let us now explain why we have reasonable ground to believe that chaos
1015 %can improve statistical properties.
1016 %We will show in this section that chaotic properties as defined in the
1017 %mathematical theory of chaos are related to some statistical tests that can be found
1018 %in the NIST battery. Furthermore, we will check that, when mixing defective PRNGs with
1019 %chaotic iterations, the new generator presents better statistical properties
1020 %(this section summarizes and extends the work of~\cite{bfg12a:ip}).
1024 %\subsection{Qualitative relations between topological properties and statistical tests}
1027 %There are various relations between topological properties that describe an unpredictable behavior for a discrete
1028 %dynamical system on the one
1029 %hand, and statistical tests to check the randomness of a numerical sequence
1030 %on the other hand. These two mathematical disciplines follow a similar
1031 %objective in case of a recurrent sequence (to characterize an intrinsically complicated behavior for a
1032 %recurrent sequence), with two different but complementary approaches.
1033 %It is true that the following illustrative links give only qualitative arguments,
1034 %and proofs should be provided later to make such arguments irrefutable. However
1035 %they give a first understanding of the reason why we think that chaotic properties should tend
1036 %to improve the statistical quality of PRNGs.
1038 %Let us now list some of these relations between topological properties defined in the mathematical
1039 %theory of chaos and tests embedded into the NIST battery. %Such relations need to be further
1040 %%investigated, but they presently give a first illustration of a trend to search similar properties in the
1041 %%two following fields: mathematical chaos and statistics.
1045 % \item \textbf{Regularity}. As stated in Section~\ref{subsec:Devaney}, a chaotic dynamical system must
1046 %have an element of regularity. Depending on the chosen definition of chaos, this element can be the existence of
1047 %a dense orbit, the density of periodic points, etc. The key idea is that a dynamical system with no periodicity
1048 %is not as chaotic as a system having periodic orbits: in the first situation, we can predict something and gain a
1049 %knowledge about the behavior of the system, that is, it never enters into a loop. A similar importance for periodicity is emphasized in
1050 %the two following NIST tests~\cite{Nist10}:
1052 % \item \textbf{Non-overlapping Template Matching Test}. Detect generators that produce too many occurrences of a given non-periodic (aperiodic) pattern.
1053 % \item \textbf{Discrete Fourier Transform (Spectral) Test}. Detect periodic features (i.e., repetitive patterns that are close one to another) in the tested sequence that would indicate a deviation from the assumption of randomness.
1056 %\item \textbf{Transitivity}. This topological property previously introduced states that the dynamical system is intrinsically complicated: it cannot be simplified into
1057 %two subsystems that do not interact, as we can find in any neighborhood of any point another point whose orbit visits the whole phase space.
1058 %This focus on the places visited by the orbits of the dynamical system takes various nonequivalent formulations in the mathematical theory
1059 %of chaos, namely: transitivity, strong transitivity, total transitivity, topological mixing, and so on~\cite{bg10:ij}. A similar attention
1060 %is brought on the states visited during a random walk in the two tests below~\cite{Nist10}:
1062 % \item \textbf{Random Excursions Variant Test}. Detect deviations from the expected number of visits to various states in the random walk.
1063 % \item \textbf{Random Excursions Test}. Determine if the number of visits to a particular state within a cycle deviates from what one would expect for a random sequence.
1066 %\item \textbf{Chaos according to Li and Yorke}. Two points of the phase space $(x,y)$ define a couple of Li-Yorke when $\limsup_{n \rightarrow +\infty} d(f^{(n)}(x), f^{(n)}(y))>0$ et $\liminf_{n \rightarrow +\infty} d(f^{(n)}(x), f^{(n)}(y))=0$, meaning that their orbits always oscillate as the iterations pass. When a system is compact and contains an uncountable set of such points, it is claimed as chaotic according
1067 %to Li-Yorke~\cite{Li75,Ruette2001}. A similar property is regarded in the following NIST test~\cite{Nist10}.
1069 % \item \textbf{Runs Test}. To determine whether the number of runs of ones and zeros of various lengths is as expected for a random sequence. In particular, this test determines whether the oscillation between such zeros and ones is too fast or too slow.
1071 % \item \textbf{Topological entropy}. The desire to formulate an equivalency of the thermodynamics entropy
1072 %has emerged both in the topological and statistical fields. Once again, a similar objective has led to two different
1073 %rewritting of an entropy based disorder: the famous Shannon definition of entropy is approximated in the statistical approach,
1074 %whereas topological entropy is defined as follows:
1075 %$x,y \in \mathcal{X}$ are $\varepsilon-$\emph{separated in time $n$} if there exists $k \leqslant n$ such that $d\left(f^{(k)}(x),f^{(k)}(y)\right)>\varepsilon$. Then $(n,\varepsilon)-$separated sets are sets of points that are all $\varepsilon-$separated in time $n$, which
1076 %leads to the definition of $s_n(\varepsilon,Y)$, being the maximal cardinality of all $(n,\varepsilon)-$separated sets. Using these notations,
1077 %the topological entropy is defined as follows: $$h_{top}(\mathcal{X},f) = \displaystyle{\lim_{\varepsilon \rightarrow 0} \Big[ \limsup_{n \rightarrow +\infty} \dfrac{1}{n} \log s_n(\varepsilon,\mathcal{X})\Big]}.$$
1078 %This value measures the average exponential growth of the number of distinguishable orbit segments.
1079 %In this sense, it measures the complexity of the topological dynamical system, whereas
1080 %the Shannon approach comes to mind when defining the following test~\cite{Nist10}:
1082 %\item \textbf{Approximate Entropy Test}. Compare the frequency of the overlapping blocks of two consecutive/adjacent lengths ($m$ and $m+1$) against the expected result for a random sequence.
1085 % \item \textbf{Non-linearity, complexity}. Finally, let us remark that non-linearity and complexity are
1086 %not only sought in general to obtain chaos, but they are also required for randomness, as illustrated by the two tests below~\cite{Nist10}.
1088 %\item \textbf{Binary Matrix Rank Test}. Check for linear dependence among fixed length substrings of the original sequence.
1089 %\item \textbf{Linear Complexity Test}. Determine whether or not the sequence is complex enough to be considered random.
1094 %We have proven in our previous works~\cite{guyeux12:bc} that chaotic iterations satisfying Theorem~\ref{Th:Caractérisation des IC chaotiques} are, among other
1095 %things, strongly transitive, topologically mixing, chaotic as defined by Li and Yorke,
1096 %and that they have a topological entropy and an exponent of Lyapunov both equal to $ln(\mathsf{N})$,
1097 %where $\mathsf{N}$ is the size of the iterated vector.
1098 %These topological properties make that we are ground to believe that a generator based on chaotic
1099 %iterations will probably be able to pass all the existing statistical batteries for pseudorandomness like
1100 %the NIST one. The following subsections, in which we prove that defective generators have their
1101 %statistical properties improved by chaotic iterations, show that such an assumption is true.
1103 %\subsection{Details of some Existing Generators}
1105 %The list of defective PRNGs we will use
1106 %as inputs for the statistical tests to come is introduced here.
1108 %Firstly, the simple linear congruency generators (LCGs) will be used.
1109 %They are defined by the following recurrence:
1111 %x^n = (ax^{n-1} + c)~mod~m,
1114 %where $a$, $c$, and $x^0$ must be, among other things, non-negative and inferior to
1115 %$m$~\cite{LEcuyerS07}. In what follows, 2LCGs and 3LCGs refer to two (resp. three)
1116 %combinations of such LCGs. For further details, see~\cite{bfg12a:ip,combined_lcg}.
1118 %Secondly, the multiple recursive generators (MRGs) which will be used,
1119 %are based on a linear recurrence of order
1120 %$k$, modulo $m$~\cite{LEcuyerS07}:
1122 %x^n = (a^1x^{n-1}+~...~+a^kx^{n-k})~mod~m .
1125 %The combination of two MRGs (referred as 2MRGs) is also used in these experiments.
1127 %Generators based on linear recurrences with carry will be regarded too.
1128 %This family of generators includes the add-with-carry (AWC) generator, based on the recurrence:
1132 %x^n = (x^{n-r} + x^{n-s} + c^{n-1})~mod~m, \\
1133 %c^n= (x^{n-r} + x^{n-s} + c^{n-1}) / m, \end{array}\end{equation}
1134 %the SWB generator, having the recurrence:
1138 %x^n = (x^{n-r} - x^{n-s} - c^{n-1})~mod~m, \\
1141 %1 ~~~~~\text{if}~ (x^{i-r} - x^{i-s} - c^{i-1})<0\\
1142 %0 ~~~~~\text{else},\end{array} \right. \end{array}\end{equation}
1143 %and the SWC generator, which is based on the following recurrence:
1147 %x^n = (a^1x^{n-1} \oplus ~...~ \oplus a^rx^{n-r} \oplus c^{n-1}) ~ mod ~ 2^w, \\
1148 %c^n = (a^1x^{n-1} \oplus ~...~ \oplus a^rx^{n-r} \oplus c^{n-1}) ~ / ~ 2^w. \end{array}\end{equation}
1150 %Then the generalized feedback shift register (GFSR) generator has been implemented, that is:
1152 %x^n = x^{n-r} \oplus x^{n-k} .
1157 %Finally, the nonlinear inversive (INV) generator~\cite{LEcuyerS07} has been studied, which is:
1164 %(a^1 + a^2 / z^{n-1})~mod~m & \text{if}~ z^{n-1} \neq 0 \\
1165 %a^1 & \text{if}~ z^{n-1} = 0 .\end{array} \right. \end{array}\end{equation}
1170 %%\renewcommand{\arraystretch}{1}
1171 %\caption{TestU01 Statistical Test Failures}
1174 % \begin{tabular}{lccccc}
1176 %Test name &Tests& Logistic & XORshift & ISAAC\\
1177 %Rabbit & 38 &21 &14 &0 \\
1178 %Alphabit & 17 &16 &9 &0 \\
1179 %Pseudo DieHARD &126 &0 &2 &0 \\
1180 %FIPS\_140\_2 &16 &0 &0 &0 \\
1181 %SmallCrush &15 &4 &5 &0 \\
1182 %Crush &144 &95 &57 &0 \\
1183 %Big Crush &160 &125 &55 &0 \\ \hline
1184 %Failures & &261 &146 &0 \\
1192 %%\renewcommand{\arraystretch}{1}
1193 %\caption{TestU01 Statistical Test Failures for Old CI algorithms ($\mathsf{N}=4$)}
1194 %\label{TestU01 for Old CI}
1196 % \begin{tabular}{lcccc}
1198 %\multirow{3}*{Test name} & \multicolumn{4}{c}{Old CI}\\
1199 %&Logistic& XORshift& ISAAC&ISAAC \\
1201 %&Logistic& XORshift& XORshift&ISAAC \\ \cmidrule(r){2-5}
1202 %Rabbit &7 &2 &0 &0 \\
1203 %Alphabit & 3 &0 &0 &0 \\
1204 %DieHARD &0 &0 &0 &0 \\
1205 %FIPS\_140\_2 &0 &0 &0 &0 \\
1206 %SmallCrush &2 &0 &0 &0 \\
1207 %Crush &47 &4 &0 &0 \\
1208 %Big Crush &79 &3 &0 &0 \\ \hline
1209 %Failures &138 &9 &0 &0 \\
1218 %\subsection{Statistical tests}
1219 %\label{Security analysis}
1221 %Three batteries of tests are reputed and regularly used
1222 %to evaluate the statistical properties of newly designed pseudorandom
1223 %number generators. These batteries are named DieHard~\cite{Marsaglia1996},
1224 %the NIST suite~\cite{ANDREW2008}, and the most stringent one called
1225 %TestU01~\cite{LEcuyerS07}, which encompasses the two other batteries.
1229 %\label{Results and discussion}
1231 %%\renewcommand{\arraystretch}{1}
1232 %\caption{NIST and DieHARD tests suite passing rates for PRNGs without CI}
1233 %\label{NIST and DieHARD tests suite passing rate the for PRNGs without CI}
1235 % \begin{tabular}{|l||c|c|c|c|c|c|c|c|c|c|}
1237 %Types of PRNGs & \multicolumn{2}{c|}{Linear PRNGs} & \multicolumn{4}{c|}{Lagged PRNGs} & \multicolumn{1}{c|}{ICG PRNGs} & \multicolumn{3}{c|}{Mixed PRNGs}\\ \hline
1238 %\backslashbox{\textbf{$Tests$}} {\textbf{$PRNG$}} & LCG& MRG& AWC & SWB & SWC & GFSR & INV & LCG2& LCG3& MRG2 \\ \hline
1239 %NIST & 11/15 & 14/15 &\textbf{15/15} & \textbf{15/15} & 14/15 & 14/15 & 14/15 & 14/15& 14/15& 14/15 \\ \hline
1240 %DieHARD & 16/18 & 16/18 & 15/18 & 16/18 & \textbf{18/18} & 16/18 & 16/18 & 16/18& 16/18& 16/18\\ \hline
1244 %Table~\ref{NIST and DieHARD tests suite passing rate the for PRNGs without CI} shows the
1245 %results on the two first batteries recalled above, indicating that all the PRNGs presented
1246 %in the previous section
1247 %cannot pass all these tests. In other words, the statistical quality of these PRNGs cannot
1248 %fulfill the up-to-date standards presented previously. We have shown in~\cite{bfg12a:ip} that the use of chaotic
1249 %iterations can solve this issue.
1250 %%More precisely, to
1251 %%illustrate the effects of chaotic iterations on these defective PRNGs, experiments have been divided in three parts~\cite{bfg12a:ip}:
1253 %% \item \textbf{Single CIPRNG}: The PRNGs involved in CI computing are of the same category.
1254 %% \item \textbf{Mixed CIPRNG}: Two different types of PRNGs are mixed during the chaotic iterations process.
1255 %% \item \textbf{Multiple CIPRNG}: The generator is obtained by repeating the composition of the iteration function as follows: $x^0\in \mathds{B}^{\mathsf{N}}$, and $\forall n\in \mathds{N}^{\ast },\forall i\in \llbracket1;\mathsf{N}\rrbracket, x_i^n=$
1260 %%x_i^{n-1}~~~~~\text{if}~S^n\neq i \\
1261 %%\forall j\in \llbracket1;\mathsf{m}\rrbracket,f^m(x^{n-1})_{S^{nm+j}}~\text{if}~S^{nm+j}=i.\end{array} \right. \end{array}
1263 %%$m$ is called the \emph{functional power}.
1266 %The obtained results are reproduced in Table
1267 %\ref{NIST and DieHARD tests suite passing rate the for single CIPRNGs}.
1268 %The scores written in boldface indicate that all the tests have been passed successfully, whereas an
1269 %asterisk ``*'' means that the considered passing rate has been improved.
1270 %The improvements are obvious for both the ``Old CI'' and the ``New CI'' generators.
1271 %Concerning the ``Xor CI PRNG'', the score is less spectacular. Because of a large speed improvement, the statistics
1272 % are not as good as for the two other versions of these CIPRNGs.
1273 %However 8 tests have been improved (with no deflation for the other results).
1277 %%\renewcommand{\arraystretch}{1.3}
1278 %\caption{NIST and DieHARD tests suite passing rates for PRNGs with CI}
1279 %\label{NIST and DieHARD tests suite passing rate the for single CIPRNGs}
1281 % \begin{tabular}{|l||c|c|c|c|c|c|c|c|c|c|c|c|}
1283 %Types of PRNGs & \multicolumn{2}{c|}{Linear PRNGs} & \multicolumn{4}{c|}{Lagged PRNGs} & \multicolumn{1}{c|}{ICG PRNGs} & \multicolumn{3}{c|}{Mixed PRNGs}\\ \hline
1284 %\backslashbox{\textbf{$Tests$}} {\textbf{$Single~CIPRNG$}} & LCG & MRG & AWC & SWB & SWC & GFSR & INV& LCG2 & LCG3& MRG2 \\ \hline\hline
1285 %Old CIPRNG\\ \hline \hline
1286 %NIST & \textbf{15/15} * & \textbf{15/15} * & \textbf{15/15} & \textbf{15/15} & \textbf{15/15} * & \textbf{15/15} * & \textbf{15/15} *& \textbf{15/15} * & \textbf{15/15} * & \textbf{15/15} \\ \hline
1287 %DieHARD & \textbf{18/18} * & \textbf{18/18} * & \textbf{18/18} * & \textbf{18/18} * & \textbf{18/18} & \textbf{18/18} * & \textbf{18/18} *& \textbf{18/18} * & \textbf{18/18} *& \textbf{18/18} * \\ \hline
1288 %New CIPRNG\\ \hline \hline
1289 %NIST & \textbf{15/15} * & \textbf{15/15} * & \textbf{15/15} & \textbf{15/15} & \textbf{15/15} * & \textbf{15/15} * & \textbf{15/15} *& \textbf{15/15} * & \textbf{15/15} * & \textbf{15/15} \\ \hline
1290 %DieHARD & \textbf{18/18} * & \textbf{18/18} * & \textbf{18/18} * & \textbf{18/18} * & \textbf{18/18} & \textbf{18/18} * & \textbf{18/18} * & \textbf{18/18} * & \textbf{18/18} *& \textbf{18/18} *\\ \hline
1291 %Xor CIPRNG\\ \hline\hline
1292 %NIST & 14/15*& \textbf{15/15} * & \textbf{15/15} & \textbf{15/15} & 14/15 & \textbf{15/15} * & 14/15& \textbf{15/15} * & \textbf{15/15} *& \textbf{15/15} \\ \hline
1293 %DieHARD & 16/18 & 16/18 & 17/18* & \textbf{18/18} * & \textbf{18/18} & \textbf{18/18} * & 16/18 & 16/18 & 16/18& 16/18\\ \hline
1298 %We have then investigated in~\cite{bfg12a:ip} if it were possible to improve
1299 %the statistical behavior of the Xor CI version by combining more than one
1300 %$\oplus$ operation. Results are summarized in Table~\ref{threshold}, illustrating
1301 %the progressive increasing effects of chaotic iterations, when giving time to chaos to get settled in.
1302 %Thus rapid and perfect PRNGs, regarding the NIST and DieHARD batteries, can be obtained
1303 %using chaotic iterations on defective generators.
1306 %%\renewcommand{\arraystretch}{1.3}
1307 %\caption{Number of $\oplus$ operations to pass the whole NIST and DieHARD batteries}
1310 % \begin{tabular}{|l||c|c|c|c|c|c|c|c|}
1312 %Inputted $PRNG$ & LCG & MRG & SWC & GFSR & INV& LCG2 & LCG3 & MRG2 \\ \hline\hline
1313 %Threshold value $m$& 19 & 7 & 2& 1 & 11& 9& 3& 4\\ \hline\hline
1317 %Finally, the TestU01 battery has been launched on three well-known generators
1318 %(a logistic map, a simple XORshift, and the cryptographically secure ISAAC,
1319 %see Table~\ref{TestU011}). These results can be compared with
1320 %Table~\ref{TestU01 for Old CI}, which gives the scores obtained by the
1321 %Old CI PRNG that has received these generators.
1322 %The obvious improvement speaks for itself, and together with the other
1323 %results recalled in this section, it reinforces the opinion that a strong
1324 %correlation between topological properties and statistical behavior exists.
1327 %The next subsection will now give a concrete original implementation of the Xor CI PRNG, the
1328 %fastest generator in the chaotic iteration based family. In the remainder,
1329 %this generator will be simply referred to as CIPRNG, or ``the proposed PRNG'', if this statement does not
1333 \section{Toward Efficiency and Improvement for CI PRNG}
1334 \label{sec:efficient PRNG}
1336 \subsection{First Efficient Implementation of a PRNG based on Chaotic Iterations}
1338 %Based on the proof presented in the previous section, it is now possible to
1339 %improve the speed of the generator formerly presented in~\cite{bgw09:ip,guyeux10}.
1340 %The first idea is to consider
1341 %that the provided strategy is a pseudorandom Boolean vector obtained by a
1343 %An iteration of the system is simply the bitwise exclusive or between
1344 %the last computed state and the current strategy.
1345 %Topological properties of disorder exhibited by chaotic
1346 %iterations can be inherited by the inputted generator, we hope by doing so to
1347 %obtain some statistical improvements while preserving speed.
1349 %%RAPH : j'ai viré tout ca
1350 %% Let us give an example using 16-bits numbers, to clearly understand how the bitwise xor operations
1353 %% Suppose that $x$ and the strategy $S^i$ are given as
1355 %% Table~\ref{TableExemple} shows the result of $x \oplus S^i$.
1358 %% \begin{scriptsize}
1360 %% \begin{array}{|cc|cccccccccccccccc|}
1362 %% x &=&1&0&1&1&1&0&1&0&1&0&0&1&0&0&1&0\\
1364 %% S^i &=&0&1&1&0&0&1&1&0&1&1&1&0&0&1&1&1\\
1366 %% x \oplus S^i&=&1&1&0&1&1&1&0&0&0&1&1&1&0&1&0&1\\
1373 %% \caption{Example of an arbitrary round of the proposed generator}
1374 %% \label{TableExemple}
1380 \lstset{language=C,caption={C code of the sequential PRNG based on chaotic iterations},label={algo:seqCIPRNG}}
1384 unsigned int CIPRNG() {
1385 static unsigned int x = 123123123;
1386 unsigned long t1 = xorshift();
1387 unsigned long t2 = xor128();
1388 unsigned long t3 = xorwow();
1389 x = x^(unsigned int)t1;
1390 x = x^(unsigned int)(t2>>32);
1391 x = x^(unsigned int)(t3>>32);
1392 x = x^(unsigned int)t2;
1393 x = x^(unsigned int)(t1>>32);
1394 x = x^(unsigned int)t3;
1402 In Listing~\ref{algo:seqCIPRNG} a sequential version of the proposed PRNG based
1403 on chaotic iterations is presented. The xor operator is represented by
1404 \textasciicircum. This function uses three classical 64-bits PRNGs, namely the
1405 \texttt{xorshift}, the \texttt{xor128}, and the
1406 \texttt{xorwow}~\cite{Marsaglia2003}. In the following, we call them ``xor-like
1407 PRNGs''. As each xor-like PRNG uses 64-bits whereas our proposed generator
1408 works with 32-bits, we use the command \texttt{(unsigned int)}, that selects the
1409 32 least significant bits of a given integer, and the code \texttt{(unsigned
1410 int)(t$>>$32)} in order to obtain the 32 most significant bits of \texttt{t}.
1412 Thus producing a pseudorandom number needs 6 xor operations with 6 32-bits numbers
1413 that are provided by 3 64-bits PRNGs. This version successfully passes the
1414 stringent BigCrush battery of tests~\cite{LEcuyerS07}.
1415 At this point, we thus
1416 have defined an efficient and statistically unbiased generator. Its speed is
1417 directly related to the use of linear operations, but for the same reason,
1418 this fast generator cannot be proven as secure.
1422 \subsection{Efficient PRNGs based on Chaotic Iterations on GPU}
1423 \label{sec:efficient PRNG gpu}
1425 In order to take benefits from the computing power of GPU, a program
1426 needs to have independent blocks of threads that can be computed
1427 simultaneously. In general, the larger the number of threads is, the
1428 more local memory is used, and the less branching instructions are
1429 used (if, while, ...), the better the performances on GPU is.
1430 Obviously, having these requirements in mind, it is possible to build
1431 a program similar to the one presented in Listing
1432 \ref{algo:seqCIPRNG}, which computes pseudorandom numbers on GPU. To
1433 do so, we must firstly recall that in the CUDA~\cite{Nvid10}
1434 environment, threads have a local identifier called
1435 \texttt{ThreadIdx}, which is relative to the block containing
1436 them. Furthermore, in CUDA, parts of the code that are executed by the GPU, are
1437 called {\it kernels}.
1440 \subsection{Naive Version for GPU}
1443 It is possible to deduce from the CPU version a quite similar version adapted to GPU.
1444 The simple principle consists in making each thread of the GPU computing the CPU version of our PRNG.
1445 Of course, the three xor-like
1446 PRNGs used in these computations must have different parameters.
1447 In a given thread, these parameters are
1448 randomly picked from another PRNGs.
1449 The initialization stage is performed by the CPU.
1450 To do it, the ISAAC PRNG~\cite{Jenkins96} is used to set all the
1451 parameters embedded into each thread.
1453 The implementation of the three
1454 xor-like PRNGs is straightforward when their parameters have been
1455 allocated in the GPU memory. Each xor-like works with an internal
1456 number $x$ that saves the last generated pseudorandom number. Additionally, the
1457 implementation of the xor128, the xorshift, and the xorwow respectively require
1458 4, 5, and 6 unsigned long as internal variables.
1463 \KwIn{InternalVarXorLikeArray: array with internal variables of the 3 xor-like
1464 PRNGs in global memory\;
1465 NumThreads: number of threads\;}
1466 \KwOut{NewNb: array containing random numbers in global memory}
1467 \If{threadIdx is concerned by the computation} {
1468 retrieve data from InternalVarXorLikeArray[threadIdx] in local variables\;
1470 compute a new PRNG as in Listing\ref{algo:seqCIPRNG}\;
1471 store the new PRNG in NewNb[NumThreads*threadIdx+i]\;
1473 store internal variables in InternalVarXorLikeArray[threadIdx]\;
1476 \caption{Main kernel of the GPU ``naive'' version of the PRNG based on chaotic iterations}
1477 \label{algo:gpu_kernel}
1482 Algorithm~\ref{algo:gpu_kernel} presents a naive implementation of the proposed PRNG on
1483 GPU. Due to the available memory in the GPU and the number of threads
1484 used simultaneously, the number of random numbers that a thread can generate
1485 inside a kernel is limited (\emph{i.e.}, the variable \texttt{n} in
1486 algorithm~\ref{algo:gpu_kernel}). For instance, if $100,000$ threads are used and
1487 if $n=100$\footnote{in fact, we need to add the initial seed (a 32-bits number)},
1488 then the memory required to store all of the internals variables of both the xor-like
1489 PRNGs\footnote{we multiply this number by $2$ in order to count 32-bits numbers}
1490 and the pseudorandom numbers generated by our PRNG, is equal to $100,000\times ((4+5+6)\times
1491 2+(1+100))=1,310,000$ 32-bits numbers, that is, approximately $52$Mb.
1493 Remark that the only requirement regarding the seed regarding the security of our PRNG is
1494 that it must be randomly picked. Indeed, the asymptotic security of BBS guarantees
1495 that, as the seed length increases, no polynomial time statistical test can
1496 distinguish the pseudorandom sequences from truly random sequences with non-negligible probability,
1497 see, \emph{e.g.},~\cite{Sidorenko:2005:CSB:2179218.2179250}.
1500 This generator is able to pass the whole BigCrush battery of tests, for all
1501 the versions that have been tested depending on their number of threads
1502 (called \texttt{NumThreads} in our algorithm, tested up to $5$ million).
1505 The proposed algorithm has the advantage of manipulating independent
1506 PRNGs, so this version is easily adaptable on a cluster of computers too. The only thing
1507 to ensure is to use a single ISAAC PRNG. To achieve this requirement, a simple solution consists in
1508 using a master node for the initialization. This master node computes the initial parameters
1509 for all the different nodes involved in the computation.
1512 \subsection{Improved Version for GPU}
1514 As GPU cards using CUDA have shared memory between threads of the same block, it
1515 is possible to use this feature in order to simplify the previous algorithm,
1516 i.e., to use less than 3 xor-like PRNGs. The solution consists in computing only
1517 one xor-like PRNG by thread, saving it into the shared memory, and then to use the results
1518 of some other threads in the same block of threads. In order to define which
1519 thread uses the result of which other one, we can use a combination array that
1520 contains the indexes of all threads and for which a combination has been
1523 In Algorithm~\ref{algo:gpu_kernel2}, two combination arrays are used. The
1524 variable \texttt{offset} is computed using the value of
1525 \texttt{combination\_size}. Then we can compute \texttt{o1} and \texttt{o2}
1526 representing the indexes of the other threads whose results are used by the
1527 current one. In this algorithm, we consider that a 32-bits xor-like PRNG has
1528 been chosen. In practice, we use the xor128 proposed in~\cite{Marsaglia2003} in
1529 which unsigned longs (64 bits) have been replaced by unsigned integers (32
1532 This version can also pass the whole {\it BigCrush} battery of tests.
1536 \KwIn{InternalVarXorLikeArray: array with internal variables of 1 xor-like PRNGs
1538 NumThreads: Number of threads\;
1539 array\_comb1, array\_comb2: Arrays containing combinations of size combination\_size\;}
1541 \KwOut{NewNb: array containing random numbers in global memory}
1542 \If{threadId is concerned} {
1543 retrieve data from InternalVarXorLikeArray[threadId] in local variables including shared memory and x\;
1544 offset = threadIdx\%combination\_size\;
1545 o1 = threadIdx-offset+array\_comb1[offset]\;
1546 o2 = threadIdx-offset+array\_comb2[offset]\;
1549 t=t\textasciicircum shmem[o1]\textasciicircum shmem[o2]\;
1550 shared\_mem[threadId]=t\;
1551 x = x\textasciicircum t\;
1553 store the new PRNG in NewNb[NumThreads*threadId+i]\;
1555 store internal variables in InternalVarXorLikeArray[threadId]\;
1558 \caption{Main kernel for the chaotic iterations based PRNG GPU efficient
1560 \label{algo:gpu_kernel2}
1563 \subsection{Chaos Evaluation of the Improved Version}
1565 A run of Algorithm~\ref{algo:gpu_kernel2} consists in an operation ($x=x\oplus t$) having
1566 the form of Equation~\ref{equation Oplus}, which is equivalent to the iterative
1567 system of Eq.~\ref{eq:generalIC}. That is, an iteration of the general chaotic
1568 iterations is realized between the last stored value $x$ of the thread and a strategy $t$
1569 (obtained by a bitwise exclusive or between a value provided by a xor-like() call
1570 and two values previously obtained by two other threads).
1571 To be certain that we are in the framework of Theorem~\ref{t:chaos des general},
1572 we must guarantee that this dynamical system iterates on the space
1573 $\mathcal{X} = \mathcal{P}\left(\llbracket 1, \mathsf{N} \rrbracket\right)^\mathds{N}\times\mathds{B}^\mathsf{N}$.
1574 The left term $x$ obviously belongs to $\mathds{B}^ \mathsf{N}$.
1575 To prevent from any flaws of chaotic properties, we must check that the right
1576 term (the last $t$), corresponding to the strategies, can possibly be equal to any
1577 integer of $\llbracket 1, \mathsf{N} \rrbracket$.
1579 Such a result is obvious, as for the xor-like(), all the
1580 integers belonging into its interval of definition can occur at each iteration, and thus the
1581 last $t$ respects the requirement. Furthermore, it is possible to
1582 prove by an immediate mathematical induction that, as the initial $x$
1583 is uniformly distributed (it is provided by a cryptographically secure PRNG),
1584 the two other stored values shmem[o1] and shmem[o2] are uniformly distributed too,
1585 (this is the induction hypothesis), and thus the next $x$ is finally uniformly distributed.
1587 Thus Algorithm~\ref{algo:gpu_kernel2} is a concrete realization of the general
1588 chaotic iterations presented previously, and for this reason, it satisfies the
1589 Devaney's formulation of a chaotic behavior.
1591 \section{Experiments}
1592 \label{sec:experiments}
1594 Different experiments have been performed in order to measure the generation
1595 speed. We have used a first computer equipped with a Tesla C1060 NVidia GPU card
1597 Intel Xeon E5530 cadenced at 2.40 GHz, and
1598 a second computer equipped with a smaller CPU and a GeForce GTX 280.
1600 cards have 240 cores.
1602 In Figure~\ref{fig:time_xorlike_gpu} we compare the quantity of pseudorandom numbers
1603 generated per second with various xor-like based PRNGs. In this figure, the optimized
1604 versions use the {\it xor64} described in~\cite{Marsaglia2003}, whereas the naive versions
1605 embed the three xor-like PRNGs described in Listing~\ref{algo:seqCIPRNG}. In
1606 order to obtain the optimal performances, the storage of pseudorandom numbers
1607 into the GPU memory has been removed. This step is time consuming and slows down the numbers
1608 generation. Moreover this storage is completely
1609 useless, in case of applications that consume the pseudorandom
1610 numbers directly after generation. We can see that when the number of threads is greater
1611 than approximately 30,000 and lower than 5 million, the number of pseudorandom numbers generated
1612 per second is almost constant. With the naive version, this value ranges from 2.5 to
1613 3GSamples/s. With the optimized version, it is approximately equal to
1614 20GSamples/s. Finally we can remark that both GPU cards are quite similar, but in
1615 practice, the Tesla C1060 has more memory than the GTX 280, and this memory
1616 should be of better quality.
1617 As a comparison, Listing~\ref{algo:seqCIPRNG} leads to the generation of about
1618 138MSample/s when using one core of the Xeon E5530.
1620 \begin{figure}[htbp]
1622 \includegraphics[scale=0.7]{curve_time_xorlike_gpu.pdf}
1624 \caption{Quantity of pseudorandom numbers generated per second with the xorlike-based PRNG}
1625 \label{fig:time_xorlike_gpu}
1632 In Figure~\ref{fig:time_bbs_gpu} we highlight the performances of the optimized
1633 BBS-based PRNG on GPU. On the Tesla C1060 we obtain approximately 700MSample/s
1634 and on the GTX 280 about 670MSample/s, which is obviously slower than the
1635 xorlike-based PRNG on GPU. However, we will show in the next sections that this
1636 new PRNG has a strong level of security, which is necessarily paid by a speed
1639 \begin{figure}[htbp]
1641 \includegraphics[scale=0.7]{curve_time_bbs_gpu.pdf}
1643 \caption{Quantity of pseudorandom numbers generated per second using the BBS-based PRNG}
1644 \label{fig:time_bbs_gpu}
1647 All these experiments allow us to conclude that it is possible to
1648 generate a very large quantity of pseudorandom numbers statistically perfect with the xor-like version.
1649 To a certain extend, it is also the case with the secure BBS-based version, the speed deflation being
1650 explained by the fact that the former version has ``only''
1651 chaotic properties and statistical perfection, whereas the latter is also cryptographically secure,
1652 as it is shown in the next sections.
1660 \section{Security Analysis}
1663 This section is dedicated to the security analysis of the
1664 proposed PRNGs, both from a theoretical and from a practical point of view.
1666 \subsection{Theoretical Proof of Security}
1667 \label{sec:security analysis}
1669 The standard definition
1670 of {\it indistinguishability} used is the classical one as defined for
1671 instance in~\cite[chapter~3]{Goldreich}.
1672 This property shows that predicting the future results of the PRNG
1673 cannot be done in a reasonable time compared to the generation time. It is important to emphasize that this
1674 is a relative notion between breaking time and the sizes of the
1675 keys/seeds. Of course, if small keys or seeds are chosen, the system can
1676 be broken in practice. But it also means that if the keys/seeds are large
1677 enough, the system is secured.
1678 As a complement, an example of a concrete practical evaluation of security
1679 is outlined in the next subsection.
1681 In this section the concatenation of two strings $u$ and $v$ is classically
1683 In a cryptographic context, a pseudorandom generator is a deterministic
1684 algorithm $G$ transforming strings into strings and such that, for any
1685 seed $s$ of length $m$, $G(s)$ (the output of $G$ on the input $s$) has size
1686 $\ell_G(m)$ with $\ell_G(m)>m$.
1687 The notion of {\it secure} PRNGs can now be defined as follows.
1690 A cryptographic PRNG $G$ is secure if for any probabilistic polynomial time
1691 algorithm $D$, for any positive polynomial $p$, and for all sufficiently
1693 $$| \mathrm{Pr}[D(G(U_m))=1]-Pr[D(U_{\ell_G(m)})=1]|< \frac{1}{p(m)},$$
1694 where $U_r$ is the uniform distribution over $\{0,1\}^r$ and the
1695 probabilities are taken over $U_m$, $U_{\ell_G(m)}$ as well as over the
1696 internal coin tosses of $D$.
1699 Intuitively, it means that there is no polynomial time algorithm that can
1700 distinguish a perfect uniform random generator from $G$ with a non negligible
1701 probability. An equivalent formulation of this well-known security property
1702 means that it is possible \emph{in practice} to predict the next bit of the
1703 generator, knowing all the previously produced ones. The interested reader is
1704 referred to~\cite[chapter~3]{Goldreich} for more information. Note that it is
1705 quite easily possible to change the function $\ell$ into any polynomial function
1706 $\ell^\prime$ satisfying $\ell^\prime(m)>m)$~\cite[Chapter 3.3]{Goldreich}.
1708 The generation schema developed in (\ref{equation Oplus}) is based on a
1709 pseudorandom generator. Let $H$ be a cryptographic PRNG. We may assume,
1710 without loss of generality, that for any string $S_0$ of size $N$, the size
1711 of $H(S_0)$ is $kN$, with $k>2$. It means that $\ell_H(N)=kN$.
1712 Let $S_1,\ldots,S_k$ be the
1713 strings of length $N$ such that $H(S_0)=S_1 \ldots S_k$ ($H(S_0)$ is the concatenation of
1714 the $S_i$'s). The cryptographic PRNG $X$ defined in (\ref{equation Oplus})
1715 is the algorithm mapping any string of length $2N$ $x_0S_0$ into the string
1716 $(x_0\oplus S_0 \oplus S_1)(x_0\oplus S_0 \oplus S_1\oplus S_2)\ldots
1717 (x_o\bigoplus_{i=0}^{i=k}S_i)$. One in particular has $\ell_{X}(2N)=kN=\ell_H(N)$.
1718 We claim now that if this PRNG is secure,
1719 then the new one is secure too.
1722 \label{cryptopreuve}
1723 If $H$ is a secure cryptographic PRNG, then $X$ is a secure cryptographic
1728 The proposition is proven by contraposition. Assume that $X$ is not
1729 secure. By Definition, there exists a polynomial time probabilistic
1730 algorithm $D$, a positive polynomial $p$, such that for all $k_0$ there exists
1731 $N\geq \frac{k_0}{2}$ satisfying
1732 $$| \mathrm{Pr}[D(X(U_{2N}))=1]-\mathrm{Pr}[D(U_{kN}=1]|\geq \frac{1}{p(2N)}.$$
1733 We describe a new probabilistic algorithm $D^\prime$ on an input $w$ of size
1736 \item Decompose $w$ into $w=w_1\ldots w_{k}$, where each $w_i$ has size $N$.
1737 \item Pick a string $y$ of size $N$ uniformly at random.
1738 \item Compute $z=(y\oplus w_1)(y\oplus w_1\oplus w_2)\ldots (y
1739 \bigoplus_{i=1}^{i=k} w_i).$
1740 \item Return $D(z)$.
1744 Consider for each $y\in \mathbb{B}^{kN}$ the function $\varphi_{y}$
1745 from $\mathbb{B}^{kN}$ into $\mathbb{B}^{kN}$ mapping $w=w_1\ldots w_k$
1746 (each $w_i$ has length $N$) to
1747 $(y\oplus w_1)(y\oplus w_1\oplus w_2)\ldots (y
1748 \bigoplus_{i=1}^{i=k_1} w_i).$ By construction, one has for every $w$,
1749 \begin{equation}\label{PCH-1}
1750 D^\prime(w)=D(\varphi_y(w)),
1752 where $y$ is randomly generated.
1753 Moreover, for each $y$, $\varphi_{y}$ is injective: if
1754 $(y\oplus w_1)(y\oplus w_1\oplus w_2)\ldots (y\bigoplus_{i=1}^{i=k_1}
1755 w_i)=(y\oplus w_1^\prime)(y\oplus w_1^\prime\oplus w_2^\prime)\ldots
1756 (y\bigoplus_{i=1}^{i=k} w_i^\prime)$, then for every $1\leq j\leq k$,
1757 $y\bigoplus_{i=1}^{i=j} w_i^\prime=y\bigoplus_{i=1}^{i=j} w_i$. It follows,
1758 by a direct induction, that $w_i=w_i^\prime$. Furthermore, since $\mathbb{B}^{kN}$
1759 is finite, each $\varphi_y$ is bijective. Therefore, and using (\ref{PCH-1}),
1761 $\mathrm{Pr}[D^\prime(U_{kN})=1]=\mathrm{Pr}[D(\varphi_y(U_{kN}))=1]$ and,
1763 \begin{equation}\label{PCH-2}
1764 \mathrm{Pr}[D^\prime(U_{kN})=1]=\mathrm{Pr}[D(U_{kN})=1].
1767 Now, using (\ref{PCH-1}) again, one has for every $x$,
1768 \begin{equation}\label{PCH-3}
1769 D^\prime(H(x))=D(\varphi_y(H(x))),
1771 where $y$ is randomly generated. By construction, $\varphi_y(H(x))=X(yx)$,
1773 \begin{equation}%\label{PCH-3} %%RAPH : j'ai viré ce label qui existe déjà, il est 3 ligne avant
1774 D^\prime(H(x))=D(yx),
1776 where $y$ is randomly generated.
1779 \begin{equation}\label{PCH-4}
1780 \mathrm{Pr}[D^\prime(H(U_{N}))=1]=\mathrm{Pr}[D(U_{2N})=1].
1782 From (\ref{PCH-2}) and (\ref{PCH-4}), one can deduce that
1783 there exists a polynomial time probabilistic
1784 algorithm $D^\prime$, a positive polynomial $p$, such that for all $k_0$ there exists
1785 $N\geq \frac{k_0}{2}$ satisfying
1786 $$| \mathrm{Pr}[D(H(U_{N}))=1]-\mathrm{Pr}[D(U_{kN}=1]|\geq \frac{1}{p(2N)},$$
1787 proving that $H$ is not secure, which is a contradiction.
1792 \subsection{Practical Security Evaluation}
1793 \label{sec:Practicak evaluation}
1795 Pseudorandom generators based on Eq.~\eqref{equation Oplus} are thus cryptographically secure when
1796 they are XORed with an already cryptographically
1797 secure PRNG. But, as stated previously,
1798 such a property does not mean that, whatever the
1799 key size, no attacker can predict the next bit
1800 knowing all the previously released ones.
1801 However, given a key size, it is possible to
1802 measure in practice the minimum duration needed
1803 for an attacker to break a cryptographically
1804 secure PRNG, if we know the power of his/her
1805 machines. Such a concrete security evaluation
1806 is related to the $(T,\varepsilon)-$security
1807 notion, which is recalled and evaluated in what
1808 follows, for the sake of completeness.
1810 Let us firstly recall that,
1812 Let $\mathcal{D} : \mathds{B}^M \longrightarrow \mathds{B}$ be a probabilistic algorithm that runs
1814 Let $\varepsilon > 0$.
1815 $\mathcal{D}$ is called a $(T,\varepsilon)-$distinguishing attack on pseudorandom
1818 $$\left| Pr[\mathcal{D}(G(k)) = 1 \mid k \in_R \{0,1\}^\ell ]\right. - \left. Pr[\mathcal{D}(s) = 1 \mid s \in_R \mathds{B}^M ]\right| \geqslant \varepsilon,$$
1819 \noindent where the probability is taken over the internal coin flips of $\mathcal{D}$, and the notation
1820 ``$\in_R$'' indicates the process of selecting an element at random and uniformly over the
1824 Let us recall that the running time of a probabilistic algorithm is defined to be the
1825 maximum of the expected number of steps needed to produce an output, maximized
1826 over all inputs; the expected number is averaged over all coin flips made by the algorithm~\cite{Knuth97}.
1827 We are now able to define the notion of cryptographically secure PRNGs:
1830 A pseudorandom generator is $(T,\varepsilon)-$secure if there exists no $(T,\varepsilon)-$distinguishing attack on this pseudorandom generator.
1839 Suppose now that the PRNG of Eq.~\eqref{equation Oplus} will work during
1840 $M=100$ time units, and that during this period,
1841 an attacker can realize $10^{12}$ clock cycles.
1842 We thus wonder whether, during the PRNG's
1843 lifetime, the attacker can distinguish this
1844 sequence from a truly random one, with a probability
1845 greater than $\varepsilon = 0.2$.
1846 We consider that $N$ has 900 bits.
1848 Predicting the next generated bit knowing all the
1849 previously released ones by Eq.~\eqref{equation Oplus} is obviously equivalent to predicting the
1850 next bit in the BBS generator, which
1851 is cryptographically secure. More precisely, it
1852 is $(T,\varepsilon)-$secure: no
1853 $(T,\varepsilon)-$distinguishing attack can be
1854 successfully realized on this PRNG, if~\cite{Fischlin}
1856 T \leqslant \dfrac{L(N)}{6 N (log_2(N))\varepsilon^{-2}M^2}-2^7 N \varepsilon^{-2} M^2 log_2 (8 N \varepsilon^{-1}M)
1857 \label{mesureConcrete}
1859 where $M$ is the length of the output ($M=100$ in
1860 our example), and $L(N)$ is equal to
1862 2.8\times 10^{-3} exp \left(1.9229 \times (N ~ln~ 2)^\frac{1}{3} \times (ln(N~ln~ 2))^\frac{2}{3}\right)
1864 is the number of clock cycles to factor a $N-$bit
1870 A direct numerical application shows that this attacker
1871 cannot achieve its $(10^{12},0.2)$ distinguishing
1872 attack in that context.
1876 \section{Cryptographical Applications}
1878 \subsection{A Cryptographically Secure PRNG for GPU}
1881 It is possible to build a cryptographically secure PRNG based on the previous
1882 algorithm (Algorithm~\ref{algo:gpu_kernel2}). Due to Proposition~\ref{cryptopreuve},
1883 it simply consists in replacing
1884 the {\it xor-like} PRNG by a cryptographically secure one.
1885 We have chosen the Blum Blum Shub generator~\cite{BBS} (usually denoted by BBS) having the form:
1886 $$x_{n+1}=x_n^2~ mod~ M$$ where $M$ is the product of two prime numbers (these
1887 prime numbers need to be congruent to 3 modulus 4). BBS is known to be
1888 very slow and only usable for cryptographic applications.
1891 The modulus operation is the most time consuming operation for current
1892 GPU cards. So in order to obtain quite reasonable performances, it is
1893 required to use only modulus on 32-bits integer numbers. Consequently
1894 $x_n^2$ need to be lesser than $2^{32}$, and thus the number $M$ must be
1895 lesser than $2^{16}$. So in practice we can choose prime numbers around
1896 256 that are congruent to 3 modulus 4. With 32-bits numbers, only the
1897 4 least significant bits of $x_n$ can be chosen (the maximum number of
1898 indistinguishable bits is lesser than or equals to
1899 $log_2(log_2(M))$). In other words, to generate a 32-bits number, we need to use
1900 8 times the BBS algorithm with possibly different combinations of $M$. This
1901 approach is not sufficient to be able to pass all the tests of TestU01,
1902 as small values of $M$ for the BBS lead to
1903 small periods. So, in order to add randomness we have proceeded with
1904 the followings modifications.
1907 Firstly, we define 16 arrangement arrays instead of 2 (as described in
1908 Algorithm \ref{algo:gpu_kernel2}), but only 2 of them are used at each call of
1909 the PRNG kernels. In practice, the selection of combination
1910 arrays to be used is different for all the threads. It is determined
1911 by using the three last bits of two internal variables used by BBS.
1912 %This approach adds more randomness.
1913 In Algorithm~\ref{algo:bbs_gpu},
1914 character \& is for the bitwise AND. Thus using \&7 with a number
1915 gives the last 3 bits, thus providing a number between 0 and 7.
1917 Secondly, after the generation of the 8 BBS numbers for each thread, we
1918 have a 32-bits number whose period is possibly quite small. So
1919 to add randomness, we generate 4 more BBS numbers to
1920 shift the 32-bits numbers, and add up to 6 new bits. This improvement is
1921 described in Algorithm~\ref{algo:bbs_gpu}. In practice, the last 2 bits
1922 of the first new BBS number are used to make a left shift of at most
1923 3 bits. The last 3 bits of the second new BBS number are added to the
1924 strategy whatever the value of the first left shift. The third and the
1925 fourth new BBS numbers are used similarly to apply a new left shift
1928 Finally, as we use 8 BBS numbers for each thread, the storage of these
1929 numbers at the end of the kernel is performed using a rotation. So,
1930 internal variable for BBS number 1 is stored in place 2, internal
1931 variable for BBS number 2 is stored in place 3, ..., and finally, internal
1932 variable for BBS number 8 is stored in place 1.
1937 \KwIn{InternalVarBBSArray: array with internal variables of the 8 BBS
1939 NumThreads: Number of threads\;
1940 array\_comb: 2D Arrays containing 16 combinations (in first dimension) of size combination\_size (in second dimension)\;
1941 array\_shift[4]=\{0,1,3,7\}\;
1944 \KwOut{NewNb: array containing random numbers in global memory}
1945 \If{threadId is concerned} {
1946 retrieve data from InternalVarBBSArray[threadId] in local variables including shared memory and x\;
1947 we consider that bbs1 ... bbs8 represent the internal states of the 8 BBS numbers\;
1948 offset = threadIdx\%combination\_size\;
1949 o1 = threadIdx-offset+array\_comb[bbs1\&7][offset]\;
1950 o2 = threadIdx-offset+array\_comb[8+bbs2\&7][offset]\;
1957 \tcp{two new shifts}
1958 shift=BBS3(bbs3)\&3\;
1960 t|=BBS1(bbs1)\&array\_shift[shift]\;
1961 shift=BBS7(bbs7)\&3\;
1963 t|=BBS2(bbs2)\&array\_shift[shift]\;
1964 t=t\textasciicircum shmem[o1]\textasciicircum shmem[o2]\;
1965 shared\_mem[threadId]=t\;
1966 x = x\textasciicircum t\;
1968 store the new PRNG in NewNb[NumThreads*threadId+i]\;
1970 store internal variables in InternalVarXorLikeArray[threadId] using a rotation\;
1973 \caption{main kernel for the BBS based PRNG GPU}
1974 \label{algo:bbs_gpu}
1977 In Algorithm~\ref{algo:bbs_gpu}, $n$ is for the quantity of random numbers that
1978 a thread has to generate. The operation t<<=4 performs a left shift of 4 bits
1979 on the variable $t$ and stores the result in $t$, and $BBS1(bbs1)\&15$ selects
1980 the last four bits of the result of $BBS1$. Thus an operation of the form
1981 $t<<=4; t|=BBS1(bbs1)\&15\;$ realizes in $t$ a left shift of 4 bits, and then
1982 puts the 4 last bits of $BBS1(bbs1)$ in the four last positions of $t$. Let us
1983 remark that the initialization $t$ is not a necessity as we fill it 4 bits by 4
1984 bits, until having obtained 32-bits. The two last new shifts are realized in
1985 order to enlarge the small periods of the BBS used here, to introduce a kind of
1986 variability. In these operations, we make twice a left shift of $t$ of \emph{at
1987 most} 3 bits, represented by \texttt{shift} in the algorithm, and we put
1988 \emph{exactly} the \texttt{shift} last bits from a BBS into the \texttt{shift}
1989 last bits of $t$. For this, an array named \texttt{array\_shift}, containing the
1990 correspondence between the shift and the number obtained with \texttt{shift} 1
1991 to make the \texttt{and} operation is used. For example, with a left shift of 0,
1992 we make an and operation with 0, with a left shift of 3, we make an and
1993 operation with 7 (represented by 111 in binary mode).
1995 It should be noticed that this generator has once more the form $x^{n+1} = x^n \oplus S^n$,
1996 where $S^n$ is referred in this algorithm as $t$: each iteration of this
1997 PRNG ends with $x = x \wedge t$. This $S^n$ is only constituted
1998 by secure bits produced by the BBS generator, and thus, due to
1999 Proposition~\ref{cryptopreuve}, the resulted PRNG is
2000 cryptographically secure.
2002 As stated before, even if the proposed PRNG is cryptocaphically
2003 secure, it does not mean that such a generator
2004 can be used as described here when attacks are
2005 awaited. The problem is to determine the minimum
2006 time required for an attacker, with a given
2007 computational power, to predict under a probability
2008 lower than 0.5 the $n+1$th bit, knowing the $n$
2009 previous ones. The proposed GPU generator will be
2010 useful in a security context, at least in some
2011 situations where a secret protected by a pseudorandom
2012 keystream is rapidly obsolete, if this time to
2013 predict the next bit is large enough when compared
2014 to both the generation and transmission times.
2015 It is true that the prime numbers used in the last
2016 section are very small compared to up-to-date
2017 security recommendations. However the attacker has not
2018 access to each BBS, but to the output produced
2019 by Algorithm~\ref{algo:bbs_gpu}, which is far
2020 more complicated than a simple BBS. Indeed, to
2021 determine if this cryptographically secure PRNG
2022 on GPU can be useful in security context with the
2023 proposed parameters, or if it is only a very fast
2024 and statistically perfect generator on GPU, its
2025 $(T,\varepsilon)-$security must be determined, and
2026 a formulation similar to Eq.\eqref{mesureConcrete}
2027 must be established. Authors
2028 hope to achieve this difficult task in a future
2032 \subsection{Toward a Cryptographically Secure and Chaotic Asymmetric Cryptosystem}
2033 \label{Blum-Goldwasser}
2034 We finish this research work by giving some thoughts about the use of
2035 the proposed PRNG in an asymmetric cryptosystem.
2036 This first approach will be further investigated in a future work.
2038 \subsubsection{Recalls of the Blum-Goldwasser Probabilistic Cryptosystem}
2040 The Blum-Goldwasser cryptosystem is a cryptographically secure asymmetric key encryption algorithm
2041 proposed in 1984~\cite{Blum:1985:EPP:19478.19501}. The encryption algorithm
2042 implements a XOR-based stream cipher using the BBS PRNG, in order to generate
2043 the keystream. Decryption is done by obtaining the initial seed thanks to
2044 the final state of the BBS generator and the secret key, thus leading to the
2045 reconstruction of the keystream.
2047 The key generation consists in generating two prime numbers $(p,q)$,
2048 randomly and independently of each other, that are
2049 congruent to 3 mod 4, and to compute the modulus $N=pq$.
2050 The public key is $N$, whereas the secret key is the factorization $(p,q)$.
2053 Suppose Bob wishes to send a string $m=(m_0, \dots, m_{L-1})$ of $L$ bits to Alice:
2055 \item Bob picks an integer $r$ randomly in the interval $\llbracket 1,N\rrbracket$ and computes $x_0 = r^2~mod~N$.
2056 \item He uses the BBS to generate the keystream of $L$ pseudorandom bits $(b_0, \dots, b_{L-1})$, as follows. For $i=0$ to $L-1$,
2059 \item While $i \leqslant L-1$:
2061 \item Set $b_i$ equal to the least-significant\footnote{As signaled previously, BBS can securely output up to $\mathsf{N} = \lfloor log(log(N)) \rfloor$ of the least-significant bits of $x_i$ during each round.} bit of $x_i$,
2063 \item $x_i = (x_{i-1})^2~mod~N.$
2066 \item The ciphertext is computed by XORing the plaintext bits $m$ with the keystream: $ c = (c_0, \dots, c_{L-1}) = m \oplus b$. This ciphertext is $[c, y]$, where $y=x_{0}^{2^{L}}~mod~N.$
2070 When Alice receives $\left[(c_0, \dots, c_{L-1}), y\right]$, she can recover $m$ as follows:
2072 \item Using the secret key $(p,q)$, she computes $r_p = y^{((p+1)/4)^{L}}~mod~p$ and $r_q = y^{((q+1)/4)^{L}}~mod~q$.
2073 \item The initial seed can be obtained using the following procedure: $x_0=q(q^{-1}~{mod}~p)r_p + p(p^{-1}~{mod}~q)r_q~{mod}~N$.
2074 \item She recomputes the bit-vector $b$ by using BBS and $x_0$.
2075 \item Alice finally computes the plaintext by XORing the keystream with the ciphertext: $ m = c \oplus b$.
2079 \subsubsection{Proposal of a new Asymmetric Cryptosystem Adapted from Blum-Goldwasser}
2081 We propose to adapt the Blum-Goldwasser protocol as follows.
2082 Let $\mathsf{N} = \lfloor log(log(N)) \rfloor$ be the number of bits that can
2083 be obtained securely with the BBS generator using the public key $N$ of Alice.
2084 Alice will pick randomly $S^0$ in $\llbracket 0, 2^{\mathsf{N}-1}\rrbracket$ too, and
2085 her new public key will be $(S^0, N)$.
2087 To encrypt his message, Bob will compute
2088 %%RAPH : ici, j'ai mis un simple $
2090 c = \left(m_0 \oplus (b_0 \oplus S^0), m_1 \oplus (b_0 \oplus b_1 \oplus S^0), \hdots, \right.
2091 \left. m_{L-1} \oplus (b_0 \oplus b_1 \hdots \oplus b_{L-1} \oplus S^0) \right)
2093 instead of $$\left(m_0 \oplus b_0, m_1 \oplus b_1, \hdots, m_{L-1} \oplus b_{L-1} \right).$$
2095 The same decryption stage as in Blum-Goldwasser leads to the sequence
2096 $$\left(m_0 \oplus S^0, m_1 \oplus S^0, \hdots, m_{L-1} \oplus S^0 \right).$$
2097 Thus, with a simple use of $S^0$, Alice can obtain the plaintext.
2098 By doing so, the proposed generator is used in place of BBS, leading to
2099 the inheritance of all the properties presented in this paper.
2101 \section{Conclusion}
2104 In this paper, a formerly proposed PRNG based on chaotic iterations
2105 has been generalized to improve its speed. It has been proven to be
2106 chaotic according to Devaney.
2107 Efficient implementations on GPU using xor-like PRNGs as input generators
2108 have shown that a very large quantity of pseudorandom numbers can be generated per second (about
2109 20Gsamples/s), and that these proposed PRNGs succeed to pass the hardest battery in TestU01,
2110 namely the BigCrush.
2111 Furthermore, we have shown that when the inputted generator is cryptographically
2112 secure, then it is the case too for the PRNG we propose, thus leading to
2113 the possibility to develop fast and secure PRNGs using the GPU architecture.
2114 An improvement of the Blum-Goldwasser cryptosystem, making it
2115 behave chaotically, has finally been proposed.
2117 In future work we plan to extend this research, building a parallel PRNG for clusters or
2118 grid computing. Topological properties of the various proposed generators will be investigated,
2119 and the use of other categories of PRNGs as input will be studied too. The improvement
2120 of Blum-Goldwasser will be deepened.
2122 Another aspect to consider might be different accelerator-based systems like
2123 Intel Xeon Phi cards and speed measurements using such cards: as heterogeneity of
2124 supercomputers tends to increase using other accelerators than GPGPUs,
2125 a Xeon Phi solution might be interesting to investigate.
2128 will try to enlarge the quantity of pseudorandom numbers generated per second either
2129 in a simulation context or in a cryptographic one.
2133 \bibliographystyle{plain}
2134 \bibliography{mabase}