[prng_gpu.git] / prng_gpu.tex

\documentclass{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{fullpage}
\usepackage{fancybox}
\usepackage{amsmath}
\usepackage{amscd}
\usepackage{moreverb}
\usepackage{commath}
\usepackage[standard]{ntheorem}

% Pour mathds : les ensembles IR, IN, etc.
\usepackage{dsfont}

% Pour avoir des intervalles d'entiers
\usepackage{stmaryrd}

\usepackage{graphicx}
% Pour faire des sous-figures dans les figures
\usepackage{subfigure}

\usepackage{color}

\newtheorem{notation}{Notation}

\newcommand{\X}{\mathcal{X}}
\newcommand{\Go}{G_{f_0}}
\newcommand{\B}{\mathds{B}}
\newcommand{\N}{\mathds{N}}
\newcommand{\BN}{\mathds{B}^\mathsf{N}}
\let\sur=\overline

\newcommand{\alert}[1]{\begin{color}{blue}\textit{#1}\end{color}}

\title{Efficient generation of pseudo random numbers based on chaotic iterations on GPU}
\begin{document}
\maketitle

\begin{abstract}
This is the abstract
\end{abstract}

\section{Introduction}

Interet des itérations chaotiques pour générer des nombre alea\\
Interet de générer des nombres alea sur GPU
...

\section{Chaotic iterations}

Présentation des itérations chaotiques

\section{blabla}

\subsection{The phase space is an interval of the real line}

\subsubsection{Toward a topological semiconjugacy}

In what follows, our intention is to establish, by using a topological semiconjugacy, that chaotic iterations over $\mathcal{X}$ can be described as iterations on a real interval. To do so, we must firstly introduce some notations and terminologies. 

Let $\mathcal{S}_\mathsf{N}$ be the set of sequences belonging into $\llbracket 1; \mathsf{N}\rrbracket$ and $\mathcal{X}_{\mathsf{N}} = \mathcal{S}_\mathsf{N} \times \B^\mathsf{N}$.


\begin{definition}
The function $\varphi: \mathcal{S}_{10} \times\mathds{B}^{10} \rightarrow \big[ 0, 2^{10} \big[$ is defined by:
$$
\begin{array}{cccl}
\varphi: & \mathcal{X}_{10} = \mathcal{S}_{10} \times\mathds{B}^{10}& \longrightarrow & \big[ 0, 2^{10} \big[ \\
 & (S,E) = \left((S^0, S^1, \hdots ); (E_0, \hdots, E_9)\right) & \longmapsto & \varphi \left((S,E)\right)
\end{array}
$$
\noindent where $\varphi\left((S,E)\right)$ is the real number:
\begin{itemize}
\item whose integral part $e$ is $\displaystyle{\sum_{k=0}^9 2^{9-k} E_k}$, that is, the binary digits of $e$ are $E_0 ~ E_1 ~ \hdots ~ E_9$.
\item whose decimal part $s$ is equal to $s = 0,S^0~ S^1~ S^2~ \hdots = \sum_{k=1}^{+\infty} 10^{-k} S^{k-1}.$ 
\end{itemize}
\end{definition}



$\varphi$ realizes the association between a point of $\mathcal{X}_{10}$ and a real number into $\big[ 0, 2^{10} \big[$. We must now translate the chaotic iterations $\Go$ on this real interval. To do so, two intermediate functions over $\big[ 0, 2^{10} \big[$ must be introduced:


\begin{definition}
\label{def:e et s}
Let $x \in \big[ 0, 2^{10} \big[$ and:
\begin{itemize}
\item $e_0, \hdots, e_9$ the binary digits of the integral part of $x$: $\displaystyle{\lfloor x \rfloor = \sum_{k=0}^{9} 2^{9-k} e_k}$.
\item $(s^k)_{k\in \mathds{N}}$ the digits of $x$, where the chosen decimal decomposition of $x$ is the one that does not have an infinite number of 9: 
$\displaystyle{x = \lfloor x \rfloor + \sum_{k=0}^{+\infty} s^k 10^{-k-1}}$.
\end{itemize}
$e$ and $s$ are thus defined as follows:
$$
\begin{array}{cccl}
e: & \big[ 0, 2^{10} \big[ & \longrightarrow & \mathds{B}^{10} \\
 & x & \longmapsto & (e_0, \hdots, e_9)
\end{array}
$$
\noindent and
$$
\begin{array}{cccl}
s: & \big[ 0, 2^{10} \big[ & \longrightarrow & \llbracket 0, 9 \rrbracket^{\mathds{N}} \\
 & x & \longmapsto & (s^k)_{k \in \mathds{N}}
\end{array}
$$
\end{definition}

We are now able to define the function $g$, whose goal is to translate the chaotic iterations $\Go$ on an interval of $\mathds{R}$.

\begin{definition}
$g:\big[ 0, 2^{10} \big[ \longrightarrow \big[ 0, 2^{10} \big[$ is defined by:
$$
\begin{array}{cccl}
g: & \big[ 0, 2^{10} \big[ & \longrightarrow & \big[ 0, 2^{10} \big[ \\
& \\
 & x & \longmapsto & g(x)
\end{array}
$$
\noindent where g(x) is the real number of $\big[ 0, 2^{10} \big[$ defined bellow:
\begin{itemize}
\item its integral part has a binary decomposition equal to $e_0', \hdots, e_9'$, with:
$$
e_i' = \left\{
\begin{array}{ll}
e(x)_i & \textrm{ if } i \neq s^0\\
e(x)_i + 1 \textrm{ (mod 2)} & \textrm{ if } i = s^0\\
\end{array}
\right.
$$
\item whose decimal part is $s(x)^1, s(x)^2, \hdots$
\end{itemize}
\end{definition}

\bigskip


In other words, if $x = \displaystyle{\sum_{k=0}^{9} 2^{9-k} e_k +  \sum_{k=0}^{+\infty} s^{k} ~10^{-k-1}}$, then: $$g(x) = \displaystyle{\sum_{k=0}^{9} 2^{9-k} (e_k + \delta(k,s^0) \textrm{ (mod 2)}) +  \sum_{k=0}^{+\infty} s^{k+1} 10^{-k-1}}.$$

\subsubsection{Defining a metric on $\big[ 0, 2^{10} \big[$}

Numerous metrics can be defined on the set $\big[ 0, 2^{10} \big[$, the most usual one being the Euclidian distance recalled bellow:

\begin{notation}
\index{distance!euclidienne}
$\Delta$ is the Euclidian distance on $\big[ 0, 2^{10} \big[$, that is, $\Delta(x,y) = |y-x|^2$.
\end{notation}

\medskip

This Euclidian distance does not reproduce exactly the notion of proximity induced by our first distance $d$ on $\X$. Indeed $d$ is finer than $\Delta$. This is the reason why we have to introduce the following metric:



\begin{definition}
Let $x,y \in \big[ 0, 2^{10} \big[$.
$D$ denotes the function from $\big[ 0, 2^{10} \big[^2$ to $\mathds{R}^+$ defined by: $D(x,y) = D_e\left(e(x),e(y)\right) + D_s\left(s(x),s(y)\right)$, where:
\begin{center}
$\displaystyle{D_e(E,\check{E}) = \sum_{k=0}^\mathsf{9} \delta (E_k, \check{E}_k)}$, ~~and~ $\displaystyle{D_s(S,\check{S}) = \sum_{k = 1}^\infty \dfrac{|S^k-\check{S}^k|}{10^k}}$.
\end{center}
\end{definition}

\begin{proposition}
$D$ is a distance on $\big[ 0, 2^{10} \big[$.
\end{proposition}

\begin{proof}
The three axioms defining a distance must be checked.
\begin{itemize}
\item $D \geqslant 0$, because everything is positive in its definition. If $D(x,y)=0$, then $D_e(x,y)=0$, so the integral parts of $x$ and $y$ are equal (they have the same binary decomposition). Additionally, $D_s(x,y) = 0$, then $\forall k \in \mathds{N}^*, s(x)^k = s(y)^k$. In other words, $x$ and $y$ have the same $k-$th decimal digit, $\forall k \in \mathds{N}^*$. And so $x=y$.
\item $D(x,y)=D(y,x)$.
\item Finally, the triangular inequality is obtained due to the fact that both $\delta$ and $\Delta(x,y)=|x-y|$ satisfy it.
\end{itemize}
\end{proof}


The convergence of sequences according to $D$ is not the same than the usual convergence related to the Euclidian metric. For instance, if $x^n \to x$ according to $D$, then necessarily the integral part of each $x^n$ is equal to the integral part of $x$ (at least after a given threshold), and the decimal part of $x^n$ corresponds to the one of $x$ ``as far as required''.
To illustrate this fact, a comparison between $D$ and the Euclidian distance is given Figure \ref{fig:comparaison de distances}. These illustrations show that $D$ is richer and more refined than the Euclidian distance, and thus is more precise.


\begin{figure}[t]
\begin{center}
  \subfigure[Function $x \to dist(x;1,234) $ on the interval $(0;5)$.]{\includegraphics[scale=.35]{DvsEuclidien.pdf}}\quad
  \subfigure[Function $x \to dist(x;3) $ on the interval $(0;5)$.]{\includegraphics[scale=.35]{DvsEuclidien2.pdf}}
\end{center}
\caption{Comparison between $D$ (in blue) and the Euclidian distane (in green).}
\label{fig:comparaison de distances}
\end{figure}




\subsubsection{The semiconjugacy}

It is now possible to define a topological semiconjugacy between $\mathcal{X}$ and an interval of $\mathds{R}$:

\begin{theorem}
Chaotic iterations on the phase space $\mathcal{X}$ are simple iterations on $\mathds{R}$, which is illustrated by the semiconjugacy of the diagram bellow:
\begin{equation*}
\begin{CD}
\left(~\mathcal{S}_{10} \times\mathds{B}^{10}, d~\right) @>G_{f_0}>> \left(~\mathcal{S}_{10} \times\mathds{B}^{10}, d~\right)\\
    @V{\varphi}VV                    @VV{\varphi}V\\
\left( ~\big[ 0, 2^{10} \big[, D~\right)  @>>g> \left(~\big[ 0, 2^{10} \big[, D~\right)
\end{CD}
\end{equation*}
\end{theorem}

\begin{proof}
$\varphi$ has been constructed in order to be continuous and onto.
\end{proof}

In other words, $\mathcal{X}$ is approximately equal to $\big[ 0, 2^\mathsf{N} \big[$.






\subsection{Study of the chaotic iterations described as a real function}


\begin{figure}[t]
\begin{center}
  \subfigure[ICs on the interval $(0,9;1)$.]{\includegraphics[scale=.35]{ICs09a1.pdf}}\quad
  \subfigure[ICs on the interval $(0,7;1)$.]{\includegraphics[scale=.35]{ICs07a95.pdf}}\\
  \subfigure[ICs on the interval $(0,5;1)$.]{\includegraphics[scale=.35]{ICs05a1.pdf}}\quad
  \subfigure[ICs on the interval $(0;1)$]{\includegraphics[scale=.35]{ICs0a1.pdf}}
\end{center}
\caption{Representation of the chaotic iterations.}
\label{fig:ICs}
\end{figure}




\begin{figure}[t]
\begin{center}
  \subfigure[ICs on the interval $(510;514)$.]{\includegraphics[scale=.35]{ICs510a514.pdf}}\quad
  \subfigure[ICs on the interval $(1000;1008)$]{\includegraphics[scale=.35]{ICs1000a1008.pdf}}
\end{center}
\caption{ICs on small intervals.}
\label{fig:ICs2}
\end{figure}

\begin{figure}[t]
\begin{center}
  \subfigure[ICs on the interval $(0;16)$.]{\includegraphics[scale=.3]{ICs0a16.pdf}}\quad
  \subfigure[ICs on the interval  $(40;70)$.]{\includegraphics[scale=.45]{ICs40a70.pdf}}\quad
\end{center}
\caption{General aspect of the chaotic iterations.}
\label{fig:ICs3}
\end{figure}


We have written a Python program to represent the chaotic iterations with the vectorial negation on the real line $\mathds{R}$. Various representations of these CIs are given in Figures \ref{fig:ICs}, \ref{fig:ICs2} and \ref{fig:ICs3}. It can be remarked that the function $g$ is \alert{affine par morceaux}: it is linear on each interval having the form $\left[ \dfrac{n}{10}, \dfrac{n+1}{10}\right[$, $n \in \llbracket 0;2^{10}\times 10 \rrbracket$ \alert{and its line has a pent equal to 10}. Let us justify these claims:

\begin{proposition}
\label{Prop:derivabilite des ICs}
Chaotic iterations $g$ defined on $\mathds{R}$ are \alert{infiniment dérivables} on $\big[ 0, 2^{10} \big[$, except on the 10241 points in $I$ defined by $\left\{ \dfrac{n}{10} ~\big/~ n \in \llbracket 0;2^{10}\times 10\rrbracket \right\}$.

Furthermore, on each interval of the form $\left[ \dfrac{n}{10}, \dfrac{n+1}{10}\right[$, with $n \in \llbracket 0;2^{10}\times 10 \rrbracket$, the function $g$ is \emph{affine}. \alert{It is a line of pent equal to 10}: $\forall x \notin I, g'(x)=10$.
\end{proposition}


\begin{proof}
Let $I_n = \left[ \dfrac{n}{10}, \dfrac{n+1}{10}\right[$, with $n \in \llbracket 0;2^{10}\times 10 \rrbracket$. All the points of $I_n$ have the same integral prat $e$ and the same decimal part $s^0$: on the set $I_n$,  functions $e(x)$ and $x \mapsto s(x)^0$ of Definition \ref{def:e et s} only depend on $n$. So all the images $g(x)$ of these points $x$:
\begin{itemize}
\item Have the same integral part, which is $e$, except probably the bit number $s^0$. In other words, this integer has approximately the same binary decomposition than $e$, the sole exception being the digit $s^0$ (this number is then either $e+2^{10-s^0}$ or $e-2^{10-s^0}$, depending on the parity of $s^0$, \emph{i.e.}, it is equal to $e+(-1)^{s^0}\times 2^{10-s^0}$).
\item A shift to the left has been applied to the decimal part $y$, losing by doing so the common first digit $s^0$. In other words, $y$ has been mapped into $10\times y - s^0$.
\end{itemize}
To sum up, the action of $g$ on the points of $I$ is as follows: first, make a multiplication by 10, and second, add the same constant to each term, which is $\dfrac{1}{10}\left(e+(-1)^{s^0}\times 2^{10-s^0}\right)-s^0$.
\end{proof}

\begin{remark}
Finally, chaotic iterations are elements of the large family of functions that are \alert{chaotiques linéaires par morceaux}, like the tent map, the \alert{doublement de l'angle}, \emph{etc.}
\end{remark}



\subsection{Comparison of the two metrics on $\big[ 0, 2^\mathsf{N} \big[$}

The two propositions bellow allow to compare our two distances on $\big[ 0, 2^\mathsf{N} \big[$:

\begin{proposition}
Id: $\left(~\big[ 0, 2^\mathsf{N} \big[,\Delta~\right) \to \left(~\big[ 0, 2^\mathsf{N} \big[, D~\right)$ is not continuous. 
\end{proposition}

\begin{proof}
The sequence $x^n = 1,999\hdots 999$ constituted by $n$ 9 as decimal part, is such that:
\begin{itemize}
\item $\Delta (x^n,2) \to 0.$
\item But $D(x^n,2) \geqslant 1$, then $D(x^n,2)$ does not converge to 0.
\end{itemize}

The sequential characterization of the continuity concludes the demonstration.
\end{proof}



A contrario:

\begin{proposition}
Id: $\left(~\big[ 0, 2^\mathsf{N} \big[,D~\right) \to \left(~\big[ 0, 2^\mathsf{N} \big[, \Delta ~\right)$ is a continuous fonction. 
\end{proposition}

\begin{proof}
If $D(x^n,x) \to 0$, then $D_e(x^n,x) = 0$ at least for $n$ larger than a given threshold, because $D_e$ only returns integers. So, after this threshold, the integral parts of all the $x^n$ are equal to the integral part of $x$. 

Additionally, $D_s(x^n, x) \to 0$, then $\forall k \in \mathds{N}^*, \exists N_k \in \mathds{N}, n \geqslant N_k \Rightarrow D_s(x^n,x) \leqslant 10^{-k}$. This means that for all $k$, an index $N_k$ can be found such that, $\forall n \geqslant N_k$, all the $x^n$ have the same $k$ firsts digits, which are the digits of $x$. We can deduce the convergence $\Delta(x^n,x) \to 0$, and thus the result.
\end{proof}

The conclusion of these propositions is that the proposed metric is more \alert{précise} than the Euclidian distance, that is:

\begin{corollary}
$D$ is finer than the Euclidian distance $\Delta$.
\end{corollary}

This corollary can be reformulated as follows:

\begin{itemize}
\item The topology produced by $\Delta$ is a subset of the topology produced by $D$.
\item $D$ has more open sets than $\Delta$.
\item It is harder to converge for the topology $\tau_D$ inherited by $D$, than to converge with the one inherited by $\Delta$, which is denoted here by $\tau_\Delta$.
\end{itemize}


\subsection{Chaos of the chaotic iterations on $\mathds{R}$}
\label{chpt:Chaos des itérations chaotiques sur R}



\subsubsection{Chaos according to Devaney}

We have recalled previously that the chaotic iterations $\left(\Go, \mathcal{X}_d\right)$ are chaotic according to the formulation of Devaney. We can deduce that they are chaotic on $\mathds{R}$ too, when considering the \alert{topology of order}, because:
\begin{itemize}
\item $\left(\Go, \mathcal{X}_d\right)$ and $\left(g, \big[ 0, 2^{10} \big[_D\right)$ are semiconjugate by $\varphi$,
\item Then $\left(g, \big[ 0, 2^{10} \big[_D\right)$ is a system chaotic according to Devaney, because the semiconjugacy preserve this character.
\item But the topology generated by $D$ is finer than the topology generated by the Euclidian distance $\Delta$ -- which is the \alert{topology of order}.
\item According to Theorem \ref{Th:chaos et finesse}, we can deduce that the chaotic iterations $g$ are indeed chaotic, as defined by Devaney, for the \alert{topology of order} on $\mathds{R}$.
\end{itemize}

This result can be formulated as follows.

\begin{theorem}
\label{th:IC et topologie de l'ordre}
The chaotic iterations $g$ on $\mathds{R}$ are chaotic according to the Devaney's formulation, when $\mathds{R}$ has his usual topology, which is the \alert{topology of order}.
\end{theorem}

Indeed this result is \alert{weaker} than the theorem establishing the chaos for the finer topology $d$. However the Theorem \ref{th:IC et topologie de l'ordre} still remains important. Indeed, we have studied in our previous works a set different from the usual set of study ($\mathcal{X}$ instead of $\mathds{R}$), in order to be as close as possible from the computer: the properties of disorder proved theoretically will then be preserved when computing. However, we could wonder whether this change does not lead to a disorder of a lower quality. In other words, have we replaced a situation of a good disorder lost when computing, to another situation of a disorder preserved but of bad quality. Theorem \ref{th:IC et topologie de l'ordre} prove exactly the contrary.
 



\section{Efficient prng based on chaotic iterations}

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
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
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}
}
\end{center}
\caption{sequential Chaotic Iteration PRNG}
\label{algo:seqCIprng}
\end{figure}

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
  R. Simard. Testu01].

\section{Efficient prng based on chaotic iterations on GPU}

On parle du passage du sequentiel au GPU

\section{Experiments}

On passe le BigCrush\\
On donne des temps de générations sur GPU/CPU\\
On donne des temps de générations de nombre sur GPU puis on rappatrie sur CPU / CPU ? bof bof, on verra


\section{Conclusion}
\bibliographystyle{plain}
\bibliography{mabase}
\end{document}