ajout de cr

[14Secrypt.git] / main.tex

\documentclass[a4paper,twoside]{article}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc} 
\usepackage[english]{babel}
%\usepackage{ntheorem}
\usepackage{alltt}
\usepackage{amsmath,amssymb,amsthm,amstext,latexsym,eufrak,euscript}
\usepackage{subfigure,pstricks,pst-node,pst-coil}
\usepackage{dsfont}
\usepackage{stmaryrd}
\usepackage{cite}
\usepackage{graphicx}
\usepackage{algorithm2e}
\usepackage[font=footnotesize]{subfig}
\usepackage{listings}
\usepackage{booktabs}
\usepackage{epsfig}
\usepackage{url}
\usepackage{calc}
\usepackage{multicol}
\usepackage{pslatex}
\usepackage{apalike}
\usepackage{SCITEPRESS}
\usepackage{caption}



\subfigtopskip=0pt
\subfigcapskip=0pt
\subfigbottomskip=0pt







\newcommand{\JFC}[1]{\begin{color}{green}\textit{}\end{color}}
\newcommand{\CG}[1]{\begin{color}{blue}\textit{}\end{color}}
\newcommand{\tv}[1] {\lVert #1 \rVert_{\rm TV}}


\newtheorem*{xpl}{Running example}
\newtheorem{theorem}{Theorem}
 



\newcommand{\vectornorm}[1]{\ensuremath{\left|\left|#1\right|\right|_2}}
\newcommand{\ie}{\textit{i.e.}}
\newcommand{\Nats}[0]{\ensuremath{\mathbb{N}}}
\newcommand{\Reels}[0]{\ensuremath{\mathbb{R}}}
\newcommand{\Bool}[0]{\ensuremath{\mathds{B}}}
\newcommand{\rel}[0]{\ensuremath{{\mathcal{R}}}}
\newcommand{\Gall}[0]{\ensuremath{\mathcal{G}}}

\begin{document}

\title{Pseudorandom Number Generators with Balanced Gray Codes}


%\author{\authorname{J.-F. Couchot\sup{1}, P.-C. Heam\sup{1}, C. Guyeux\sup{1}, Q. Wang\sup{2},  and J. M. Bahi\sup{1}}
%\affiliation{\sup{1}FEMTO-ST Institute, University of Franche-Comt\'{e}, France}
%\affiliation{\sup{2}College of Automation, Guangdong University of Technology, China }
%\email{\{jean-francois.couchot,pierre-cyrille.heam,christophe.guyeux,jacques.bahi\}@univ-fcomte.fr, wangqianxue@gdut.edu.cn}
%}

\keywords{Pseudorandom Number Generator, Theory of Chaos, Markov Matrices, Hamiltonian Paths, Mixing Time, Gray Codes, Statistical Tests}

\abstract{In this article,
it is shown that a large class of truly chaotic Pseudorandom Number Generators
can be constructed. The generators are based on iterating Boolean maps, which are computed 
using balanced Gray codes. The number of such Gray codes gives the size of the class.
The construction of such generators is automatic for small number of bits, but remains an open problem when this 
number becomes large. A running example is used throughout the paper. Finally, first statistical experiments of these generators are presented, they show how efficient and promising the proposed
approach seems.}


\onecolumn \maketitle \normalsize \vfill






\section{\uppercase{Introduction}}

Many fields of research or applications like numerical simulations, 
stochastic optimization, or information security are highly dependent on
the use of fast and unbiased random number generators. Depending on the
targeted application, reproducibility must be either required, leading
to deterministic algorithms that produce numbers as close as possible
to random sequences, or refused, which implies to use an external
physical noise. The formers are called pseudorandom number generators
(PRNGs) while the laters are designed by truly random number generators (TRNGs).
TRNGs are used for instance in cypher keys generation, or in hardware based
simulations or security devices. Such TRNGs are often based on a chaotic
physical signal, may be quantized depending on the application.
This quantization however raises the problem of the degradation of chaotic 
properties.

The use of PRNGs, for its part, is a necessity in a large variety 
of numerical simulations, in which
responses of devices under study must be compared using the same ``random''
stream. This reproducibility is required too for symmetric encryption like
one-time pad, as sender and receiver must share the same pad. However, in
that situation, security of the pseudorandom stream must be mathematically
proven: an attacker must not be able to computationally distinguish a
pseudorandom sequence generated by the considered PRNG with a really random
one. Such cryptographically secure pseudorandom number generators are
however only useful in cryptographic contexts, due to their slowness
resulting from their security. 


Other kind of properties are desired for PRNGs used in numerical simulations
or in programs that embed a Monte-Carlo algorithm. In these situations,
required properties are speed and random-like profiles of the generated 
sequences. The fact that a given PRNG is unbiased and behaves as a white
noise is thus verified using batteries of statistical tests on a large amount
of pseudorandom numbers. Reputed and up-to-date batteries are
currently the NIST suite~\cite{Nist10}, DieHARD~\cite{Marsaglia1996}, 
and TestU01~\cite{Simard07testu01:a}, this latter being the
most stringent one. Finally, chaotic properties can be desired when
simulating a chaotic physical phenomenon or in hardware security, in 
which cryptographical proofs are not realizable.
In both truly and pseudorandom number generation, there is thus a need
to mathematically guarantee the presence of chaos, and to show that 
a post-treatment on a given secure and/or unbiased generator can be realized,
which adds chaos without deflating these desired properties.


This work takes place in this domain with the desire of automatically 
generating a large class of PRNGs with chaos and statistical properties.
In a sense, it is close to~\cite{BCGR11} where the authors shown that 
some Boolean maps may be embedded into an algorithm to provide a PRNG that has both
the theoretical Devaney's chaos property and the practical 
property of succeeding NIST statistical battery of tests.  
To achieve this, it has been proven in 
this article that it is sufficient   
for the iteration graph to be strongly connected,
and it is necessary and sufficient for its Markov probability matrix to be doubly stochastic.
However, they do not purpose conditions 
to provide such Boolean maps.
Admittedly, sufficient conditions
to retrieve Boolean maps whose graphs are 
strongly connected are given, but it remains to further filter those whose 
Markov matrix is doubly stochastic.
This approach suffers from delaying the second requirement to a final step
whereas this is a necessary condition. 
In this position article, we provide a completely new approach
to generate Boolean functions, whose Markov matrix is doubly stochastic and whose
graph of iterations is strongly connected. 
Furthermore the rate of convergence is always taken into consideration to provide 
PRNG with good statistical properties.


This research work is organized as follows.
It firstly recall some preliminaries that make the document self-contained (Section~\ref{sec:preliminaries}),
The next section (Section~\ref{sec:DSSC}) shows how this problem can be theoretically solved with a classical constraint logic programming.
Section~\ref{sec:hamiltonian} is the strongest contribution of this work.
It presents the main algorithm to generate Boolean maps with all the required properties and 
proves that such a construction is correct. 
Statistical evaluations are then presented in Section~\ref{sec:exp}. 
Conclusive remarks, open problems, and perspectives are 
finally provided.



 \section{\uppercase{Preliminaries}}\label{sec:preliminaries}

In what follows, we consider the Boolean algebra on the set 
$\Bool=\{0,1\}$ with the classical operators of conjunction '.', 
of disjunction '+', of negation '$\overline{~}$', and of 
disjunctive union $\oplus$. 

Let $n$ be a positive integer. A  {\emph{Boolean map} $f$ is 
a function from the Boolean domain 
 to itself 
such that 
$x=(x_1,\dots,x_n)$ maps to $f(x)=(f_1(x),\dots,f_n(x))$.
Functions are iterated as follows. 
At the $t^{th}$ iteration, only the $s_{t}-$th component is
``iterated'', where $s = \left(s_t\right)_{t \in \mathds{N}}$ is a sequence of indices taken in $\llbracket 1;n \rrbracket$ called ``strategy''. Formally,
let $F_f: \llbracket1;n\rrbracket\times \Bool^{n}$ to $\Bool^n$ be defined by
\[
F_f(i,x)=(x_1,\dots,x_{i-1},f_i(x),x_{i+1},\dots,x_n).
\]
Then, let $x^0\in\Bool^n$ be an initial configuration
and $s\in
\llbracket1;n\rrbracket^\Nats$ be a strategy, 
the dynamics are described by the recurrence
\begin{equation}\label{eq:asyn}
x^{t+1}=F_f(s_t,x^t).
\end{equation}


Let be given a Boolean map $f$. Its associated   
{\emph{iteration graph}}  $\Gamma(f)$ is the
directed graph such that  the set of vertices is
$\Bool^n$, and for all $x\in\Bool^n$ and $i\in \llbracket1;n\rrbracket$,
the graph $\Gamma(f)$ contains an arc from $x$ to $F_f(i,x)$. 

\begin{xpl}
Let us consider for instance $n=3$. Let 
$f^*: \Bool^3 \rightarrow \Bool^3$ be defined by
$f^*(x_1,x_2,x_3)=(x_2 \oplus x_3, x_1 \oplus \overline{x_3},\overline{x_3})$.
The iteration graph $\Gamma(f^*)$ of this function is given in 
Figure~\ref{fig:iteration:f*}.


\begin{figure}[ht]

\begin{center}
\includegraphics[scale=0.5]{iter_f0b.eps}
\end{center}
\caption{Iteration Graph $\Gamma(f^*)$ of the function $f^*$}\label{fig:iteration:f*}

\end{figure}
\end{xpl}

It is easy to associate a Markov Matrix $M$ to such a graph $G(f)$
as follows: 
\begin{itemize}
\item $M_{ij} = \frac{1}{n}$ if there is an edge from $i$ to $j$ in $\Gamma(f)$ and $i \neq j$;
\item $M_{ii} = 1 - \sum\limits_{j=1, j\neq i}^n M_{ij}$;
\item $M_{ij} = 0$ otherwise.
\end{itemize}

\begin{xpl}
The Markov matrix associated to the function $f^*$ is 
$$
\dfrac{1}{3}
\left(
\begin{array}{llllllll}
1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \\ 
1 & 1 & 0 & 0 & 0 & 1 & 0 & 0 \\ 
0 & 0 & 1 & 1 & 0 & 0 & 1 & 0 \\ 
0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 \\ 
1 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\ 
0 & 0 & 0 & 0 & 1 & 1 & 0 & 1 \\ 
0 & 0 & 0 & 0 & 1 & 0 & 1 & 1 \\ 
0 & 0 & 0 & 1 & 0 & 0 & 1 & 1 
\end{array}
\right)
$$
 

\end{xpl}


It is usual to check whether rows of such kind of matrices
converge to a specific 
distribution. 
Let us first recall the  \emph{Total Variation} distance $\tv{\pi-\mu}$,
which is defined for two distributions $\pi$ and $\mu$ on the same set 
$\Omega$  by:
$$\tv{\pi-\mu}=\max_{A\subset \Omega} |\pi(A)-\mu(A)|.$$ 
% It is known that
% $$\tv{\pi-\mu}=\frac{1}{2}\sum_{x\in\Omega}|\pi(x)-\mu(x)|.$$

Let then $M(x,\cdot)$ be the
distribution induced by the $x$-th row of $M$. If the Markov chain
induced by
$M$ has a stationary distribution $\pi$, then we define
$$d(t)=\max_{x\in\Omega}\tv{M^t(x,\cdot)-\pi}.$$
Intuitively $d(t)$ is the largest deviation between
the distribution $\pi$ and $M^t(x,\cdot)$, which 
is the result of iterating $t$ times the function.

Finally, let $\varepsilon$ be a positive number, the \emph{mixing time} 
with respect to $\varepsilon$ is given by
$$t_{\rm mix}(\varepsilon)=\min\{t \mid d(t)\leq \varepsilon\}.$$
It defines the smallest iteration number 
that is sufficient to obtain a deviation lesser than $\varepsilon$.
Notice that the upper and lower bounds of mixing times cannot    
directly be computed with eigenvalues formulae as expressed 
in~\cite[Chap. 12]{LevinPeresWilmer2006}. The authors of this latter work  
only consider reversible Markov matrices whereas we do no restrict our 
matrices to such a form.



Let us finally present the pseudorandom number generator $\chi_{\textit{14Secrypt}}$
which is based on random walks in $\Gamma(f)$. 
More precisely, let be given a Boolean map $f:\Bool^n \rightarrow \Bool^n$,
a PRNG \textit{Random},
an integer $b$ that corresponds to an awaited mixing time, and 
an initial configuration $x^0$. 
Starting from $x^0$, the algorithm repeats $b$ times 
a random choice of which edge to follow and traverses this edge.
The final configuration is thus outputted.
This PRNG is formalized in Algorithm~\ref{CI Algorithm}.




\begin{algorithm}[ht]
%\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$\;
\For{$i=0,\dots,b-1$}
{
$s\leftarrow{\textit{Random}(n)}$\;
$x\leftarrow{F_f(s,x)}$\;
}
return $x$\;
%\end{scriptsize}
\caption{Pseudo Code of the $\chi_{\textit{14Secrypt}}$ PRNG}
\label{CI Algorithm}
\end{algorithm}


This PRNG is a particularized version of Algorithm given in~\cite{BCGR11}.
Compared to this latter, the length of the random 
walk of our algorithm is always constant (and is equal to $b$) whereas it 
was given by a second PRNG in this latter.
However, all the theoretical results that are given in~\cite{BCGR11} remain
true since the proofs do not rely on this fact. 
Let us recall the following theorem.



\begin{theorem}[{\cite[Th. 4, p. 135]{BCGR11}}]
  Let $f: \Bool^{n} \rightarrow \Bool^{n}$, $\Gamma(f)$ its
  iteration graph, and $M$ its Markov matrix .
  If $\Gamma(f)$ is strongly connected, then 
  the output of the PRNG follows 
  a law that tends to the uniform distribution 
  if and only if $M$ is a doubly stochastic matrix.
\end{theorem}
  

With all this material, 
we can present  an efficient method to
generate Boolean functions
with Doubly Stochastic matrix and Strongly Connected iteration graph,
further (abusively) denoted as DSSC matrix.   



\section{\uppercase{Generation of DSSC Matrices}}\label{sec:DSSC}

This aim of this section is to show 
that finding DSSC matrices from a hypercube
is a typical finite domain satisfaction 
problem, classically denoted as 
Constraint Logic Programming on Finite Domains (CLPFD). 
This part is addressed in the first section. Next, we analyse the first
results to provide a generation of DSSC matrices with small mixing times. 

\subsection{Constraint Logic Programming on Finite Domains}

First of all,  
let $n$ be the number of elements. In order to avoid fractions 
in this article, we consider here that the 
sum of each column and each row is $n$. It  can easily be normalized to 1.  
The goal is thus to find all the $2^n\times 2^n$ matrices $M$ 
such that: 
  
\begin{enumerate}
\item $M_{ij}$ is 0 if $j$ is not a neighbor of $i$, \textit{i.e.},
  there is no edge from $i$ to $j$ in the $n$-cube.
\item $0 \le M_{ii} \le n$: the number of loops around $i$ is lesser than $n$  
\item Otherwise $0 \le M_{ij} \le 1$: if the edge from $i$ to $j$ is kept, 
$M_{ij}$ is 1, and 0 otherwise.
\item For any index of line $i$, $1 \le i\le 2^n$, $n = \sum_{1 \le j\le 2^n} M_{ij}$: 
  the matrix is right stochastic.
\item For any index of column $j$, 
  $1 \le j\le 2^n$, $n = \sum_{1 \le i\le 2^n} M_{ij}$: 
  the matrix is left stochastic. 
\item All the values of $\sum_{1\le k\le 2^n}M^k$ are strictly positive, \textit{i.e.}, the induced graph is strongly connected.
\end{enumerate}


Since these variables range into finite integer domains with sum and 
product operations, this problem can be theoretically handled by  
Constraint Logic Programming on Finite Domains (CLPFD), as implemented
in prolog.
The algorithm given in Figure~\ref{fig:prolog}
is indeed the core of the prolog file that is used to instantiate all the solutions
when $n$ is 2. In this code, 
\verb+mmult(+$X,Y,R$\verb+)+ and 
\verb+summ(+$X,Y,R$\verb+)+ 
 are true if and only if $R$ is the matrix product 
 or the matrix sum between 
 $X$ and $Y$ respectively. It is not hard to adapt such a code to any 
value of positive integer $n$.  

\begin{figure}[ht]
\begin{scriptsize}
\lstset{language=prolog}
\begin{lstlisting}
bistoc(X):-
  M=[[M0_0, M0_1, 0, M0_3], [M1_0, M1_1, 0, M1_3], 
     [M2_0, 0, M2_2, M2_3], [0, M3_1, M3_2, M3_3]],
  [M0_0, M1_1, M2_2, M3_3] ins 0..2,
  [M0_1, M0_3, M1_0, M1_3, M2_0, M2_3, M3_1, M3_2] 
     ins 0..1,
  M0_0+ M0_1+ M0_2 #=2, M1_0+ M1_1+ M1_3 #=2,
  M2_0+ M2_2+ M2_3 #=2, M3_1+ M3_2+ M3_3 #=2,
  M0_0+ M1_0+ M2_0 #=2, M0_1+ M1_1+ M3_1 #=2, 
  M0_2+ M2_2+ M3_2 #=2, M1_3+ M2_3+ M3_3 #=2,
  mmult(M,M,M2),
  mmult(M,M2,M3),
  mmult(M,M3,M4),
  summ(M,M2,S2),
  summ(S2,M3,S3),
  summ(S3,M4,S4),
  allpositive(S4).
\end{lstlisting}
\end{scriptsize}
\caption{Prolog Problem to Find DSSC Matrix when $n=2$}\label{fig:prolog}
\end{figure}

Finally, we define the relation $\mathcal{R}$, which is established on the 
two functions
$f$ and $g$ if the iteration graphs $\Gamma(f)$ and $\Gamma(g)$ 
are isomorphic. 
Obviously, this is an equivalence relation. 



\subsection{Analysis of the  Approach}

When executed on a personal computer, prolog finds 
in less than 1 second the 49 solutions when $n$ is 2, where only 2 are not equivalent, and in less than 1 minute the 27642 solutions where only 111 are
not equivalent when $n$ is 3. But it does not achieve the generation of 
all the solutions when $n$ is 4.
This approach suffers from not being efficient enough for large $n$ due to 
a \emph{generate and test} approach, despite the efficiency of the native 
backtrack of in CLP.

However, first results for small values of $n$ have been evaluated.  
More precisely, non equivalent generated
functions have been compared according to their 
ability to efficiently produce uniformly distributed outputs, \textit{i.e.},
to have the smallest mixing time.

\begin{xpl}
Table~\ref{table:mixing:3} gives the 5 best Boolean functions
ordered by their mixing times (MT) in the third column for $\varepsilon=10^{-5}$. 
\begin{table}[ht]
\begin{small}
$$
\begin{array}{|l|l|l|}
\hline
\textrm{Name} & \textrm{Image} & \textrm{MT}\\
\hline
f^* &  (x_2 \oplus x_3, x_1 \oplus \overline{x_3},\overline{x_3})  &  16   \\
\hline
f^a  &  (x_2 \oplus x_3, \overline{x_1}\overline{x_3} + x_1\overline{x_2},
\overline{x_1}\overline{x_3} + x_1x_2)  &  17   \\
\hline
f^b  &  (\overline{x_1}(x_2+x_3) + x_2x_3,\overline{x_1}(\overline{x_2}+\overline{x_3}) + \overline{x_2}\overline{x_3}, &  \\
& \qquad \overline{x_3}(\overline{x_1}+x_2) + \overline{x_1}x_2)  &  26   \\
\hline
f^c  &  (\overline{x_1}(x_2+x_3) + x_2x_3,\overline{x_1}(\overline{x_2}+\overline{x_3}) + \overline{x_2}\overline{x_3}, & \\
& \overline{x_3}(\overline{x_1}+x_2) + \overline{x_1}x_2)  &  29   \\
\hline
f^d  &  (x_1\oplus x_2,x_3(\overline{x_1}+\overline{x_2}),\overline{x_3})  &  30   \\
\hline
% [3, 4, 7, 0, 1, 6, 5, 2], #16
% [3, 4, 7, 0, 2, 6, 5, 1], #17
% [3, 6, 7, 5, 2, 0, 1, 4], #26 
% [3, 4, 7, 5, 2, 6, 1, 0], #29
% [3, 2, 5, 6, 7, 0, 1, 4]] #30
\end{array}
$$
\end{small}
\caption{The 5 Boolean functions with smallest MT when $n=3$.}\label{table:mixing:3}
\end{table}
\end{xpl}

A syntactical analysis of the functions for $n$=3 
does not help to understand
how to build a Boolean map with a small mixing time. 
However the iteration graph of $f^*$ 
(given in Figure~\ref{fig:iteration:f*}) 
is the $3$-cube in which the cycle 
$000,100,101,001,011,111,110,010,000$ 
has been removed.
This cycle that visits each vertex exactly once is usually denotes as 
\emph{Hamiltonian cycle}.
We are now focusing on the generation of DSSC matrices by  
removing Hamiltonian cycles of the $n$-cube. This is the aims of the next 
section.


\section{\uppercase{Removing Hamiltonian Cycles}}\label{sec:hamiltonian}
This section addresses the problem of removing an Hamiltonian path in the $n$-cube.
The first theoretical section shows that this approach produces DSSC matrix, as wished.
The motivation to focus on balanced Gray code is then given in Sec.~\ref{sub:gray}.
We end this section by recalling an algorithm that aims at computing such codes (Section~\ref{sub:algo:gray}).

\subsection{Theoretical Aspects of Removing Hamiltonian Cycles}

We first have the following result on stochastic matrix and $n$-cube without
Hamiltonian cycle.

\begin{theorem}
The Markov Matrix $M$ resulting from the $n$-cube in
which an Hamiltonian 
cycle is removed, is doubly stochastic.
\end{theorem}
\begin{proof}
An Hamiltonian cycle visits each vertex exactly once. 
For each vertex $v$ in the $n$-cube,
one ongoing edge $(o,v)$ and one outgoing edge $(v,e)$ 
are thus removed.

Let us consider the Markov matrix $M$ of the $n$-cube. 
It is obviously doubly stochastic.
Since we exactly remove one outgoing edge, the value of $M_{ve}$ decreases
from $\frac{1}{n}$  to 0 and
$M_{vv}$ is $\frac{1}{n}$. The $M$ matrix is stochastic again.
Similarly for ongoing edge, since one ongoing edge is dropped for each 
vertex, the value of $M_{ov}$ decreases
from $\frac{1}{n}$  to $0$. Moreover, since $M_{vv}$ is $\frac{1}{n}$,
the sum of values in column $v$ is $1$, and $M$ is doubly stochastic. 
\end{proof}

The following result states that the $n$-cube without one
Hamiltonian cycle 
has the awaited property with regard to the connectivity.

\begin{theorem}
The iteration graph issued from the $n$-cube where an Hamiltonian 
cycle is removed is strongly connected.
\end{theorem}
\begin{proof}
We consider the reverse cycle $r$ of the Hamiltonian cycle $c$.
There is no edge that belongs to both $r$ and $c$: otherwise $c$
would contain one vertex twice. Thus, no edges of $r$ has been removed.
The cycle $r$ is obviously an Hamiltonian cycle and contains all the nodes.
Any node of the $n$-cube where $c$ has been removed  can thus reach any node.  
The  iteration graph is thus strongly connected.  
\end{proof}

Removing an Hamiltonian cycle in the $n$-cube solves thus the DSSC constraint.
We are then left to focus on the generation of Hamiltonian cycles in the 
$n$-cube. Such a problem is equivalent to find Gray codes, \textit{i.e.},
to find a sequence of $2^n$ codewords ($n$-bits strings) 
where two successive elements differ in only one bit position. 
The next section is dedicated to these codes.

\subsection{Linking to (Totally) Balanced Gray Codes}\label{sub:gray}
Many research works~\cite{journals/combinatorics/BhatS96,VanSup04,journals/combinatorics/FlahiveB07} 
have addressed the subject of finding Gray  codes.
Since our approach is based on removing a cycle, we are concerned 
with cyclic Gray codes, \textit{i.e.}, sequences where the last codeword   
differs in only one bit position from the first one.

Let $n$ be a given integer. As far as we know, 
the exact number of Gray codes in $\Bool^n$ is not known but 
a lower bound, $$\left(\frac{n*\log2}{e \log \log n}*(1 - o(1))\right)^{2^n}$$
has been given in~\cite{Feder:2009:NTB:1496397.1496515}.
For example, when $n$ is $6$, such a number is larger than $10^{13}$.

To avoid this combinatorial explosion, we want to
restrict the generation to any Gray code 
such that the induced graph of iteration $\Gamma(f)$ is 
``uniform''. In other words, if we count in $\Gamma(f)$ 
the number of edges that modify the bit $i$, for $1 \le i \le n$,
all these values have to be close to each other.
Such an approach is equivalent to restrict the search of cyclic Gray codes
which are uniform too.

This notion can be formalized as follows. Let  $L = w_1, w_2, \dots, w_{2^n}$ be the sequence 
of a $n$-bits cyclic Gray code. Let $S = s_1, s_2, \dots, s_{2^n}$ be the 
transition sequence where $s_i$, $1 \le i \le 2^n$ indicates which bit position changes between 
codewords at index $i$ and $i+1$ modulo $2^n$. 
Let $\textit{TC}_n : \{1,\dots, n\} \rightarrow \{0, \ldots, 2^n\}$ the \emph{transition count} function
that counts the number of times $i$ occurs in $S$, \textit{i.e.}, the number of times 
the bit $i$ has been switched in $L$.   
The Gray code is \emph{totally balanced} if $\textit{TC}_n$ is constant (and equal to $\frac{2^n}{n}$).
It is \emph{balanced} if for any two bit indices $i$ and $j$, 
$|\textit{TC}_n(i) - \textit{TC}_n(j)| \le  2$.
   
\begin{xpl}
Let $L^*=000,100,101,001,011,111,110,010$ be the Gray code that corresponds to 
the Hamiltonian cycle that has been removed in $f^*$.
Its transition sequence is $S=3,1,3,2,3,1,3,2$ and its transition count function is 
$\textit{TC}_3(1)= \textit{TC}_3(2)=2$ and  $\textit{TC}_3(3)=4$. Such a Gray code is balanced. 

Let now  
$L^4=0000, 0010, 0110, 1110, 1111, 0111, 0011, 0001, 0101,$
$0100, 1100, 1101, 1001, 1011, 1010, 1000$
be a cyclic Gray code. Since $S=2,3,4,1,4,3,2,3,1,4,1,3,2,1,2,4$ $\textit{TC}_4$ is equal to 4 everywhere, this code
is thus totally balanced.

On the contrary, for the standard $4$-bits Gray code  
$L^{\textit{st}}=0000,0001,0011,0010,0110,0111,0101,0100,1100,$
$1101,1111,1110,1010,1011,1001,1000$,
we have $\textit{TC}_4(1)=8$ $\textit{TC}_4(2)=4$ $\textit{TC}_4(3)=\textit{TC}_4(4)=2$ and
the code is neither balanced nor totally balanced.
\end{xpl}

\subsection{Induction-Based Generation of  Balanced Gray Codes}\label{sub:algo:gray} 
The algorithm we adapt is based on the ``Construction B''~\cite{VanSup04} we recall now.
This method inductively constructs $n$-bits Gray code given a $n-2$-bit Gray code.

It starts with the transition sequence $S_{n-2}$ of such code.

\begin{enumerate}
\item \label{item:nondet}Let $l$ be an even positive integer. Find 
$u_1, u_2, \dots , u_{l-2}, v$ (maybe empty) subsequences of $S_{n-2}$ 
such that $S_{n-2}$ is the concatenation of 
$$
s_{i_1}, u_0, s_{i_2}, u_1, s_{i_3}, u_2, . . . , s_{i_l-1}, u_{l-2}, s_{i_l}, v
$$
where $i_1 = 1$, $i_2 = 2$, and $u_0 = \emptyset$ (the empty sequence).
\item Replace in $S_{n-2}$ the sequences $u_0, u_1, u_2, \ldots, u_{l-2}$ 
  by 
  $n - 1,  u'(u_1,n - 1, n) , u'(u_2,n, n - 1), u'(u_3,n - 1,n), \dots, u'(u_{l-2},n, n - 1)$
  respectively, where $u'(u,x,y)$ is the sequence $u,x,u^R,y,u$ such that 
  $u^R$ is $u$ in reversed order. The obtained sequence is further denoted as $U$.
\item Construct the sequences $V=v^R,n,v$, $W=n-1,S_{n-2},n$, and let $W'$ be $W$ where the first 
two elements have been exchanged.
\item The transition sequence $S_{n}$ is thus the concatenation $U^R, V, W'$.
\end{enumerate} 

It has been proven in~\cite{VanSup04} that 
$S_{n}$ is transition sequence of a cyclic $n$-bits Gray code 
if $S_{n-2}$ is. 
However, the step~\ref{item:nondet} is not a constructive 
step that precises how to select the subsequences which ensures that 
yielded Gray code is balanced.

Let us now  evaluate the number of subsequences $u$  
than can be produced. Since $s_{i_1}$ and $s_{i_2}$ are well
defined, we have to chose the $l-2$ elements of $s_3,s_4,\dots,s_{2^{n-2}}$
that become $s_{i_3},\dots, s_{i_l}$. Let $l = 2l'$.
There are thus 
$$
\#_n = \sum_{l'=1}^{2^{n-3}} {2^{n-2}-2 \choose 2l'-2}
$$
distinct subsequences $u$.
Numerical values of $\#_n$ are given in table~\ref{table:combinations}.
Even for small values of $n$, it is not reasonable to hope to evaluate the whole 
set of subsequences.
\begin{table}[ht]
\begin{center}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|}
\hline
$n$ &  5 & 6 & 7 & 8 & 9 \\  
\hline
$\#_n$ & 
31 & 8191 & 5.3e8 & 2.3e18 & 4.2e37\\
\hline
$\#'_n$ & 
15 & 3003 & 1.4e8 & 4.5e17 & 1.6e36 \\
\hline
\end{tabular}
\end{center}
\caption{Number of distinct $u$ subsequences.}\label{table:combinations}
\end{table}

However, it is  shown in the article that $\textit{TC}_n(n-1)$ and $\textit{TC}_n(n)$ are
equal to $l$. Since this step aims at generating (totally) balanced Gray codes,  
we have set $l$ to be the largest even integer less or equal than $\frac{2^n}{n}$.
This improvement allows to reduce the number of subsequences to study.
Examples of such cardinalities are given in table~\ref{table:combinations} and are refered as 
$\#'_n$.

Finally, the table~\ref{table:nbFunc} gives the number of non-equivalent functions issued from 
(totally) balanced Gray codes that can be generated
with the approach presented in this article with respect to the number of bits. In other words, it corresponds to the size of the class of generators that can be produced.

\begin{table}[ht]
\begin{center}
\begin{tabular}{|l|c|c|c|c|}
\hline
$n$              & 3 & 4 & 5 & 6  \\
\hline
nb. of functions & 1 & 1 & 2 & 1332   \\
\hline
\end{tabular}
\end{center}
\caption{Number of Generators w.r.t. the number of bits.}
\label{table:nbFunc}
\end{table}


\section{\uppercase{Experiments}}\label{sec:exp}
We present in Algorithm~\ref{CI} the method
that allows to take any chaotic function as the core of a pseudo
random number generator. Among the parameters, it takes the
number b of minimal iterations that have to be executed to
get a uniform like distribution. For function $f$ and our experiments $b$ is set
with the value given in the fourth column of Table~\ref{table:nc}.
\begin{algorithm}[ht]
%\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$\;
\For{$i=0,\dots,b-1$}
{
$s\leftarrow{\textit{Random}~mod~n}$\;
$x\leftarrow{(x-(x\&(1<<s))+f(x)\&(1<<s))}$\;
}
return $x$\;
%\end{scriptsize}
\caption{Pseudo Code of the $\chi_{\textit{14Secrypt}}$ PRNG}
\label{CI}
\end{algorithm}

%$$ 
%\begin{array}{l}
%(\overline{x_1}(\overline{x_2}.x_4)+\overline{x_2}.\overline{x_3}.\overline{x_4}+ x_2.x_3%.x_4,\\
%x_1.\overline{x_2} + x_1.\overline{x_3}.\overline{x_4}+\overline{x_1}x_3x_4 + \overline{x_%2}.\overline{x_3}.\overline{x_4},\\
% x_2.\overline{x_3}+\overline{x_2}(x_1 \oplus x_3),\\
%\overline{x_2}.\overline{x_3}+x_1.\overline{x_2}x_3+\overline{x_1}x_2(\overline{x_3}
%+\overline{x_4}))
%\end{array}
%$$


\begin{table*}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
Function $f$ & $f(x)$, for $x$ in $(0,1,2,\hdots,2^n-1)$ & $n$ & $b$ \\ 
\hline
$\textcircled{a}$&[13, 10, 9, 14, 3, 11, 1, 12, 15, 4, 7, 5, 2, 6, 0, 8]&4&32\\
% (
%\overline{x_1}(\overline{x_2}x_4)+\overline{x_2}\overline{x_3}\overline{x_4}+ x_2x_3x_4,
%x_1\overline{x_2} + x_1\overline{x_3}\overline{x_4}+\overline{x_1}x_3x_4 + \overline{x_2}\overline{x_3}\overline{x_4},
% x_2\overline{x_3}+\overline{x_2}(x_1 \oplus x_3),
%\overline{x_2}\overline{x_3}+x_1\overline{x_2}x_3+\overline{x_1}x_2(\overline{x_3}+\overline{x_4}))
\hline
&[55, 60, 45, 56, 58, 62, 61, 40, 53, 38, 52, 54, 35, 51, 33, 49, 39, 14, 
&&\\
$\textcircled{b}$&
47, 46, 59, 43, 57, 44, 37, 6, 36, 4, 3, 50, 1, 48, 63, 26, 25, 30, 19, 
&6&49\\
&
27, 17, 28, 31, 20, 23, 21, 18, 22, 16, 24, 13, 12, 29, 8, 10, 42, 41, 
&&\\
&
0, 15, 2, 7, 5, 11, 34, 9, 32]&&\\

\hline
&[223, 250, 249, 254, 187, 234, 241, 252, 183, 230, 229, 180, 227, 178, &&\\
&
240, 248, 237, 236, 253, 172, 251, 238, 201, 224, 247, 166, 165, 244, 
&&\\
&
163, 242, 161, 225, 215, 220, 205, 216, 218, 222, 221, 208, 213, 210, 
&&\\
&
135, 196, 199, 132, 194, 130, 129, 200, 159, 186, 185, 190, 59, 170, 
&&\\
&
177, 188, 191, 246, 245, 52, 243, 50, 176, 184, 173, 46, 189, 174, 235, 
&&\\
&
42, 233, 232, 231, 38, 37, 228, 35, 226, 33, 168, 151, 156, 141, 152, 
&&\\
&
154, 158, 157, 144, 149, 146, 148, 150, 155, 147, 153, 145, 175, 14, 
&&\\
$\textcircled{c}$&
143, 204, 11, 202, 169, 8, 7, 198, 197, 4, 195, 2, 1, 192, 255, 124, 
&8&75\\
&
109, 120, 107, 126, 125, 112, 103, 114, 116, 100, 123, 98, 121, 113, 79, 
&&\\
&
106, 111, 110, 75, 122, 97, 108, 71, 118, 117, 68, 115, 66, 96, 104, 
&&\\
&
127, 90, 89, 94, 83, 91, 81, 92, 95, 84, 87, 85, 82, 86, 80, 88, 77, 76, 
&&\\
&
93, 72, 74, 78, 105, 64, 69, 102, 101, 70, 99, 67, 73, 65, 55, 60, 45, 
&&\\
& 
56, 51, 62, 61, 48, 119, 182, 181, 53, 179, 54, 57, 49, 15, 44, 47, 40, 
&&\\
&
171, 58, 9, 32, 167, 6, 5, 164, 3, 162, 41, 160, 63, 26, 25, 30, 19, 27, 
&&\\
&
17, 28, 31, 20, 23, 21, 18, 22, 16, 24, 13, 10, 29, 140, 43, 138, 137, 
&&\\
&
12, 39, 134, 133, 36, 131, 34, 0, 128]&&\\
\hline

\end{tabular}

\end{center}
\caption{Functions with DSCC Matrix and smallest MT\label{table:nc}}
\end{table*}


For each number $n=4, 6, 8$ of bits, we have generated all the functions according the method 
given in Section~\ref{sec:hamiltonian}.
For each $n$, we have then restricted this evaluation to the function 
whose Markov Matrix has the smallest mixing time. Such functions are 
given in Table~\ref{table:nc}.
In this table, let us consider for instance 
the function $\textcircled{a}$ from $\Bool^4$ to $\Bool^4$
defined by the following images : 
$[13, 10, 9, 14, 3, 11, 1, 12, 15, 4, 7, 5, 2, 6, 0, 8]$.
In other words,  the image of $3~(0011)$ by $\textcircled{a}$ is $14~(1110)$: it is obtained as  the  binary  value  of  the  fourth element  in  the  second  list
(namely~14).  


\subsection{NIST}
\label{NIST}
In our experiments, 100 sequences (s = 100) of 1,000,000 bits are generated and tested.
If the value $\mathbb{P}_T$ of any test is smaller than 0.0001, the sequences are considered to be not good enough
and the generator is unsuitable. Table~\ref{The passing rate} shows $\mathbb{P}_T$ of sequences based on discrete
chaotic iterations using different schemes. If there are at least two statistical values in a test, this test is
marked with an asterisk and the average value is computed to characterize the statistics.
We can see in Table \ref{The passing rate} that all the rates are greater than 97/100, \textit{i. e.}, all the generators pass the NIST test.
\begin{table*} 
\renewcommand{\arraystretch}{1.3}
\caption{NIST SP 800-22 test results ($\mathbb{P}_T$)}
\label{The passing rate}
\begin{center}
  \begin{tabular}{|l||c||c||c|}
    \hline
Method &$\textcircled{a}$& $\textcircled{b}$ & $\textcircled{c}$   \\ \hline\hline
Frequency (Monobit)& 0.678 (0.99)& 0.574 (0.99)& 0.699 (0.96)\\ \hline 
Frequency  within a Block& 0.102 (0.99)& 0.816 (0.98)& 0.419 (0.99)\\ \hline 
Runs& 0.171 (0.98)& 0.657 (0.98)& 0.554 (0.99)\\ \hline 
Longest Run of Ones in a Block& 0.115 (0.98)& 0.534 (1.0)& 0.534 (0.98)\\ \hline 
Binary Matrix Rank& 0.401 (0.97)& 0.554 (0.99)& 0.911 (0.99)\\ \hline 
Discrete Fourier Transform (Spectral)& 0.554 (0.98)& 0.350 (0.98)& 0.080 (0.99)\\ \hline 
Non-overlapping Template Matching*& 0.509 (0.990)& 0.443 (0.990)& 0.499 (0.989)\\ \hline 
Overlapping Template Matching& 0.437 (0.99)& 0.699 (1.0)& 0.236 (0.99)\\ \hline 
Maurer's "Universal Statistical"& 0.171 (0.99)& 0.000 (0.97)& 0.657 (0.99)\\ \hline 
Linear Complexity& 0.171 (1.0)& 0.637 (0.99)& 0.834 (1.0)\\ \hline 
Serial* (m=10)& 0.435 (1.0)& 0.565 (0.98)& 0.592 (1.0)\\ \hline 
Approximate Entropy (m=10)& 0.137 (0.98)& 0.867 (0.99)& 0.062 (1.0)\\ \hline 
Cumulative Sums (Cusum) *& 0.580 (0.99)& 0.368 (0.99)& 0.569 (0.97)\\ \hline 
Random Excursions *& 0.245 (0.980)& 0.421 (0.991)& 0.656 (0.985)\\ \hline 
Random Excursions Variant *& 0.292 (0.991)& 0.450 (0.990)& 0.520 (0.996)\\ \hline
Success & 15/15 & 15/15 & 15/15 \\ 
\hline
\end{tabular}
  \end{center}
\end{table*}

\subsection{DieHARD}
\label{Subsec:DieHARD}

Table~\ref{Results of DieHARD battery of tests} gives the results derived from applying the DieHARD battery~\cite{Marsaglia1996} of tests
to the PRNGs considered in this work.
As it can be observed, all the generator presented in this document can pass the DieHARD battery of tests.

\begin{table}[htc]
\renewcommand{\arraystretch}{1.3}
\caption{Results of DieHARD battery of tests}
\label{Results of DieHARD battery of tests}
\begin{center}
\begin{tabular}{@{}l@{}l@{}ccc@{}} \toprule
\textbf{No.} &\textbf{Test name} &\multicolumn{3}{c}{\textbf{Generators}} \\ \cmidrule(r){3-5}
& & $\textcircled{a}$ &  $\textcircled{b}$ & $\textcircled{c}$ \\ \midrule
1 & Overlapping Sum &Pass &Pass &Pass \\
2 & Runs Up 1 &Pass & Pass &Pass \\
&Runs Down 1 &Pass &Pass &Pass \\
&Runs Up 2 &Pass &Pass &Pass \\
&Runs Down 2 &Pass & Pass &Pass \\
3 & 3D Spheres &Pass &Pass &Pass \\
4 & Parking Lot &Pass &Pass &Pass \\
5 & Birthday Spacing &Pass &Pass &Pass \\
6 & Count the ones 1 &Pass &Pass &Pass \\
7 &Binary Rank $6 \times 8$ &Pass & Pass &Pass \\
8 &Binary Rank $31 \times 31$ &Pass &Pass &Pass \\
9 &Binary Rank $32 \times 32$ &Pass &Pass &Pass \\
10 &Count the ones 2 &Pass &Pass&Pass \\
11 &Bit Stream &Pass &Pass&Pass  \\
12 &Craps Wins &Pass &Pass&Pass  \\
&Throws &Pass &Pass &Pass  \\
13 &Minimum Distance &Pass &Pass &Pass  \\
14 &Overlapping Perm. &Pass &Pass &Pass  \\
15 &Squeeze &Pass &Pass&Pass \\
16 &OPSO &Pass &Pass&Pass \\
17 &OQSO &Pass &Pass&Pass  \\
18 &DNA &Pass &Pass&Pass  \\
&Number of tests passed &18 &18 &18 \\\bottomrule
\end{tabular}
\end{center}
\end{table}


\subsection{TestU01}
\label{Subsec:TestU01}

TestU01~\cite{Simard07testu01:a} is a software library that 
provides general implementations of the classical statistical tests for 
random number generators, as well as several others proposed in the 
literature, and some original ones. This library is currently the
most reputed and stringent one for testing the randomness
profile of a given sequence. TestU01 encompasses the NIST and DieHARD
tests suites with 2 batteries specific to hardware based generators 
(namely, Rabbit and Alphabit). Its three core batteries of tests are
however SmallCrush, Crush, and BigCrush, classified according to
their difficulty.

To date, we can claim after experiments that $\textcircled{a}$
generator is able to pass the 15 tests embedded into the 
SmallCrush battery and it succedded too to pass the 144 tests of 
the Crush one. BigCrush results on $\textcircled{a}$ are
expected soon, while TestU01 has been launched too on generators
having the other iteration functions detailed in this article.


\section{\uppercase{Conclusion}}\label{sec:conclusion}
This article has presented a method to compute a large class of truly chaotic PRNGs.
First experiments through the batteries of NIST, DieHard, and TestU01
have shown that the statistical properties are almost established for $n=4,6,8$.
The iterated map inside the generator is built by removing 
from a $n$-cube an Hamiltonian path that corresponds to a (totally) balanced Gray code.  
The number of balanced gray code is large and each of them can be 
considered as a key of the PRNG.
However, many problems still remain open, most important ones being listed 
thereafter. 

The first one involves the function to iterate. Producing a DSSC matrix is indeed necessary and sufficient
but is not linked with the convergence rate to the uniform distribution. 
To solve this problem, we have proposed to remove from the $n$-cube an Hamiltonian path that 
is a (totally) balanced Gray code. We do not have proven that this 
proposal is the one that minimizes the  mixing time. This optimization task is an open problem we plan to study.

Secondly, the approach depends on finding (totally) balanced Gray 
codes. Even if such codes exist for all even numbers, there is no constructive method to built them
when $n$ is large, as far as we know. These two open problems will
be investigated in a future work.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{apalike}
{\small
\bibliography{markov}}
\end{document}