1 \chapterauthor{Raphaël Couturier}{Femto-ST Institute, University of Franche-Comt\'{e}}
2 \chapterauthor{Christophe Guyeux}{Femto-ST Institute, University of Franche-Comt\'{e}}
5 \chapter{Pseudorandom Number Generator on GPU}
10 Randomness is of importance in many fields such as scientific
11 simulations or cryptography. ``Random numbers'' can mainly be
12 generated either by a deterministic and reproducible algorithm called
13 a pseudorandom number generator (PRNG), or by a physical non-deterministic
14 process having all the characteristics of a random noise, called a truly random number
15 generator (TRNG). In this chapter, we focus on
16 reproducible generators, useful for instance in Monte-Carlo based
17 simulators. These domains need PRNGs that are statistically
18 irreproachable. In some fields such as in numerical simulations,
19 speed is a strong requirement that is usually attained by using
20 parallel architectures. In that case, a recurrent problem is that a
21 deflation of the statistical qualities is often reported, when the
22 parallelization of a good PRNG is realized. In some fields such as in
23 numerical simulations, speed is a strong requirement that is usually
24 attained by using parallel architectures. In that case, a recurrent
25 problem is that a deflation of the statistical qualities is often
26 reported, when the parallelization of a good PRNG is realized. This
27 is why ad-hoc PRNGs for each possible architecture must be found to
28 achieve both speed and randomness. On the other side, speed is not
29 the main requirement in cryptography: the great need is to
30 define \emph{secure} generators able to withstand malicious
31 attacks. Roughly speaking, an attacker should not be able in practice
32 to make the distinction between numbers obtained with the secure
33 generator and a true random sequence. Or, in an equivalent
34 formulation, he or she should not be able (in practice) to predict the
35 next bit of the generator, having the knowledge of all the binary
36 digits that have been already released~\cite{Goldreich}. ``Being able in practice''
37 refers here to the possibility to achieve this attack in polynomial
38 time, and to the exponential growth of the difficulty of this
39 challenge when the size of the parameters of the PRNG increases.
42 Finally, a small part of the community working in this domain focuses on a
43 third requirement, that is to define chaotic generators~\cite{kellert1994wake, Wu20051195,gleick2011chaos}.
44 The main idea is to take benefits from a chaotic dynamical system to obtain a
45 generator that is unpredictable, disordered, sensible to its seed, or in other word chaotic.
46 Their desire is to map a given chaotic dynamics into a sequence that seems random
47 and unassailable due to chaos.
48 However, the chaotic maps used as a pattern are defined in the real line
49 whereas computers deal with finite precision numbers.
50 This distortion leads to a deflation of both chaotic properties and speed.
51 Furthermore, authors of such chaotic generators often claim their PRNG
52 as secure due to their chaos properties, but there is no obvious relation
53 between chaos and security as it is understood in cryptography.
54 This is why the use of chaos for PRNG still remains marginal and disputable.
55 However, we have established in previous researches that these flaws can
56 be circumvented by using a tool called choatic iterations.
57 Such investigations have led to the definition of a new
58 family of PRNGs that are chaotic while being fast and statistically perfect,
59 or cryptographically secure~\cite{bgw09:ip,bgw10:ip,bfgw11:ij,bfg12a:ip}. This family is improved and adapted for GPU
60 architectures in this chapter.
63 Let us finish this introduction by noticing that, in this paper,
64 statistical perfection refers to the ability to pass the whole
65 {\it BigCrush} battery of tests, which is widely considered as the most
66 stringent statistical evaluation of a sequence claimed as random.
67 This battery can be found in the well-known TestU01 package~\cite{LEcuyerS07}.
68 More precisely, each time we performed a test on a PRNG, we ran it
69 twice in order to observe if all $p-$values are inside [0.01, 0.99]. In
70 fact, we observed that few $p-$values (less than ten) are sometimes
71 outside this interval but inside [0.001, 0.999], so that is why a
72 second run allows us to confirm that the values outside are not for
73 the same test. With this approach all our PRNGs pass the {\it
74 BigCrush} successfully and all $p-$values are at least once inside
76 Chaos, for its part, refers to the well-established definition of a
77 chaotic dynamical system defined by Devaney~\cite{Devaney}.
79 The remainder of this chapter is organized as follows.
80 Basic definitions and terminologies about both topological chaos
81 and chaotic iterations are provided in the next section. Various
82 chaotic iterations based pseudorandom number generators are then
83 presented in Section~\ref{sec:efficient PRNG}. They encompass
84 naive and improved efficient generators for CPU and for GPU.
85 These generators are finally experimented in Section~\ref{sec:experiments}.
88 \section{Basic Recalls}
89 \label{section:BASIC RECALLS}
91 This section is devoted to basic definitions and terminologies in the fields of
92 topological chaos and chaotic iterations. We assume the reader is familiar
93 with basic notions on topology (see for instance~\cite{Devaney}).
96 \subsection{A Short Presentation of Chaos}
99 Chaos theory studies the behavior of dynamical systems that are perfectly predictable, yet appear to be wildly amorphous and without meaningful.
100 Chaotic systems are highly sensitive to initial conditions,
101 which is popularly referred to as the butterfly effect.
102 In other words, small differences in initial conditions (such as those due to rounding errors in numerical computation) yield widely diverging outcomes,
103 rendering long-term prediction impossible in general \cite{kellert1994wake}. This happens even though these systems are deterministic, meaning that their future behavior is fully determined by their initial conditions, with no random elements involved \cite{kellert1994wake}. That is, the deterministic nature of these systems does not make them predictable \cite{kellert1994wake,Werndl01032009}. This behavior is known as deterministic chaos, or simply chaos. It has been well-studied in mathematics and
104 physics, leading among other things to the well-established definition of Devaney
111 \subsection{On Devaney's Definition of Chaos}
113 Consider a metric space $(\mathcal{X},d)$ and a continuous function $f:\mathcal{X}\longrightarrow \mathcal{X}$, for one-dimensional dynamical systems of the form:
115 x^0 \in \mathcal{X} \textrm{ and } \forall n \in \mathds{N}^*, x^n=f(x^{n-1}),
118 the following definition of chaotic behavior, formulated by Devaney~\cite{Devaney}, is widely accepted.
121 A dynamical system of Form~(\ref{Devaney}) is said to be chaotic if the following conditions hold.
123 \item Topological transitivity:
126 \forall U,V \textrm{ open sets of } \mathcal{X},~\exists k>0, f^k(U) \cap V \neq \varnothing .
129 Intuitively, a topologically transitive map has points that eventually move under iteration from one arbitrarily small neighborhood to any other. Consequently, the dynamical system cannot be decomposed into two disjoint open sets that are invariant under the map. Note that if a map possesses a dense orbit, then it is clearly topologically transitive.
130 \item Density of periodic points in $\mathcal{X}$.
132 Let $P=\{p\in \mathcal{X}|\exists n \in \mathds{N}^{\ast}:f^n(p)=p\}$ the set of periodic points of $f$. Then $P$ is dense in $\mathcal{X}$:
135 \overline{P}=\mathcal{X} .
138 Density of periodic orbits means that every point in the space is approached arbitrarily closely by periodic orbits. Topologically mixing systems failing this condition may not display sensitivity to initial conditions presented below, and hence may not be chaotic.
139 \item Sensitive dependence on initial conditions:
141 $\exists \varepsilon>0,$ $\forall x \in \mathcal{X},$ $\forall \delta >0,$ $\exists y \in \mathcal{X},$ $\exists n \in \mathbb{N},$ $d(x,y)<\delta$ and $d\left(f^n(x),f^n(y)\right) \geqslant \varepsilon.$
143 Intuitively, a map possesses sensitive dependence on initial conditions if there exist points arbitrarily close to $x$ that eventually separate from $x$ by at least $\varepsilon$ under iteration of $f$. Not all points near $x$ need eventually separate from $x$ under iteration, but there must be at least one such point in every neighborhood of $x$. If a map possesses sensitive dependence on initial conditions, then for all practical purposes, the dynamics of the map defy numerical computation. Small errors in computation that are introduced by round-off may become magnified upon iteration. The results of numerical computation of an orbit, no matter how accurate, may bear no resemblance whatsoever with the real orbit.
148 When $f$ is chaotic, then the system $(\mathcal{X}, f)$ is chaotic and quoting Devaney: ``it is unpredictable because of the sensitive dependence on initial conditions. It cannot be broken down or decomposed into two subsystems which do not interact because of topological transitivity. And, in the midst of this random behavior, we nevertheless have an element of regularity.'' Fundamentally different behaviors are consequently possible and occur in an unpredictable way.
157 \subsection{Chaotic iterations}
158 \label{subsection:Chaotic iterations}
160 Let us now introduce an example of a dynamical systems family that has
161 the potentiality to become chaotic, depending on the choice of the iteration
162 function. This family is the basis of the PRNGs we will
163 develop during this chapter.
166 \label{Chaotic iterations}
167 The set $\mathds{B}$ denoting $\{0,1\}$, let $f:\mathds{B}^{\mathsf{N}%
168 }\longrightarrow \mathds{B}^{\mathsf{N}}$ be an ``iteration'' function and $S$ be a
169 sequence of subsets of $\llbracket 1, \mathsf{N}\rrbracket$ called a chaotic strategy. Then, the so-called \emph{chaotic iterations} are defined by~\cite{Robert1986}:
173 \left\{\begin{array}{l}
174 x^0\in \mathds{B}^{\mathsf{N}}, \\
175 \forall n\in \mathds{N}^{\ast },\forall i\in \llbracket1;\mathsf{N}\rrbracket%
179 x_i^{n-1} & \text{if}~ i\notin S^n \\
180 f(x^{n-1})_{i} & \text{if}~ i \in S^n.
188 In other words, at the $n^{th}$ iteration, only the cells of $S^{n}$ are
190 Chaotic iterations generate a set of vectors;
191 they are defined by an initial state $x^{0}$, an iteration function $f$, and a chaotic strategy $S$~\cite{bg10:ij}.
192 These ``chaotic iterations'' can behave chaotically as defined by Devaney,
193 depending on the choice of $f$~\cite{bg10:ij}. For instance,
194 chaos is obtained when $f$ is the vectorial negation.
195 Note that, with this example of function, chaotic iterations
196 defined above can be rewritten as:
198 \label{equation Oplus}
199 x^0 \in \llbracket 0,2^\mathsf{N}-1\rrbracket,~\mathcal{S}^n \in \mathcal{P}\left(\llbracket 1,2^\mathsf{N}-1\rrbracket\right)^\mathds{N},~~ x^{n+1}=x^n \oplus \mathcal{S}^n,
201 where $\mathcal{P}(X)$ stands for the set of subsets of $X$, whereas
202 $a\oplus b$ is the bitwise exclusive or operation between the binary
203 representation of the integers $a$ and $b$. Note that the term
204 $\mathcal{S}^n$ is directly and obviously linked to the $S^n$ of
205 Eq.\ref{eq:generalIC}. As recalled above, such an iterative sequence
206 satisfies the Devaney's definition of chaos.
208 \section{Toward Efficiency and Improvement for CI PRNG}
209 \label{sec:efficient PRNG}
211 \subsection{First Efficient Implementation of a PRNG based on Chaotic Iterations}
213 %Based on the proof presented in the previous section, it is now possible to
214 %improve the speed of the generator formerly presented in~\cite{bgw09:ip,guyeux10}.
215 %The first idea is to consider
216 %that the provided strategy is a pseudorandom Boolean vector obtained by a
218 %An iteration of the system is simply the bitwise exclusive or between
219 %the last computed state and the current strategy.
220 %Topological properties of disorder exhibited by chaotic
221 %iterations can be inherited by the inputted generator, we hope by doing so to
222 %obtain some statistical improvements while preserving speed.
224 %%RAPH : j'ai viré tout ca
225 %% Let us give an example using 16-bits numbers, to clearly understand how the bitwise xor operations
228 %% Suppose that $x$ and the strategy $S^i$ are given as
230 %% Table~\ref{TableExemple} shows the result of $x \oplus S^i$.
233 %% \begin{scriptsize}
235 %% \begin{array}{|cc|cccccccccccccccc|}
237 %% x &=&1&0&1&1&1&0&1&0&1&0&0&1&0&0&1&0\\
239 %% S^i &=&0&1&1&0&0&1&1&0&1&1&1&0&0&1&1&1\\
241 %% x \oplus S^i&=&1&1&0&1&1&1&0&0&0&1&1&1&0&1&0&1\\
248 %% \caption{Example of an arbitrary round of the proposed generator}
249 %% \label{TableExemple}
255 \lstset{language=C,caption={C code of the sequential PRNG based on chaotic iterations},label={algo:seqCIPRNG}}
259 unsigned int CIPRNG() {
260 static unsigned int x = 123123123;
261 unsigned long t1 = xorshift();
262 unsigned long t2 = xor128();
263 unsigned long t3 = xorwow();
264 x = x^(unsigned int)t1;
265 x = x^(unsigned int)(t2>>32);
266 x = x^(unsigned int)(t3>>32);
267 x = x^(unsigned int)t2;
268 x = x^(unsigned int)(t1>>32);
269 x = x^(unsigned int)t3;
277 In Listing~\ref{algo:seqCIPRNG} a sequential version of the proposed PRNG based
278 on chaotic iterations is presented, which extends the generator family
279 formerly presented in~\cite{bgw09:ip,guyeux10}. The xor operator is represented by
280 \textasciicircum. This function uses three classical 64-bits PRNGs, namely the
281 \texttt{xorshift}, the \texttt{xor128}, and the
282 \texttt{xorwow}~\cite{Marsaglia2003}. In the following, we call them ``xor-like
283 PRNGs''. As each xor-like PRNG uses 64-bits whereas our proposed generator
284 works with 32-bits, we use the command \texttt{(unsigned int)}, that selects the
285 32 least significant bits of a given integer, and the code \texttt{(unsigned
286 int)(t$>>$32)} in order to obtain the 32 most significant bits of \texttt{t}.
288 Thus producing a pseudorandom number needs 6 xor operations with 6 32-bits numbers
289 that are provided by 3 64-bits PRNGs. This version successfully passes the
290 stringent BigCrush battery of tests~\cite{LEcuyerS07}.
291 At this point, we thus
292 have defined an efficient and statistically unbiased generator. Its speed is
293 directly related to the use of linear operations, but for the same reason,
294 this fast generator cannot be proven as secure.
298 \subsection{Efficient PRNGs based on Chaotic Iterations on GPU}
299 \label{sec:efficient PRNG gpu}
301 In order to take benefits from the computing power of GPU, a program
302 needs to have independent blocks of threads that can be computed
303 simultaneously. In general, the larger the number of threads is, the
304 more local memory is used, and the less branching instructions are
305 used (if, while, ...), the better the performances on GPU is.
306 Obviously, having these requirements in mind, it is possible to build
307 a program similar to the one presented in Listing
308 \ref{algo:seqCIPRNG}, which computes pseudorandom numbers on GPU. To
309 do so, we must firstly recall that in the CUDA~\cite{Nvid10}
310 environment, threads have a local identifier called
311 \texttt{ThreadIdx}, which is relative to the block containing
312 them. Furthermore, in CUDA, parts of the code that are executed by the GPU, are
313 called {\it kernels}.
316 \subsection{Naive Version for GPU}
319 It is possible to deduce from the CPU version a quite similar version adapted to GPU.
320 The simple principle consists in making each thread of the GPU computing the CPU version of our PRNG.
321 Of course, the three xor-like
322 PRNGs used in these computations must have different parameters.
323 In a given thread, these parameters are
324 randomly picked from another PRNGs.
325 The initialization stage is performed by the CPU.
326 To do it, the ISAAC PRNG~\cite{Jenkins96} is used to set all the
327 parameters embedded into each thread.
329 The implementation of the three
330 xor-like PRNGs is straightforward when their parameters have been
331 allocated in the GPU memory. Each xor-like works with an internal
332 number $x$ that saves the last generated pseudorandom number. Additionally, the
333 implementation of the xor128, the xorshift, and the xorwow respectively require
334 4, 5, and 6 unsigned long as internal variables.
339 \KwIn{InternalVarXorLikeArray: array with internal variables of the 3 xor-like
340 PRNGs in global memory\;
341 NumThreads: number of threads\;}
342 \KwOut{NewNb: array containing random numbers in global memory}
343 \If{threadIdx is concerned by the computation} {
344 retrieve data from InternalVarXorLikeArray[threadIdx] in local variables\;
346 compute a new PRNG as in Listing\ref{algo:seqCIPRNG}\;
347 store the new PRNG in NewNb[NumThreads*threadIdx+i]\;
349 store internal variables in InternalVarXorLikeArray[threadIdx]\;
352 \caption{Main kernel of the GPU ``naive'' version of the PRNG based on chaotic iterations}
353 \label{algo:gpu_kernel}
358 Algorithm~\ref{algo:gpu_kernel} presents a naive implementation of the proposed PRNG on
359 GPU. Due to the available memory in the GPU and the number of threads
360 used simultaneously, the number of random numbers that a thread can generate
361 inside a kernel is limited (\emph{i.e.}, the variable \texttt{n} in
362 algorithm~\ref{algo:gpu_kernel}). For instance, if $100,000$ threads are used and
363 if $n=100$\footnote{in fact, we need to add the initial seed (a 32-bits number)},
364 then the memory required to store all of the internals variables of both the xor-like
365 PRNGs\footnote{we multiply this number by $2$ in order to count 32-bits numbers}
366 and the pseudorandom numbers generated by our PRNG, is equal to $100,000\times ((4+5+6)\times
367 2+(1+100))=1,310,000$ 32-bits numbers, that is, approximately $52$Mb.
369 This generator is able to pass the whole BigCrush battery of tests, for all
370 the versions that have been tested depending on their number of threads
371 (called \texttt{NumThreads} in our algorithm, tested up to $5$ million).
374 \subsection{Improved Version for GPU}
376 As GPU cards using CUDA have shared memory between threads of the same block, it
377 is possible to use this feature in order to simplify the previous algorithm,
378 i.e., to use less than 3 xor-like PRNGs. The solution consists in computing only
379 one xor-like PRNG by thread, saving it into the shared memory, and then to use the results
380 of some other threads in the same block of threads. In order to define which
381 thread uses the result of which other one, we can use a combination array that
382 contains the indexes of all threads and for which a combination has been
385 In Algorithm~\ref{algo:gpu_kernel2}, two combination arrays are used. The
386 variable \texttt{offset} is computed using the value of
387 \texttt{combination\_size}. Then we can compute \texttt{o1} and \texttt{o2}
388 representing the indexes of the other threads whose results are used by the
389 current one. In this algorithm, we consider that a 32-bits xor-like PRNG has
390 been chosen. In practice, we use the xor128 proposed in~\cite{Marsaglia2003} in
391 which unsigned longs (64 bits) have been replaced by unsigned integers (32
394 This version can also pass the whole {\it BigCrush} battery of tests.
398 \KwIn{InternalVarXorLikeArray: array with internal variables of 1 xor-like PRNGs
400 NumThreads: Number of threads\;
401 array\_comb1, array\_comb2: Arrays containing combinations of size combination\_size\;}
403 \KwOut{NewNb: array containing random numbers in global memory}
404 \If{threadId is concerned} {
405 retrieve data from InternalVarXorLikeArray[threadId] in local variables including shared memory and x\;
406 offset = threadIdx\%combination\_size\;
407 o1 = threadIdx-offset+array\_comb1[offset]\;
408 o2 = threadIdx-offset+array\_comb2[offset]\;
411 t=t\textasciicircum shmem[o1]\textasciicircum shmem[o2]\;
412 shared\_mem[threadId]=t\;
413 x = x\textasciicircum t\;
415 store the new PRNG in NewNb[NumThreads*threadId+i]\;
417 store internal variables in InternalVarXorLikeArray[threadId]\;
420 \caption{Main kernel for the chaotic iterations based PRNG GPU efficient
422 \label{algo:gpu_kernel2}
425 \subsection{Chaos Evaluation of the Improved Version}
427 A run of Algorithm~\ref{algo:gpu_kernel2} consists in an operation ($x=x\oplus t$) having
428 the form of Equation~\ref{equation Oplus}, which is equivalent to the iterative
429 system of Eq.~\ref{eq:generalIC}. That is, an iteration of the general chaotic
430 iterations is realized between the last stored value $x$ of the thread and a strategy $t$
431 (obtained by a bitwise exclusive or between a value provided by a xor-like() call
432 and two values previously obtained by two other threads).
433 To be certain that such iterations corresponds to the chaotic one recalled at the
434 end of the Section~\ref{sec:dev},
435 we must guarantee that this dynamical system iterates on the space
436 $\mathcal{X} =\mathds{B}^\mathsf{N} \times \mathcal{P}\left(\llbracket 1, 2^\mathsf{N} \rrbracket\right)^\mathds{N}$.
437 The left term $x$ obviously belongs to $\mathds{B}^ \mathsf{N}$.
438 To prevent from any flaws of chaotic properties, we must check that the right
439 term (the last $t$), corresponding to the strategies, can possibly be equal to any
440 integer of $\llbracket 1, 2^\mathsf{N} \rrbracket$.
442 Such a result is obvious, as for the xor-like(), all the
443 integers belonging into its interval of definition can occur at each iteration, and thus the
444 last $t$ respects the requirement. Furthermore, it is possible to
445 prove by an immediate mathematical induction that, supposes that the initial $x$
446 is uniformly distributed, %(it is provided by a cryptographically secure PRNG),
447 the two other stored values shmem[o1] and shmem[o2] are uniformly distributed too,
448 (this is the induction hypothesis), and thus the next $x$ is finally uniformly distributed.
450 Thus Algorithm~\ref{algo:gpu_kernel2} is a concrete realization of the general
451 chaotic iterations presented previously, and for this reason, it satisfies the
452 Devaney's formulation of a chaotic behavior.
454 \section{Experiments}
455 \label{sec:experiments}
457 Different experiments have been performed in order to measure the generation
458 speed. We have used a first computer equipped with a Tesla C1060 NVidia GPU card
460 Intel Xeon E5530 cadenced at 2.40 GHz, and
461 a second computer equipped with a smaller CPU and a GeForce GTX 280.
463 cards have 240 cores.
465 In Figure~\ref{fig:time_xorlike_gpu} we compare the quantity of pseudorandom numbers
466 generated per second with various xor-like based PRNGs. In this figure, the optimized
467 versions use the {\it xor64} described in~\cite{Marsaglia2003}, whereas the naive versions
468 embed the three xor-like PRNGs described in Listing~\ref{algo:seqCIPRNG}. In
469 order to obtain the optimal performances, the storage of pseudorandom numbers
470 into the GPU memory has been removed. This step is time consuming and slows down the numbers
471 generation. Moreover this storage is completely
472 useless, in case of applications that consume the pseudorandom
473 numbers directly after generation. We can see that when the number of threads is greater
474 than approximately 30,000 and lower than 5 million, the number of pseudorandom numbers generated
475 per second is almost constant. With the naive version, this value ranges from 2.5 to
476 3GSamples/s. With the optimized version, it is approximately equal to
477 20GSamples/s. Finally we can remark that both GPU cards are quite similar, but in
478 practice, the Tesla C1060 has more memory than the GTX 280, and this memory
479 should be of better quality.
480 As a comparison, Listing~\ref{algo:seqCIPRNG} leads to the generation of about
481 138MSample/s when using one core of the Xeon E5530.
485 \includegraphics[scale=0.65]{Chapters/chapter18/figures/curve_time_xorlike_gpu.pdf}
487 \caption{Quantity of pseudorandom numbers generated per second with the xorlike-based PRNG}
488 \label{fig:time_xorlike_gpu}
495 In Figure~\ref{fig:time_bbs_gpu} we highlight the performances of the optimized
496 BBS-based PRNG on GPU. On the Tesla C1060 we obtain approximately 700MSample/s
497 and on the GTX 280 about 670MSample/s, which is obviously slower than the
498 xorlike-based PRNG on GPU. However, we will show in the next sections that this
499 new PRNG has a strong level of security, which is necessarily paid by a speed
504 \includegraphics[scale=0.65]{Chapters/chapter18/figures/curve_time_bbs_gpu.pdf}
506 \caption{Quantity of pseudorandom numbers generated per second using the BBS-based PRNG}
507 \label{fig:time_bbs_gpu}
510 All these experiments allow us to conclude that it is possible to
511 generate a very large quantity of pseudorandom numbers statistically perfect with the xor-like version.
512 To a certain extend, it is also the case with the secure BBS-based version, the speed deflation being
513 explained by the fact that the former version has ``only''
514 chaotic properties and statistical perfection, whereas the latter is also cryptographically secure,
515 as it is shown in the next sections.
527 \putbib[Chapters/chapter18/biblio18]