X-Git-Url: https://bilbo.iut-bm.univ-fcomte.fr/and/gitweb/prng_gpu.git/blobdiff_plain/12410996f49b94978103d0ba32d17f85c89fa51a..a595bc795f31e05fc7fcc8415e1549bdda84b076:/prng_gpu.tex

diff --git a/prng_gpu.tex b/prng_gpu.tex
index 082c379..730b620 100644
--- a/prng_gpu.tex
+++ b/prng_gpu.tex
@@ -7,6 +7,8 @@
 \usepackage{amscd}
 \usepackage{moreverb}
 \usepackage{commath}
+\usepackage{algorithm2e}
+\usepackage{listings}
 \usepackage[standard]{ntheorem}
 
 % Pour mathds : les ensembles IR, IN, etc.
@@ -34,6 +36,9 @@
 
 \title{Efficient generation of pseudo random numbers based on chaotic iterations on GPU}
 \begin{document}
+
+\author{Jacques M. Bahi, Rapha\"{e}l Couturier, and Christophe Guyeux\thanks{Authors in alphabetic order}}
+
 \maketitle
 
 \begin{abstract}
@@ -44,12 +49,245 @@ This is the abstract
 
 Interet des itérations chaotiques pour générer des nombre alea\\
 Interet de générer des nombres alea sur GPU
+\alert{RC, un petit state-of-the-art sur les PRNGs sur GPU ?}
 ...
 
-\section{Chaotic iterations}
 
-Présentation des itérations chaotiques
+\section{Basic Recalls}
+\label{section:BASIC RECALLS}
+This section is devoted to basic definitions and terminologies in the fields of topological chaos and chaotic iterations.
+\subsection{Devaney's chaotic dynamical systems}
+
+In the sequel $S^{n}$ denotes the $n^{th}$ term of a sequence $S$ and $V_{i}$ denotes the $i^{th}$ component of a vector $V$. $f^{k}=f\circ ...\circ f$ denotes the $k^{th}$ composition of a function $f$. Finally, the following notation is used: $\llbracket1;N\rrbracket=\{1,2,\hdots,N\}$.
+
+
+Consider a topological space $(\mathcal{X},\tau)$ and a continuous function $f : \mathcal{X} \rightarrow \mathcal{X}$.
+
+\begin{definition}
+$f$ is said to be \emph{topologically transitive} if, for any pair of open sets $U,V \subset \mathcal{X}$, there exists $k>0$ such that $f^k(U) \cap V \neq \varnothing$.
+\end{definition}
+
+\begin{definition}
+An element $x$ is a \emph{periodic point} for $f$ of period $n\in \mathds{N}^*$ if $f^{n}(x)=x$.% The set of periodic points of $f$ is denoted $Per(f).$
+\end{definition}
+
+\begin{definition}
+$f$ is said to be \emph{regular} on $(\mathcal{X}, \tau)$ if the set of periodic points for $f$ is dense in $\mathcal{X}$: for any point $x$ in $\mathcal{X}$, any neighborhood of $x$ contains at least one periodic point (without necessarily the same period).
+\end{definition}
+
+
+\begin{definition}
+$f$ is said to be \emph{chaotic} on $(\mathcal{X},\tau)$ if $f$ is regular and topologically transitive.
+\end{definition}
+
+The chaos property is strongly linked to the notion of ``sensitivity'', defined on a metric space $(\mathcal{X},d)$ by:
+
+\begin{definition}
+\label{sensitivity} $f$ has \emph{sensitive dependence on initial conditions}
+if there exists $\delta >0$ such that, for any $x\in \mathcal{X}$ and any neighborhood $V$ of $x$, there exist $y\in V$ and $n > 0$ such that $d\left(f^{n}(x), f^{n}(y)\right) >\delta $.
+
+$\delta$ is called the \emph{constant of sensitivity} of $f$.
+\end{definition}
 
+Indeed, Banks \emph{et al.} have proven in~\cite{Banks92} that when $f$ is chaotic and $(\mathcal{X}, d)$ is a metric space, then $f$ has the property of sensitive dependence on initial conditions (this property was formerly an element of the definition of chaos). To sum up, quoting Devaney in~\cite{Devaney}, a chaotic dynamical system ``is unpredictable because of the sensitive dependence on initial conditions. It cannot be broken down or simplified 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.
+
+
+
+\subsection{Chaotic iterations}
+\label{sec:chaotic iterations}
+
+
+Let us consider  a \emph{system} with a finite  number $\mathsf{N} \in
+\mathds{N}^*$ of elements  (or \emph{cells}), so that each  cell has a
+Boolean  \emph{state}. Having $\mathsf{N}$ Boolean values for these
+ cells  leads to the definition of a particular \emph{state  of the
+system}. A sequence which  elements belong to $\llbracket 1;\mathsf{N}
+\rrbracket $ is called a \emph{strategy}. The set of all strategies is
+denoted by $\mathbb{S}.$
+
+\begin{definition}
+\label{Def:chaotic iterations}
+The      set       $\mathds{B}$      denoting      $\{0,1\}$,      let
+$f:\mathds{B}^{\mathsf{N}}\longrightarrow  \mathds{B}^{\mathsf{N}}$ be
+a  function  and  $S\in  \mathbb{S}$  be  a  strategy.  The  so-called
+\emph{chaotic      iterations}     are     defined      by     $x^0\in
+\mathds{B}^{\mathsf{N}}$ and
+$$
+\forall    n\in     \mathds{N}^{\ast     },    \forall     i\in
+\llbracket1;\mathsf{N}\rrbracket ,x_i^n=\left\{
+\begin{array}{ll}
+  x_i^{n-1} &  \text{ if  }S^n\neq i \\
+  \left(f(x^{n-1})\right)_{S^n} & \text{ if }S^n=i.
+\end{array}\right.
+$$
+\end{definition}
+
+In other words, at the $n^{th}$ iteration, only the $S^{n}-$th cell is
+\textquotedblleft  iterated\textquotedblright .  Note  that in  a more
+general  formulation,  $S^n$  can   be  a  subset  of  components  and
+$\left(f(x^{n-1})\right)_{S^{n}}$      can     be      replaced     by
+$\left(f(x^{k})\right)_{S^{n}}$, where  $k<n$, describing for example,
+delays  transmission~\cite{Robert1986,guyeux10}.  Finally,  let us  remark that
+the term  ``chaotic'', in  the name of  these iterations,  has \emph{a
+priori} no link with the mathematical theory of chaos, recalled above.
+
+
+Let us now recall how to define a suitable metric space where chaotic iterations are continuous. For further explanations, see, e.g., \cite{guyeux10}.
+
+Let $\delta $ be the \emph{discrete Boolean metric}, $\delta (x,y)=0\Leftrightarrow x=y.$ Given a function $f$, define the function:
+\begin{equation*}
+\begin{array}{lrll}
+F_{f}: & \llbracket1;\mathsf{N}\rrbracket\times \mathds{B}^{\mathsf{N}} &
+\longrightarrow & \mathds{B}^{\mathsf{N}} \\
+& (k,E) & \longmapsto & \left( E_{j}.\delta (k,j)+f(E)_{k}.\overline{\delta
+(k,j)}\right) _{j\in \llbracket1;\mathsf{N}\rrbracket},%
+\end{array}%
+\end{equation*}%
+\noindent where + and . are the Boolean addition and product operations.
+Consider the phase space:
+\begin{equation*}
+\mathcal{X} = \llbracket 1 ; \mathsf{N} \rrbracket^\mathds{N} \times
+\mathds{B}^\mathsf{N},
+\end{equation*}
+\noindent and the map defined on $\mathcal{X}$:
+\begin{equation}
+G_f\left(S,E\right) = \left(\sigma(S), F_f(i(S),E)\right), \label{Gf}
+\end{equation}
+\noindent where $\sigma$ is the \emph{shift} function defined by $\sigma (S^{n})_{n\in \mathds{N}}\in \mathbb{S}\longrightarrow (S^{n+1})_{n\in \mathds{N}}\in \mathbb{S}$ and $i$ is the \emph{initial function}  $i:(S^{n})_{n\in \mathds{N}} \in \mathbb{S}\longrightarrow S^{0}\in \llbracket 1;\mathsf{N}\rrbracket$. Then the chaotic iterations defined in (\ref{sec:chaotic iterations}) can be described by the following iterations:
+\begin{equation*}
+\left\{
+\begin{array}{l}
+X^0 \in \mathcal{X} \\
+X^{k+1}=G_{f}(X^k).%
+\end{array}%
+\right.
+\end{equation*}%
+
+With this formulation, a shift function appears as a component of chaotic iterations. The shift function is a famous example of a chaotic map~\cite{Devaney} but its presence is not sufficient enough to claim $G_f$ as chaotic. 
+
+To study this claim, a new distance between two points $X = (S,E), Y = (\check{S},\check{E})\in
+\mathcal{X}$ has been introduced in \cite{guyeux10} as follows:
+\begin{equation*}
+d(X,Y)=d_{e}(E,\check{E})+d_{s}(S,\check{S}),
+\end{equation*}
+\noindent where
+\begin{equation*}
+\left\{
+\begin{array}{lll}
+\displaystyle{d_{e}(E,\check{E})} & = & \displaystyle{\sum_{k=1}^{\mathsf{N}%
+}\delta (E_{k},\check{E}_{k})}, \\
+\displaystyle{d_{s}(S,\check{S})} & = & \displaystyle{\dfrac{9}{\mathsf{N}}%
+\sum_{k=1}^{\infty }\dfrac{|S^k-\check{S}^k|}{10^{k}}}.%
+\end{array}%
+\right.
+\end{equation*}
+
+
+This new distance has been introduced to satisfy the following requirements.
+\begin{itemize}
+\item When the number of different cells between two systems is increasing, then their distance should increase too.
+\item In addition, if two systems present the same cells and their respective strategies start with the same terms, then the distance between these two points must be small because the evolution of the two systems will be the same for a while. Indeed, the two dynamical systems start with the same initial condition, use the same update function, and as strategies are the same for a while, then components that are updated are the same too.
+\end{itemize}
+The distance presented above follows these recommendations. Indeed, if the floor value $\lfloor d(X,Y)\rfloor $ is equal to $n$, then the systems $E, \check{E}$ differ in $n$ cells. In addition, $d(X,Y) - \lfloor d(X,Y) \rfloor $ is a measure of the differences between strategies $S$ and $\check{S}$. More precisely, this floating part is less than $10^{-k}$ if and only if the first $k$ terms of the two strategies are equal. Moreover, if the $k^{th}$ digit is nonzero, then the $k^{th}$ terms of the two strategies are different.
+
+Finally, it has been established in \cite{guyeux10} that,
+
+\begin{proposition}
+Let $f$ be a map from $\mathds{B}^n$ to itself. Then $G_{f}$ is continuous in the metric space $(\mathcal{X},d)$.
+\end{proposition}
+
+The chaotic property of $G_f$ has been firstly established for the vectorial Boolean negation \cite{guyeux10}. To obtain a characterization, we have secondly introduced the notion of asynchronous iteration graph recalled bellow.
+
+Let $f$ be a map from $\mathds{B}^n$ to itself. The
+{\emph{asynchronous iteration graph}} associated with $f$ is the
+directed graph $\Gamma(f)$ defined by: the set of vertices is
+$\mathds{B}^n$; for all $x\in\mathds{B}^n$ and $i\in \llbracket1;n\rrbracket$,
+the graph $\Gamma(f)$ contains an arc from $x$ to $F_f(i,x)$. 
+The relation between $\Gamma(f)$ and $G_f$ is clear: there exists a
+path from $x$ to $x'$ in $\Gamma(f)$ if and only if there exists a
+strategy $s$ such that the parallel iteration of $G_f$ from the
+initial point $(s,x)$ reaches the point $x'$.
+
+We have finally proven in \cite{FCT11} that,
+
+
+\begin{theorem}
+\label{Th:Caractérisation   des   IC   chaotiques}  
+Let $f:\mathds{B}^n\to\mathds{B}^n$. $G_f$ is chaotic  (according to  Devaney) 
+if and only if $\Gamma(f)$ is strongly connected.
+\end{theorem}
+
+This result of chaos has lead us to study the possibility to build a pseudo-random number generator (PRNG) based on the chaotic iterations. 
+As $G_f$, defined on the domain   $\llbracket 1 ;  n \rrbracket^{\mathds{N}}  \times \mathds{B}^n$, is build from Boolean networks $f : \mathds{B}^n \rightarrow \mathds{B}^n$, we can preserve the theoretical properties on $G_f$ during implementations (due to the discrete nature of $f$). It is as if $\mathds{B}^n$ represents the memory of the computer whereas $\llbracket 1 ;  n \rrbracket^{\mathds{N}}$ is its input stream (the seeds, for instance).
+
+\section{Application to Pseudo-Randomness}
+
+We have proposed in~\cite{bgw09:ip} a new family of generators that receives 
+two PRNGs as inputs. These two generators are mixed with chaotic iterations, 
+leading thus to a new PRNG that improves the statistical properties of each
+generator taken alone. Furthermore, our generator 
+possesses various chaos properties
+that none of the generators used as input present.
+
+\begin{algorithm}[h!]
+%\begin{scriptsize}
+\KwIn{a function $f$, an iteration number $b$, an initial configuration $x^0$ ($n$ bits)}
+\KwOut{a configuration $x$ ($n$ bits)}
+$x\leftarrow x^0$\;
+$k\leftarrow b + \textit{XORshift}(b+1)$\;
+\For{$i=0,\dots,k-1$}
+{
+$s\leftarrow{\textit{XORshift}(n)}$\;
+$x\leftarrow{F_f(s,x)}$\;
+}
+return $x$\;
+%\end{scriptsize}
+\caption{PRNG with chaotic functions}
+\label{CI Algorithm}
+\end{algorithm}
+
+\begin{algorithm}[h!]
+%\SetAlgoLined                        %%RAPH: cette ligne provoque une erreur chez moi
+\KwIn{the internal configuration $z$ (a 32-bit word)}
+\KwOut{$y$ (a 32-bit word)}
+$z\leftarrow{z\oplus{(z\ll13)}}$\;
+$z\leftarrow{z\oplus{(z\gg17)}}$\;
+$z\leftarrow{z\oplus{(z\ll5)}}$\;
+$y\leftarrow{z}$\;
+return $y$\;
+\medskip
+\caption{An arbitrary round of \textit{XORshift} algorithm}
+\label{XORshift}
+\end{algorithm}
+
+
+
+
+
+This generator is synthesized in Algorithm~\ref{CI Algorithm}.
+It takes as input: a function $f$;
+an integer $b$, ensuring that the number of executed iterations is at least $b$ and at most $2b+1$; and an initial configuration $x^0$.
+It returns the new generated configuration $x$.  Internally, it embeds two
+\textit{XORshift}$(k)$ PRNGs \cite{Marsaglia2003} that returns integers uniformly distributed
+into $\llbracket 1 ; k \rrbracket$.
+\textit{XORshift} is a category of very fast PRNGs designed by George Marsaglia, which repeatedly uses the transform of exclusive or (XOR, $\oplus$) on a number with a bit shifted version of it. This PRNG, which has a period of $2^{32}-1=4.29\times10^9$, is summed up in Algorithm~\ref{XORshift}. It is used in our PRNG to compute the strategy length and the strategy elements.
+
+
+We have proven in \cite{FCT11} that,
+
+\begin{theorem}
+  Let $f: \mathds{B}^{n} \rightarrow \mathds{B}^{n}$, $\Gamma(f)$ its
+  iteration graph, $\check{M}$ its adjacency
+  matrix and $M$ a $n\times n$ matrix defined as in the previous lemma.
+  If $\Gamma(f)$ is strongly connected, then 
+  the output of the PRNG detailed in Algorithm~\ref{CI Algorithm} follows 
+  a law that tends to the uniform distribution 
+  if and only if $M$ is a double stochastic matrix.
+\end{theorem} 
+
+
+
+\alert{Mettre encore un peu de blabla sur le PRNG, puis enchaîner en disant que, ok, on peut préserver le chaos quand on passe sur machine, mais que le chaos dont il s'agit a été prouvé pour une distance bizarroïde sur un espace non moins hémoroïde, d'où ce qui suit}
 
 
 
@@ -431,52 +669,121 @@ Indeed this result is weaker than the theorem establishing the chaos for the fin
 
 On parle du séquentiel avec des nombres 64 bits\\
 
-Faire le lien avec le paragraphe précédent (je considère que la stratégie s'appelle $S^i$\\
 
 In  order to  implement efficiently  a PRNG  based on  chaotic iterations  it is
 possible to improve  previous works [ref]. One solution  consists in considering
-that the  strategy used $S^i$  contains all the  bits for which the  negation is
+that the  strategy used contains all the  bits for which the  negation is
 achieved out. Then instead of applying  the negation on these bits we can simply
-apply the  xor operator between  the current number  and the strategy  $S^i$. In
+apply the  xor operator between  the current number  and the strategy. In
 order to obtain the strategy we also use a classical PRNG.
 
-\begin{figure}[htbp]
-\begin{center}
-\fbox{
-\begin{minipage}{14cm}
-unsigned int CIprng() \{\\
-  static unsigned int x = 123123123;\\
-  unsigned long t1 = xorshift();\\
-  unsigned long t2 = xor128();\\
-  unsigned long t3 = xorwow();\\
-  x = x\textasciicircum (unsigned int)t1;\\
-  x = x\textasciicircum (unsigned int)(t2$>>$32);\\
-  x = x\textasciicircum (unsigned int)(t3$>>$32);\\
-  x = x\textasciicircum (unsigned int)t2;\\
-  x = x\textasciicircum (unsigned int)(t1$>>$32);\\
-  x = x\textasciicircum (unsigned int)t3;\\
-  return x;\\
-\}
-\end{minipage}
+%% \begin{figure}[htbp]
+%% \begin{center}
+%% \fbox{
+%% \begin{minipage}{14cm}
+%% unsigned int CIprng() \{\\
+%%   static unsigned int x = 123123123;\\
+%%   unsigned long t1 = xorshift();\\
+%%   unsigned long t2 = xor128();\\
+%%   unsigned long t3 = xorwow();\\
+%%   x = x\textasciicircum (unsigned int)t1;\\
+%%   x = x\textasciicircum (unsigned int)(t2$>>$32);\\
+%%   x = x\textasciicircum (unsigned int)(t3$>>$32);\\
+%%   x = x\textasciicircum (unsigned int)t2;\\
+%%   x = x\textasciicircum (unsigned int)(t1$>>$32);\\
+%%   x = x\textasciicircum (unsigned int)t3;\\
+%%   return x;\\
+%% \}
+%% \end{minipage}
+%% }
+%% \end{center}
+%% \caption{sequential Chaotic Iteration PRNG}
+%% \label{algo:seqCIprng}
+%% \end{figure}
+
+
+
+\lstset{language=C,caption={C code of the sequential chaotic iterations based PRNG},label=algo:seqCIprng}
+\begin{lstlisting}
+unsigned int CIprng() {
+  static unsigned int x = 123123123;
+  unsigned long t1 = xorshift();
+  unsigned long t2 = xor128();
+  unsigned long t3 = xorwow();
+  x = x^(unsigned int)t1;
+  x = x^(unsigned int)(t2>>32);
+  x = x^(unsigned int)(t3>>32);
+  x = x^(unsigned int)t2;
+  x = x^(unsigned int)(t1>>32);
+  x = x^(unsigned int)t3;
+  return x;
 }
-\end{center}
-\caption{sequential Chaotic Iteration PRNG}
-\label{algo:seqCIprng}
-\end{figure}
+\end{lstlisting}
+
 
-In Figure~\ref{algo:seqCIprng}  a sequential  version of our  chaotic iterations
-based PRNG  is presented.  This version  uses three classical 64  bits PRNG: the
-\texttt{xorshift},  the \texttt{xor128}  and the  \texttt{xorwow}.   These three
-PRNGs  are presented  in~\cite{Marsaglia2003}.   As each  PRNG  used works  with
-64-bits and as  our PRNG works with 32 bits, the  use of \texttt{(unsigned int)}
-selects the 32 least  significant bits whereas \texttt{(unsigned int)(t3$>>$32)}
-selects the 32  most significants bits of the  variable \texttt{t}. This version
-sucesses   the   BigCrush   of    the   TestU01   battery   [P.   L’ecuyer   and
+
+
+
+In listing~\ref{algo:seqCIprng}  a sequential  version of our  chaotic iterations
+based PRNG  is presented.  This version  uses three classical  64-bits PRNG: the
+\texttt{xorshift},   the  \texttt{xor128}  and   the  \texttt{xorwow}.   In  the
+following,  we  call them  xor-like  PRNGSs.   These  three PRNGs  are  presented
+in~\cite{Marsaglia2003}.  As each xor-like  PRNG used works with 64-bits  and as our PRNG
+works  with 32-bits, the  use of  \texttt{(unsigned int)}  selects the  32 least
+significant bits  whereas \texttt{(unsigned int)(t3$>>$32)} selects  the 32 most
+significants  bits of the  variable \texttt{t}.  So to  produce a  random number
+realizes 6  xor operations with  6 32-bits numbers  produced by 3  64-bits PRNG.
+This version  successes the  BigCrush of the  TestU01 battery [P.   L’ecuyer and
   R. Simard. Testu01].
 
 \section{Efficient prng based on chaotic iterations on GPU}
 
-On parle du passage du sequentiel au GPU
+In  order to benefit  from computing  power of  GPU, a  program needs  to define
+independent blocks of threads which  can be computed simultaneously. In general,
+the larger the number of threads is,  the more local memory is used and the less
+branching  instructions are  used (if,  while, ...),  the better  performance is
+obtained  on  GPU.  So  with  algorithm  \ref{algo:seqCIprng}  presented in  the
+previous section, it is possible to  build a similar program which computes PRNG
+on  GPU. The principe  consists in  assigning the  computation of  a PRNG  as in
+sequential to each thread of the GPU.  Of course, it is essential that the three
+xor-like PRNGs used  for our computation have different  parameters. So we chose
+them randomly with another PRNG. As  the initialisation is performed by the CPU,
+we have chosen to use the ISAAC PRNG to initalize all the parameters for the GPU
+version  of  our  PRNG.  The  implementation  of  the  three xor-like  PRNGs  is
+straightforward  as soon  as their  parameters have  been allocated  in  the GPU
+memory. Each xor-like  PRNGs used works with an internal  number $x$ which keeps
+the last generated random numbers. Other internal variables are also used by the
+xor-like PRNGs. More  precisely, the implementation of the  xor128, the xorshift
+and  the xorwow  respectively  require 4,  5  and 6  unsigned  long as  internal
+variables.
+
+\begin{algorithm}
+
+\KwIn{InternalVarXorLikeArray: array with internal variables of the 3 xor-like PRNGs in global memory\;
+NumThreads: Number of threads\;}
+\KwOut{NewNb: array containing random numbers in global memory}
+\If{threadId is concerned} {
+  retrieve data from InternalVarXorLikeArray[threadId] in local variables\;
+  \For{i=1 to n} {
+    compute a new PRNG as in Listing\ref{algo:seqCIprng}\;
+    store the new PRNG in NewNb[NumThreads*threadId+i]\;
+  }
+  store internal variables in InternalVarXorLikeArray[threadId]\;
+}
+
+\caption{main kernel for the chaotic iterations based PRNG GPU version}
+\label{algo:gpu_kernel}
+\end{algorithm}
+
+According to  the available  memory in the  GPU and  the number of  threads used
+simultenaously, the number of random numbers that a thread can generate inside a
+kernel     is     limited,      i.e.      the     variable     \texttt{n}     in
+algorithm~\ref{algo:gpu_kernel}. For example, if  $100,000$ threads are used and
+if $n=100$\footnote{in fact, we need to add the initial seed (a 32-bits number)}
+then   the  memory   required   to  store   internals   variables  of   xor-like
+PRNGs\footnote{we multiply this number by $2$ in order to count 32-bits numbers}
+and  random  number of  our  PRNG  is  equals to  $100,000\times  ((4+5+6)\times
+2+(1+100))=1,310,000$ 32-bits numbers, i.e. about $52$Mb.
 
 \section{Experiments}