Jonathan Passerat-Palmbach

Claude Mazel

Bruno Bachelet

David R C Hill

ShoveRand: a Model-Driven Framework to Easily Generate Random Numbers on GP-GPU

Keywords: General-Purpose computation on Graphics Processing Units (GPGPU), Libraries and Programming Environments, Parallelization of Simulation, Object Oriented Programming & Design

Stochastic simulations are often sensitive to the randomness source that characterizes the statistical quality of their results. Consequently, we need highly reliable Random Number Generators (RNGs) to feed such applications. Recent developments try to shrink the computation time by using more and more General Purpose Graphics Processing Units (GP-GPUs) to speed-up stochastic simulations. Such devices bring new parallelization possibilities, but they also introduce new programming difficulties. Since RNGs are at the base of any stochastic simulation, they also need to be ported to GP-GPU. There is still a lack of well-designed implementations of qualityproven RNGs on GP-GPU platforms. In this paper, we introduce ShoveRand, a framework defining common rules to generate random numbers uniformly on GP-GPU. Our framework is designed to cope with any GPU-enabled development platform and to expose a straightforward interface to users. We also provide an existing RNG implementation with this framework to demonstrate its efficiency in both development and ease of use.

INTRODUCTION

Random number generation is a sensitive component for stochastic simulation models. For computationally eager applications, High Performance Computing (HPC) devices such as General-Purpose Graphics Processing Units (GP-GPU) are now widely used to reduce computation time. However, we still need better software tools to develop applications on such architectures. CPUenabled simulations can take advantage of a wide range of statistically sound Pseudo-Random Number Generators (PRNGs) and libraries, but we still need quality random number generators for GP-GPU architectures, and more precisely, a particular care has to be given to their parallelization.

Many tries have been done to propose GPU-enabled implementations of PRNGs. Previous studies like [START_REF] Passerat-Palmbach | Reliable initialization of gpuenabled parallel stochastic simulations using mersenne twister for graphics processors[END_REF][START_REF] Passerat-Palmbach | Pseudo-random number generation on gpgpu[END_REF] bring up random number generation on GPU. Many strategies were adopted, but few of them distinguish particularly from the others. It is also very difficult to select the right PRNG to use since you need to consider both the statistical quality of the PRNG and the GPU implementation quality if you do not want to see the overall performance of your application seriously drop.

Recent development frameworks like CUDA and OpenCL enabled much more ambitious developments on GPU. We can now make use of some software engineering paradigms, such as Object-Oriented programming or generic programming, and this introduces a new scope in GPU-enabled applications development. As GPU-enabled software become more complex, they also need enhanced tools to rely on. Thus, PRNGs implementations on GP-GPU cannot be achieved in non-structured software development environments anymore. We need a framework to smoothly integrate PRNGs in existing developments.

In this article, we propose a solution merging two main ideas. First, we want to enforce the use of good quality PRNGs on GPU. To do so, this work intends to aggregate well-designed random number generation solutions in order to distribute them in a packaged form. Second, we propose reflections on the form this package will take. To secure this aspect, we rely on sound software engineering principles to provide a modern framework. The latter must display a common interface to its user, regardless of the RNG he has chosen. Another important aspect, in our mind, is that this framework must evolve easily thanks to external developments. It is designed in such a way, to quickly integrate any already developed GPU-enabled RNG.

In this study, we will:

• Propose a meta-model for RNG libraries on GPU;

• Define implementations constraints to build such a dedicated framework;

• Implement a generic declination of a quality-proven GPU-enabled PRNG.

GPU-ENABLED RNG LIBRARIES

This section will present two libraries dedicated to random numbers generation on GP-GPU architectures: NVIDIA's CURAND [13] and Thrust::random [START_REF] Hoberock | Thrust: A parallel template library[END_REF]. They are, to our knowledge, the most relevant libraries to obtain random sequences on GP-GPU, at the time of writing.

Next, we will describe their intrinsic properties, considering both the RNGs they embed and the quality of their implementations.

NVIDIA's CURAND

Presentation

Introduced in the 3.2 version of the CUDA framework, CURAND has been designed to generate random numbers in a straightforward way on CUDA-enabled GP-GPUs. The main advantage of CURAND is that it does not target any platform or application kind. It is indeed able to produce both quasi-random and pseudo-random sequences either on GP-GPU or on CPU. Except from the initialization part, the API (Application Programming Interface) of the library stays globally the same, no matter which kind of random sequence you select or the platform you run your application on.

Characteristics

The current version of CURAND, shipped with CUDA SDK 3.2, comes with two kinds of random numbers (pseudo-random and quasi-random). This leads to two different implementations cohabitating in the library. On the one hand, we find a Sobol quasi-random numbers generator (QRNG), a very good choice that has also been done for the GNU Scientific Library (GSL).

The "quasi" modifier denotes that we do not deal with random or pseudo-random numbers. Instead, we are in presence of low-discrepancy deterministic sequences, which are interesting to improve numerical integration or other applications with better space filling and density properties. Toughly speaking, these generators seek to uniformly fill an n-dimensional space with points, without creating clusters of points. They are mainly used with quasi-Monte-Carlo computations (in physics, finance, etc.). On the other hand, the pseudo-random side, a XORWOW [START_REF] Marsaglia | Xorshift rngs[END_REF] PRNG has been implemented. Direct applications of PRNGs are stochastic simulations, which require fast randomness sources but also need to be reproducible as any scientific experiment and also for debugging purposes. Since we focus on simulations fed by pseudo-random sequences, we will only discuss the XORWOW algorithm. It is reckoned for its speed and small memory footprint. The latter property is a significant advantage when you need to implement this PRNG on a GP-GPU. The quicker the memory areas are, the fewer amount of data they are able to store. Thus, a PRNG saving memory is more likely supposed to be implemented using quick memories on GP-GPU. It is also fairly unpleasant to waste large amounts of bytes to store data for a PRNG when you still have to place your algorithm's data in the quickest memory areas you have at your disposal.

However, the main drawback we found with this algorithm is that it has shown limited statistical qualities in the literature [START_REF] Panneton | Construction d'ensembles de points basée sur des récurrences linéaires dans un corps fini de caractéristique 2 pour la simulation monte carlo et l'intégration quasi-monte carlo[END_REF]. The CURAND documentation even specifies that problems could be encountered when using this PRNG [13] (pp. 38-39).

Thrust::random

Presentation

Thrust::random is part of a GP-GPU-enabled general purpose library called Thrust. This opensource project intends to provide a GP-GPU-enabled library equivalent to some classic generalpurpose C++ libraries, such as STL or Boost. Classes are split through several namespaces, which Thrust::random is an example. The latter contains all classes and methods related to random numbers generation on GP-GPU. Thrust::random implements three PRNGs, each through a different C++ template class. We find a Linear Congruential Generator (LCG), a Linear Feedback Shift (LFS) and a Subtract With Carry (SWC), which we will all describe later and that are, in our opinion, not adapted to High Performance Computing. These classes are the base of this module that reuses them to define combined PRNGs or API adaptors of these original classes.

Characteristics

We need to distinguish two aspects when attempting to characterize this library. First, the template form of Thrust::random enables it to execute quickly on the host. Generic programming is used here to bypass the fact that we do not have true object-oriented features, such as inheritance and polymorphism, which would however imply an overhead at runtime. Polymorphism is in fact not available on every GP-GPU architecture yet, as we will see in Section 4.1. Even if the three PRNGs implementations are not directly related through a class hierarchy, they display an unchanged API. This way, PRNG intrinsic particularities are hidden to users who can simply call identical methods through every object issued by these classes. Furthermore, the library makes a clever use of the Adapter design pattern [START_REF] Gamma | Design patterns: elements of reusable object-oriented software[END_REF] to furnish combined PRNGs or to return extracts of an original random sequence. Nonetheless, Thrust::random does not provide a unified framework to add new PRNGs so that a potential newcomer would have to respect implicitly library conventions in order to integrate his own PRNG implementation in it.

Second, the substance of Thrust::random, formed by the underlying PRNGs, may also be discussed. We previously mentioned three PRNGs implementations. Let us now focus on their statistical qualities according to the literature. The first to be considered is a LCG, but we will only present it briefly, because it has several times been spotted as a flawed PRNG for both sequential [START_REF] Ecuyer | Testing random number generators[END_REF][START_REF] Ecuyer | Testu01: A c library for empirical testing of random number generators[END_REF] and parallel [START_REF] Hill | Practical distribution of random streams for stochastic high performance computing[END_REF][START_REF]Encyclopedia of Quantitative Finance[END_REF] applications. Therefore, we discourage users to choose it, even through Thrust's implementation. The LFS PRNG also implemented in Thrust::random is derived from the Boost.Random library. The latter implements in fact a Tausworthe PRNG introduced in [START_REF] Tausworthe | Random numbers generated by linear recurrence modulo two[END_REF]. Tausworthe generator alone is known to have bad empirical performance, as stated in [START_REF] Hellekalek | Good random number generators are (not so) easy to find[END_REF]. That is why it is commonly coupled with other generators to form combined generators. Hopefully, generators can be easily combined thanks to Thrust::random design. The third PRNG proposed is a Subtract With Carry PRNG announced to be based upon [START_REF] Marsaglia | A random number generator for pc's[END_REF][START_REF] Marsaglia | A new class of random number generators[END_REF]. Please note that Marsaglia rather proposed a Subtract With Borrow (SWB). . . Nevertheless, this generator does not seem flawless. It misses some classics statistical tests such as BirthdaySpacing according to [START_REF] Ecuyer | Testing random number generators[END_REF]. A more recent study from the same authors [START_REF] Ecuyer | Testu01: A c library for empirical testing of random number generators[END_REF] even likens SWB to LCG, which are structurally weak for modern applications as we previously stated.

A MODEL-DRIVEN LIBRARY TO OVERCOME KNOWN ISSUES

In the light of the previous descriptions, we spotted lacks in different domains involved in RNG libraries design: on the one hand, the embedded RNGs choice and on the other hand, the API design. These shortages showed the importance of defining a novel framework for RNG implementations on GP-GPU. To enable developers to use RNGs easily, we still trust the library approach is the best option. In this section we will introduce the choices we made to prevent our model-driven framework to display such misses. Therefore, we consider two axes in this presentation: the PRNG that will be embedded in ShoveRand and the meta-model which our library relies on.

Default shipped PRNG

Libraries are a reliable way to spread software elements. Considering our previous assertions on PRNGs embedded in other libraries, we took care to issue quality proven generators with our proposal. In this work, we have selected Mersenne Twister for Graphics Processors (MTGP) [START_REF] Saito | A variant of mersenne twister suitable for graphics processors[END_REF], benchmarked in [START_REF] Passerat-Palmbach | Reliable initialization of gpuenabled parallel stochastic simulations using mersenne twister for graphics processors[END_REF] and recently improved by Matsumoto and Saito. Briefly, this PRNG inherits from former Mersenne Twister generators. It is based upon the same kind of algorithm that allows its elders to display huge periods. MTGP can bear periods up to 2 110503 . This feature is not mandatory for modern applications as said in [START_REF] Hill | Practical distribution of random streams for stochastic high performance computing[END_REF][START_REF]Encyclopedia of Quantitative Finance[END_REF] and can be a drawback when we consider the memory footprint of the generator. As a matter of fact, the higher the period is, the more memory is consumed.

Moreover, as for the initial Mersenne Twister generator, MTGP is shipped with an algorithm called Dynamic Creator (DC) that issues independent parameter sets, called parameterized statuses hereafter. DC enables parallelization through the parameterization technique, which is supposed to furnish independent random sequences thanks to an upstream step providing a distinct parameterized status to each parallel element. A unique identifier is directly integrated as a part of the characteristic polynomial of the matrix that defines the recurrence, and belongs to the parameterized status of MTGP. Two identifiers will consequently lead to two different parameterized statuses. Furthermore, DC ensures that the characteristic polynomials we get are mutually prime. Its authors assert that the random sequences generated with such distinct parameterized statuses will be highly independent. Even if this fact cannot be mathematically proven, it is widely admitted in the scientific community.

The couple formed by MTGP and DC fulfills the requirements of a parallel PRNG, thanks to its ability to produce independent random sequences. Details of the integration process of this generator in ShoveRand are provided in Section 4.3.

Inputs of the meta-model

To improve the software design, we propose hereafter a meta-model to represent RNGs running on any platform. As we have seen with the previous examples, there is currently no RNG framework on GP-GPU. This leads to heterogeneous GP-GPU-enabled RNG propositions. The latter can even sometimes be of very good quality but a bit tough to use for non-RNG experts. Thus, our goal is to design this generic RNG framework and furnish an implementation embedding quality-proven GP-GPU-enabled RNGs. If we tried to avoid the previously noticed flaws in the literature propositions, we also acknowledge good design elements that we reused without any hesitation in our model.

Once again, ShoveRand is designed as a base framework for other developers to integrate their GP-GPU-enabled RNGs implementations. We intend to furnish an empty shell that will require a small amount of work to embed a new RNG. By doing this, we detach our production from potential deprecation of the RNGs it proposes. Newer implementations will always be available at the smallest effort of integration, in the direct line of our primarily announced goal: simplicity of use.

The previous conclusions led us to work on a meta-model matching the evocated needs. Hopefully, most of the RNGs share common characteristics that we were able to factorize in few classes. Before describing the machinery of our library model, let us introduce its schematized version through the UML class diagram in Figure 1:

We have chosen to handle each part of a classical RNG through different classes. This distinction enables our model to represent various kinds of RNGs, from the most naive to the most complicated one, regardless of the type of the generated random numbers. We obtain a resulting hierarchy constituted by four base classes for the library: RNG, Algorithm, ParameterizedStatus For understanding purpose, parameters representation needs to be detailed. As it is described on the class diagram, we split them in two classes. First, ParameterizedStatus contains the effective input parameters that determine the whole random sequence. It is highly tied to the Algorithm using it, so that a couple of instances from these two classes define a random sequence. Second, SeedStatus mirrors the internal state vector that stores the current state of the generator. Please note that the concrete implementations of these two Status classes are bound to the Algorithm employing them. You can notice that this relation is explicitly drawn for PRNG2 only on the class diagram in Figure 1, in order not to overload the diagram. Still, every RNG is intended to do so. It is also compulsory for each Algorithm implementation to furnish a representation of each of the two introduced classes, even if one is empty. For instance, some PRNGs do not employ parameters, so such generators will furnish an empty implementation of the ParameterizedStatus class.

Finally, every member of this meta-model is generic. This feature allows us to abstract the algorithms from the data type (noted T on the figures) of the generated random numbers, which are commonly 32-bit or 64-bit integers, or floating-point numbers of single or double precision. The library is also abstracted from the underlying programming framework used to develop applications that will run on GP-GPU. Indeed, the core of the library is not bound to a particular GP-GPU device code technology, allowing developers to write RNGs implementations in their favorite language (CUDA or OpenCL), as long as it uses a C++ compliant compiler. For instance, the PRNG shipped with our library is a CUDA implementation perfectly compiling with NVIDIA's nvcc CUDA compiler.

META-MODEL IMPLEMENTATION

We intend to offer a straightforward API, enabling one to call a RNG without wondering indefinitely which parameters he should set to use it correctly. As it has been done in the libraries we studied, setting RNGs parameters to default values prevents users to set parameters to wrong values, without understanding their meaning.

Other libraries already proposed a common API regardless of the employed RNG. This is in our mind the key feature of these libraries. Indeed, it can be interesting to rapidly switch between two or more RNGs in order to compare the output results of the same application depending on its randomness source. With a uniform API, enhanced by generic programming, the only thing you have to do is to change a single parameter. Consequently, we obviously reused this design element in our meta-model. This section will show how ShoveRand takes advantage of the introduction of generic programming to propose even more advanced features.

Policy-based class design

Several pseudo-random number generation algorithms can be implemented in ShoveRand through different classes. Indeed, the RNG class gets its behavior from these classes. This situation largely recalls the Strategy design pattern [START_REF] Gamma | Design patterns: elements of reusable object-oriented software[END_REF], which intends to dynamically modify the behavior of the class instances according to the Strategy they are combined with. This feature simply relies on virtual functions calls. Unfortunately, even cutting-edge GP-GPUs prevent us from using strong polymorphism, i.e. virtual functions. As a matter of fact, GP-GPU architectures do not support the implementation of virtual functions tables, at the time of writing. Yet, this latter capacity is mandatory to enable polymorphism. Hopefully, [START_REF] Alexandrescu | Modern C++ Design[END_REF] introduced a compile time equivalent to the Strategy design pattern, which are called policies.

A policy defines a software component designed to be a unit of behavior. These components are basic classes that can potentially be combined to form complex classes. Policies intend to be an efficient design element, in this way they are passed as template parameters to the classes using them: the host classes. Now, remember that the dynamic Strategy solution presented the aspect to display a uniform API, thanks to inheritance and virtual functions. Policies provide an equivalent to this notion, given that they are a set of rules defining how a class should look like in order to furnish the features they propose. These rules constrain the interface of a policy's implementations classes, called policy classes, to guide their behavior towards the policy's expectations. Finally, all the policy classes implementing the same policy must display an identical interface, so that they can be swapped without any problem in the class using them.

In our case, the policy idiom is applied to the random number generation algorithm. Concretely, Algorithm appears as a template parameter of the RNG class. Each policy class implements a different kind of generator that can be interfaced with RNG and used identically. Thus, policies helped us to overcome hardware limitations and to set up an equivalent to the polymorphic Strategy design pattern. The Algorithm policy can be formulated as follows: Algorithm prescribes a class template of one type T representing the type of the generated values. Algorithm exposes at least two member-functions init and next that respectively initializes the generator and generates the next random value. The implementation must own representations of both a SeedStatus and a ParameterizedStatus, in order to handle its internal parameters.

Strong static typing

Template meta-programming presents other advantages in our case than genericity facilities. To structure our GP-GPU-enabled RNG framework, it is important to define some rules that will guide the future RNG implementations. The previous model fulfills this need by two aspects.

First, when RNG is defined for a particular (T; Algorithm) couple, it defines a new type allowing us to guess statically, understand at compilation time, whether two RNG instances possess the same template parameters. This feature can be very useful when using RNGs in parallel since we need to ensure the independence between random streams produced by the RNGs. For statistical analysis, it would be a nonsense to mix random sources during successive executions of a simulation. Such a misconception could wreak havoc on a stochastic simulation, for example by introducing undesired correlations in the random streams allocated to its replications. When CURAND exposes functions that take any of their implemented generators' status as parameter, our strong typing feature, coupled with an object-oriented design, prevents users to feed their applications with different RNGs for the same execution.

Second, we set up the concept checking mechanism at the heart of the library, in order to prevent users to provide any class as an Algorithm to RNG. Concept checking is a generic programming feature available through different implementations with C++. Introduced by [START_REF] Siek | Concept checking: Binding parametric polymorphism in c++[END_REF] and implemented in the Boost Concept Check Library (BCCL), concept checking enables us to statically check the interface of a given class. This mechanism enables us to ensure of the correctness of the interfaces of the policy classes. In the previous part, we have enounced a constraints list to define the Algorithm policy. We have translated this list into concept checking rules thanks to BCCL. For instance, the Algorithm policy can be basically checked with the few lines exposed in Figure 2.

// c o n c e p t r e q u i r e m e n t s BOOST_CONCEPT_USAGE( RNGAlgorithm ) { // r e q u i r e Algo<T> : : i n i t ( ) al_ . i n i t ( ) ; // r e q u i r e T Algo<T> : : next ( ) value_ = al_ . next ( ) ; // r e q u i r e Algo<T> : : ss_ t o be o f // S e e d S t a t u s <Algo> type same_type ( ss_ , al_ . ss_ ) ; // r e q u i r e Algo<T> : : ps_ t o be o f // P a r a m e t e r i z e d S t a t u s <Algo> type same_type ( ps_ , al_ . ps_ ) ; } As a result, when users try to provide an Algorithm implementation to instantiate a particular RNG, the input Algorithm class is checked for the expected requirements list. Beyond the fact of proposing a generic uniform interface to RNGs, we ensure here that future implementations respect the standard we propose. You need to puzzle out that these constraints are here for information purpose. Users become aware of some lacks in their implementations with clear compiling errors, rather than letting misconception errors lurk in the shadows until someone calls a missing method. To conclude, this feature is a full part of the developer-friendly aspect of our library, and its straightforward uniform API is likewise.

Generic programming inputs from the two last studied points have been integrated in the meta-model. Figure 3 shows a class diagram based upon Figure 1, but where policies and concept checking have replaced the Strategy design pattern:

MTGP integration

MTGP is the PRNG we chose to integrate as an example for our framework. Its authors provide a CUDA implementation of their work. The latter intends to accelerate the generation of pseudorandom numbers on GP-GPU to get them back to the host side afterwards. However, we want to generate numbers that can be directly consumed by the application currently running on the GP-GPU. Thus, we had to slightly modify Saito's original proposition to change its behavior regarding generated numbers. The original version is able to generate three numbers per GP-GPU thread, reducing this way memory transfer latency between host and GP-GPU by making each thread compute three times more. In our case, threads draw numbers in order to directly feed their next operation. We cannot make them draw several numbers and then use only one: this would introduce an overtime to finally get numbers wasted. Thus, we slimly changed the original code to make it fit our utilization constraints. Please note that these changes have been done upstream from the framework integration and are not tied to this operation. Natively, MTGP issues three random number data types: integers, single precision floats and single precision floats contained in [0;1[. As we said previously, our framework implies to separate Parameters and Algorithm into different meta-classes. Algorithm depends on a data type, whereas Parameters depends on an Algorithm. This template machinery allows us to write MTGP code once by abstracting the output data type. We will then be able to use the MTGP algorithm through the RNG meta-class, without paying attention to the involved data type. Outside our framework, the same code would have been duplicated three times, once per supported data type. Parameters sometimes depend on the handled data type. Luckily, templates enable us to cover this kind of situation. Thanks to our meta-model that manages separately parameters and algorithms, and to the template specialization mechanism, we are able to treat particular parameters by specializing the template ParameterizedStatus class to provide an implementation dedicated to the Algorithm and the involved data type.

Considering the MTGP template class that implements the aforementioned RNG, the following code snippet gives an insight of ShoveRand's capacities and ease of use. A simple CUDA parallel program, a kernel, declares and uses MTGP to compute the product of two random numbers:

CONCLUSION

We benchmarked two widespread libraries, furnished by NVIDIA research teams. The first one, CURAND, is dedicated to random number generation on GP-GPU, whereas the second, Thrust::random, is a part of the well known Thrust, a global purpose GP-GPU library.

However, according to several statistical studies, the encapsulated RNGs are mostly stated as bad quality randomness sources, in spite of their quickness and small memory footprint. Thus, we proposed a PRNG framework designed 'à la Boost', so as Thrust::random, but encapsulating quality-proven RNGs. In this way, we integrated Saito's work MTGP introduced in [START_REF] Saito | A variant of mersenne twister suitable for graphics processors[END_REF], which is a GP-GPU-enabled parallel PRNG. It was originally written to use the GP-GPU as an hardware accelerator to furnish large amounts of random numbers in a little time. The purpose here was to use it as a regular generator for stochastic computing on GP-GPU. Part of this work implied to modify the behaviour of this PRNG to have it generate a single random number per thread that could be directly consumed by the program currently running on the device. We intend to implement now other fine generators from L'Ecuyer's team, such as MRG32k3a [START_REF] L'ecuyer | An object-oriented random-number package with many long streams and substreams[END_REF].

Another major input of ShoveRand is its concept-checking mechanism. Thanks to this feature, we ensure that every RNG integrated in ShoveRand will display a common structure. This enables us to conceive stochastic applications independently of the underlying generator. Other future components will also take advantage of this statically checked skeleton. Indeed, we plan to implement the distributions features that come with the Boost.Random library. The interesting point here is that ShoveRand's will simplify the implementation of distributions because of its concept checking capabilities. In fact, constraints added to RNG structures will force every RNG implementation to very specific class members used by distributions.

Currently developed in C++, our proposition must be compiled for every new host it will run on. This solution supposes that the user owns the development toolkit to compile our sources to his own configuration. To bypass this constraint, we wish to implement our library in an abstract language such as Java or Scala, which are both executed in the Java Virtual Machine. Thanks to the OpenCL bindings proposed for the latter languages, our library could run in any environment as long as it supports Java and OpenCL.

Interested readers can obtain source code and documentation relative to ShoveRand on our university's software forge, at: http://forge.clermont-universite.fr/projects/shoverand. Vandereruch, Clément Beitone, Sébastien Cipière, Florian Guillochon and Nicolas Jean, for their careful reading and useful suggestions on the draft.

2 Figure 1 :

 21 Figure 1: Polymorphic Version of the Model

Figure 2 :

 2 Figure 2: Policy's Interface is Checked Through Boost Concept Check Library

Figure 3 :

 3 Figure 3: Meta-Model Describing the ShoveRand Framework

Figure 4 :

 4 Figure 4: MTGP Integrated in ShoveRand

ACKNOWLEDGMENTS

This work was financially supported by the regional council of Auvergne (France) in the In-novaPole/Auvergrid project. The authors would like to thank to Hélène Beaumont, Natacha

funded by the Auvergne Regional Council and the FEDER


Jonathan Passerat-Palmbach

Claude Mazel

David R C Hill

Pseudo-Random Number Generation on GP-GPU

Keywords: Stochastic Simulations, GP-GPU, Pseudo-Random Number Generators, Random Stream Parallelization

Random number generation is a key element of stochastic simulations. It has been widely studied for sequential applications purposes, enabling us to reliably use pseudo-random numbers in this case. Unfortunately, we cannot be so enthusiastic when dealing with parallel stochastic simulations. Many applications still neglect random stream parallelization, leading to potentially biased results. Particular parallel execution platforms, such as Graphics Processing Units (GPUs), add their constraints to those of Pseudo-Random Number Generators (PRNGs) used in parallel. It results in a situation where potential biases can be combined to performance drops when parallelization of random streams has not been carried out rigorously. Here, we propose criteria guiding the design of good GPU-enabled PRNGs. We enhance our comments with a study of the techniques aiming to correctly parallelize random streams, in the context of GPU-enabled stochastic simulations.

Introduction

Stochastic simulations are computationally greedy applications. Consequently, domain experts, trying to decrease over and over again their simulation execution time, look for new solutions. For about five years now, we have seen a growing interest in General-Purpose computation on Graphics Processing Units (GP-GPU) through the simulation community. GP-GPU offer a startling amount of computational power at an incredibly low cost, compared to supercomputers. The introduction of this new kind of architecture implies not only the rewriting of the simulation algorithms but also the tools they use. In this article, we have chosen to focus on Pseudo-Random Number Generators (PRNGs), which are the basis of any stochastic simulation.

Sequential PRNGs have been studied for a long time and finding a good quality PRNG to use in a sequential application has not been a problem for more than a decade. A lot of studies have been accomplished to characterize the statistical quality of PRNGs, leading to several tests batteries software. Nowadays, reference test suites are well-known and used to assess PRNGs. But, just like a prime PRNG should not be considered outside the scope of the application it feeds, we should consider different studies according to the domain implied. For instance, cryptographic applications developers ought to base their choice on the NIST test battery [START_REF] Rukhin | A statistical test suite for random and pseudorandom number generators for cryptographic applications[END_REF], whereas simulationists should use TestU01 [START_REF] Ecuyer | Testu01: A c library for empirical testing of random number generators[END_REF]. These two test batteries can be considered as a standard for their respective domains at the time of writing.

According to [START_REF] Coddington | Random number generators for parallel computers[END_REF][START_REF] Hellekalek | Good random number generators are (not so) easy to find[END_REF], a PRNG should perform well on a single processor before being parallelized. Yet, statistical quality is a necessary but not sufficient condition when selecting a PRNG to use in a parallel context: indeed, parallel streams should be independent. Thus, providing high quality random numbers becomes even more difficult when dealing with parallel architectures. We have to take into account the parallelization technique: how will we partition random streams among parallel processing elements (threads or processors)? How will we ensure the independence between parallel streams in order to prevent the involved simulations from producing biased results? The major problem concerning independence between random streams is that no mathematical proof exists to ensure it. However, some studies lay out well-known techniques to spread random streams through parallel applications [START_REF] Coddington | Random number generators for parallel computers[END_REF][START_REF]Don't trust parallel monte carlo![END_REF][START_REF] Hill | Practical distribution of random streams for stochastic high performance computing[END_REF]. They try to ensure the highest independence between random streams using different strategies. We will consider in this article how and whether these techniques could be implemented on GPU.

Apart from the parallelization technique, another point relies directly on the architecture where the involved stochastic simulations run. If we consider a GP-GPU environment, a new difficulty comes into play: harnessing the power of the device requires a rather good knowledge of GPUs. With recent programming frameworks like CUDA (Compute Unified Device Architecture) or OpenCL (Open Computing Language), merely anyone can develop applications for GPUs, but obtaining the announced performance gain implies a higher level of understanding. Most of the work and considerations exposed here rely principally on NVIDIA CUDA solutions. We are also working with the emerging OpenCL standard [START_REF]The opencl specification[END_REF], but the latter is still not robust enough in our opinion. Its current performances are slower than what you could obtain with CUDA [START_REF] Karimi | A performance comparison of cuda and opencl[END_REF]. However, we feel that this standard deserves our interest, and should take on an important part of our further work.

In this article, we will name as processing elements the elements effectively computing data in parallel. In a CUDA GP-GPU environment, threads will represent these elements, since this framework relies on a thread level logic referred to as SIMT (Single Instruction, Multiple Threads). The latter abstracts a much more standard designation known as SIMD (Single Instruction, Multiple Data). GPUs are based on this parallel architecture kind. Here, each thread must be given different data that will be computed by an identical operation. In the case of parallel stochastic simulations, we need to furnish an independent stochastic stream to each thread, in order to prevent potential bias that could be introduced otherwise. The questions are then: how are these random streams produced on GP-GPU? And more importantly: how can we ensure that they are independent?

These two main problems led us to think about design guidelines that will hopefully help developers to design better-shaped PRNGs for GPUs. Throughout these few lines, we will focus on random numbers generated and directly consumed on the GPU. In the next sections, some of the cited works have attempted to speed-up generation using the GPU before retrieving random numbers back onto the host. Current CPU-running PRNGs display fairly good performances thanks to dedicated compiler optimizations. For instance, Mersenne Twister for Graphics Processors (MTGP) [START_REF] Saito | A variant of mersenne twister suitable for graphics processors[END_REF], the recent GPU implementation of the well-known Mersenne Twister [START_REF] Matsumoto | Mersenne twister: A 623-dimensionnally equidistributed uniform pseudorandom number generator[END_REF], is only announced as being 6 times faster than the CPU reference SFMT (SIMD-oriented Fast Mersenne Twister) [START_REF] Saito | Simd-oriented fast mersenne twister: a 128-bit pseudorandom number generator[END_REF]. Some studies show that the time spent in generating pseudo-random numbers consumes at most 30% of CPU time for some "stochastic-intensive" nuclear simulations [START_REF] Maigne | Parallelization of monte carlo simulations and submission to a grid environment[END_REF], but they are very scarce. In view of the small part of the execution time used by most of the stochastic simulations to generate random numbers, it is not worth limiting GPUs usage at the generation task. Therefore, our study is not tied to the parallelization of random number generation algorithms only, and we do not propose any new PRNG in this paper. On the other hand, we share our experience of using stochasticity in GP-GPU-enabled simulations. Considering the previously evoked SIMT model, stochastic simulations with fewer divergent branches in their source code will behave better on GPU. Still, this work does not intend to list examples of simulations that will run efficiently on GPU. Instead, we focus on the parallelization techniques of pseudo-random streams used to directly feed parallel simulation programs running on GPU (applications are called kernels in the CUDA language).

In this study, we will:

• propose GP-GPU specific criteria for PRNGs design;

• suggest requirements for parallelization techniques of random streams on GP-GPU;

• study, according to the previously introduced requirements, the suitability of well-known PRNGs and random streams parallelization techniques for GP-GPU architectures.

Now, let us start with a brief reminder of GPU architecture essentials, so that we can introduce PRNGs implementation on GP-GPU.

Random number generation on GP-GPU

The purpose of this section is not to carry out a full survey of GPU-enabled PRNGs propositions in the literature. An attempt to do so can actually be found in [START_REF] Passerat-Palmbach | Reliable initialization of gpuenabled parallel stochastic simulations using mersenne twister for graphics processors[END_REF]. We focus here on PRNGs implementation level on the GPU. The latter cited survey distinguishes three implementation levels, mapped on the CUDA thread hierarchy: threads, blocks of threads and grid of blocks. These strategies directly impact the PRNG implementation, so as the parallelization technique coupled with it. Let us introduce them to understand the technical prerequisites about the two subjects we are tackling in this study. Obviously, new PRNG algorithms have to take advantage of GPU intrinsic properties. For simplicity purposes, we will briefly introduce the major concepts that rule GPU architectures, that is to say: heterogeneous memory hierarchy and thread organization. Notions described in this paragraph are represented in Figure 1. On the one hand, the thread organization is here to maximize the performances of the application. Threads are bundled into blocks of threads to be affected to one of the Streaming Multiprocessor (SM) of the GPU. SMs can mostly be considered as GPU equivalents to CPU cores. The important point here is that they have their own thread scheduler, championing threads which data are available. Selected threads then run on Streaming Processors (SP): the processing units of SMs. It is here that the matching notion, namely heterogeneous memories, comes into play. Threads can access several memories displaying various capacities, which we will discuss further in Section 3. For now, we just need to keep in mind memory areas hierarchy: i.e. the classification of memories based on their quickness and visibility from threads. From the quickest to the slowest, threads can access: registers, shared memory, local memory and global memory. When registers and local memory are dedicated to a single thread, shared memory is visible to all the threads within a block, while every thread can reach global memory. This memory organization highly affects the PRNG performances. With the first considered implementation level, using a generator per thread, the internal state of the PRNG has to be saved in each thread's local memory. Most of the time, internal states are formed by several elements contained in arrays. Now, CUDA related works, like [START_REF] Kirk | Programming Massively Parallel Processors[END_REF], specify that arrays declared within a thread are stored in the local memory. Although its name indicates a thread scope, please note that local memory is in fact a subset of global memory, and suffers consequently from the same slowness. This area is also allocated to threads when their register set is spent, since registers are available in limited quantities on GPUs.

Equivalently, with the last implementation level, a PRNG for all threads, that is to say a grid-level scope, the global memory is solicited to store the state of the unique PRNG of the application. Each thread draws a number and updates its component of the state in global memory. These two approaches make heavy use of global memory, which has the advantage of being persistent across kernel launches within the same application. Yet, this area is quite slow: it implies a 400 to 800 clock cycle latency because it is not cached [15]. The second implementation scope, the block level, is the only one left to discuss. Every thread in a block can access a shared memory area. PRNG algorithms can consequently store their internal state in this space, enabling every thread of the block to update it. Shared memory is implemented on-chip and is consequently as fast as registers. Thus, PRNGs implemented at a block level will not suffer from the memory latency induced by slow global memory accesses. The three implementation scopes detailed in this section are sketched in Figure 2:

So far, we have described PRNG algorithms, but we also found interesting libraries propositions. We noted two main proposals in this domain: CURAND and Thrust::random. They both aim to provide a straightforward interface to generate random numbers on GPU. Intro-

!"#$%&"'()" *+)(,-" !"#$%&"'()" *+)(,-" !"#$%&"./)",00"*+)(,-1" !"#$%&"'()"20/34" /."*+)(,-1"

!"# The main advantage of CURAND is that it is able to produce both quasi-random and pseudorandom sequences, either on GPU or on CPU. The quasi-RNG and pseudo-RNG implemented are respectively Sobol and XorShift [START_REF] Marsaglia | Xorshift rngs[END_REF]. The Application Programming Interface (API) of the library stays the same no matter which kind of RNG is selected and the platform on which the application is run on.

$ %# $ &# $ '# $ (# $ )# $ *# $ +# $ ,# !"#$%& ' ("%#$&' )$)*#+ ' ,*-%,' )$)*#+ ' .,*/%,' )$)*#+ '
Thrust::random is part of a GPU-enabled general purpose library called Thrust [START_REF] Hoberock | Thrust: A parallel template library[END_REF]. This open-source project intends to provide a GPU-enabled library equivalent to standard generalpurpose C++ libraries, such as STL or Boost. Classes are split through several namespaces, of which Thrust::random is an example. The latter contains all classes and methods related to random numbers generation on GP-GPU. Thrust::random implements three PRNGs, each through a different C++ template class. We find a Linear Congruential Generator (LCG), a Linear Feedback Shift (LFS) [START_REF] Tausworthe | Random numbers generated by linear recurrence modulo two[END_REF] and a Subtract With Borrow (SWB) [START_REF] Marsaglia | A random number generator for pc's[END_REF][START_REF] Marsaglia | A new class of random number generators[END_REF]. Although the latter PRNG is mentioned as Subtract With Carry in Thrust::random documentation, Marsaglia's original proposition is known as SWB. In spite of the known flaws laid out by all these generators, the library offers simple ways to combine them into better quality randomness sources, like L'Ecuyer's Tausworthe combined generators [START_REF] Ecuyer | Maximally equidistributed combined tausworthe generators[END_REF].

GP-GPU specific requirements for PRNGs

As a result of the SIMD parallelism between threads and of their graphics processors legacy, which we superficially discussed in Section 2, GP-GPUs are not equivalent to a set of standard processors used in parallel. These particularities introduce some constraints that need to be satisfied if we do not want to see the overall simulation performance dramatically drop. Thus, we will introduce in this section the requirements we find compulsory for a PRNG to run efficiently on a GP-GPU architecture.

The main goal targeted when using GP-GPUs is to obtain greater speed-ups. But this kind of device has not been primarily designed to support general computations and is more inclined to realize some arithmetic operations. Nowadays GPUs display different performances with single and double precision floating point numbers. For instance, the previous generation of NVIDIA supercomputing-dedicated GPU, the Tesla T10, was known to deliver ten times as much computation power when dealing with single precision floating point numbers rather than double. Even if the actual cutting-edge GPU, the NVIDIA new generation Fermi T20, has considerably reduced the gap between these two precisions, it would be wise to remain cautious before using double precision operations on GPU. Most of the time, single precision floats are sufficient enough to handle random values contained in [0 ; 1[ and should consequently be favoured. This proposition leads us to our first criterion: single precision floating point numbers should be preferred throughout the GP-GPU random number generation algorithm.

Another legacy of graphics processors is the heterogeneous memory organization. To complete what has previously been said on this subject, let us recall that several memory areas are reachable by threads running on a GPU. These memories capacities, some quick and large, others less so, depend on two characteristics. First, the more threads can reach a memory area, the slower it is. Indeed, registers allocated to a single thread are the quickest memory this thread will be able to communicate with. At the same efficiency level, we find shared memory, reachable by a relatively small amount of threads, all belonging to the same block, and so running on the same core of the GPU. On the other hand, every thread running on the GPU, regardless of which core they are located on, can access a wide memory area, commonly called global memory. Therefore, this latter area is far slower than its previous counterparts. Second, read-only memories are quicker since they can fully benefit from cache mechanisms, contrary to read-write memories. In the light of these memory constraints, it is obvious that GPU PRNGs should be designed with particular attention to sparing costly memory accesses. Commonly, static parameters will take place in read-only areas, whereas dynamic elements such as state vectors will be handled at the thread or thread group level. In a more formal way, the following is another criterion of good design: the algorithm should be designed in a way to avoid global memory accesses.

Taking into account the particularities of GPUs and their architecture, we have proposed two new requirements for PRNGs to run efficiently on such devices. They are summed up hereafter:

1. single precision floating point numbers should be preferred throughout the GP-GPU random number generation algorithm 2. the algorithm should be designed in a way to avoid global memory accesses

GP-GPU specific requirements for random streams parallelization

The literature is full of references describing the profile of what a good usage of parallel PRNGs should be. For example, [START_REF] Coddington | Random number generators for parallel computers[END_REF][START_REF] Hellekalek | Good random number generators are (not so) easy to find[END_REF] list requirements that should meet any sequential or parallel PRNG. GPUs are particular parallel architectures. So any PRNG running on this kind of device should, at least, meet the requirements enumerated in the previous references. In this section, we will successively check how these criteria can be adapted to GPU architectures. Emphasizing parallel PRNG performances, [START_REF] Coddington | Random number generators for parallel computers[END_REF] noticed that each processor should generate its sequence independently of the other processors. We consider, indeed, that every processing element should have its own stochastic stream at its disposal. This condition must be satisfied first, not only for efficiency, but especially because GPU-enabled stochastic simulation parallelization principles rely on it. First, it is a necessary, but not sufficient, condition to fulfil to ensure a higher independence of stochastic streams feeding different replications of a simulation. Second, considering a single replication, the SIMT parallelism level leads threads to compute their own data sets, including their own stochastic stream. Thus, our first requirement concerning random streams parallelization can be expressed as follows: each thread should dispose of its own random sequence.

As we explained previously, GP-GPU programming frameworks offer a thread scope rather than a processor one. The threads in use for GP-GPU propose an abstraction of the underlying architecture. They are concurrently running on the same device and handle their own local memory area. Thread scheduling is at the basis of GPU performances. Memory accesses are the well-known bottleneck of this kind of device. Indeed, running a large amount of threads in turn allows GPUs to bypass memory latency. There should always be runnable threads while others are waiting for their input data. Disregarding the effective number of processors, we theoretically say that the more threads you have, the better your application will leverage the device. Applications need to be written to use the maximum number of threads, but also to scale up transparently when the next GPU generation will be able to run twice as many threads as today. So, in accordance with [START_REF] Coddington | Random number generators for parallel computers[END_REF] who advocates that the generator should work for any number of processors, our second GP-GPU specific requirement for parallelization techniques of random streams is that it must be usable for any number of GP-GPU threads.

Returning to the original requirements, we then find [START_REF] Coddington | Random number generators for parallel computers[END_REF] states that parallel random streams produced should be uncorrelated. This criterion is related to both PRNG intrinsic properties and to the parallelization technique set up. We previously stated that a PRNG candidate to parallelization should first perform well on a single processor. Thus, we will not take its intrinsic qualities into account here. However, no matter the worth of the used PRNG, the parallelization techniques must be used carefully, as will be shown in the next section. Please note that this requirement is neither affected by GP-GPU architectures nor by programming frameworks. As a result, we will just recall it without modifying its expression.

[3] also noted that the same sequence of random numbers should be produced for different numbers of processors, and for the special case of a single processor. Here, we understand that the PRNG output on each processor should not depend on the number of processors used. This point is very important and must be treated carefully when choosing a parallelization technique. For example, in a distributed environment containing several processors, a scheduler can govern execution. Depending on the scheduler algorithm and on the global system charge, parallel executions of different parts of a simulation might not execute in the same order. In such a case, it is compulsory for the PRNG output to be independent from the order in which simulations parts may run. If this requirement is not met, reproducibility of simulations executions is no longer ensured. We can do this for games, but not for scientific applications. Reproducibility is needed when dealing with stochastic simulations, in order to debug a problematic case raised by a particular random stream for instance. We also think about design of experiments for simulations, where reproducibility is mandatory to isolate parameters variations impact in the results.

In the case of GP-GPU, we find exactly the same problem at the thread level. These entities are also scheduled, not atomically but by bundles. Fortunately, both threads and their bundles own a unique identifier allowing us to distinguish them among executions. Thus, if a parallel random stream is only bound to the unique identifier of a thread, according to our first requirement, output will be reproducible through multiple executions. Therefore, we can obtain the necessary bijective relation between a T i thread and an SS i stochastic stream, as stated in Figure 3:

Finally, we can write our last requirement in such a way: when the status of the PRNG is not modified, the sequence of random numbers generated for a given thread must be the same no matter the number of threads and no matter of threads scheduling.

Let us now sum up the requirements targeting GPUs we highlighted in this part, continuing the enumeration started in Section 3:

each thread should dispose of its own random sequence

! "# ! $# %# ! &# %# ! '# (( "# (( $# %# (( &# %# (( '#

the parallelization technique must be usable for any number of GP-GPU threads

5. the parallel random streams produced should be uncorrelated 6. when the status of the PRNG is not modified, the sequence of random numbers generated for a given thread must be the same no matter the number of threads and no matter of threads scheduling

Do classical random streams parallelization techniques fit GP-GPU well?

This section presents the main techniques used to distribute random streams between processing elements. A more exhaustive and dedicated survey can be found in [START_REF] Hill | Practical distribution of random streams for stochastic high performance computing[END_REF]. We will also study in this section how random streams parallelization techniques can be adapted to GPU architectures, depending on their ability to fulfil the previously introduced requirements (Sections 3 and 4). The Leap Frog is the way to partition a random stream like a deck of cards. Random numbers are allocated in turn to processors like cards would be dealt to players. Pragmatically, let each processor hold an i identifier. The latter processing element will build a Y i substream from an X original random stream such as Y i = {X i , X i+N , . . . , X i+kN }, with N equal to the number of processors. This technique is not quite adapted to split random streams on GP-GPU since it does not satisfy the last constraint expressed in Section 4. In fact, if the number of threads changes, the subsequence affected to each thread will be different. This situation is shown in Figure 4. Now, the number of threads for an application is bound to the underlying device: GPUs can run a different number of threads concurrently, depending on their architecture generation. As a result, we would not be able to ensure the reproducibility of a simulation from a GPU to another, even if we initialized the PRNG with the same parameters. Furthermore, some bugs induced by random draws might not appear on developer's device, while they would on the user's.

!"#$%&'()" * +" * ," * -" * ." * /" * 0" * 1" * 2" * 3" * 4" * ,+" * ,," -"56('78)" /"56('78)" * +" * ," * -" * ." * /" * 0" * 1" * 2" * 3" * 4" * ,+" * ,," The solution would be to implement the PRNG at a level with a constant thread number. The CUDA framework handles constant-sized bundles called warps, which always contain 32 threads at the time of writing. They are in fact a subdivision of blocks of threads, and access consequently the same shared memory area. Nowadays, a great number of GPUs do not have enough shared memory to store a PRNG status per warp. However, the newest GPU generations offer larger shared memory spaces that could potentially enable us to use Leap Frog, following this idea.

Apart from any architecture consideration, please note that Leap Frog is not adapted to every kind of PRNGs. [START_REF] Hellekalek | Good random number generators are (not so) easy to find[END_REF] has established that this technique could be flawed in some cases, because of spectral tests. We particularly insist on what an unsafe technique it is when used with LCGs. Serious stochastic simulations should ban simple LCGs since they have structural flaws [START_REF]Encyclopedia of Quantitative Finance[END_REF].

The second technique to be referenced in [START_REF] Coddington | Random number generators for parallel computers[END_REF] is Sequence Splitting (also known as fixed spacing or blocking). It consists in allocating continuous and equally sized blocks from the original random stream to form substreams. This technique implies the knowledge of how many numbers each thread will consume at most. Indeed, knowing that each thread consumes at most L random numbers, then the first L numbers will be attributed to the first thread, the L following to the second thread and so forth. Following the previous formalism, we have Y i = {X iL , X iL+1 , . . . , X (i+1)L-1 }. This numbers repartition is described in Figure 5.

!"#$%&'(')'

!"#$%&'('*' !'#+,-"./' 0 1' 0 2' 0 3' 0 4' 0 )' 0 5' 0 *' 0 6' 0 7' 0 8' 0 21' 0 22' 0 1' 0 2' 0 3' 0 4' 0 )' 0 5' 0 *' 0 6' 0 7' 0 8' 0 21' 0 22' Efficient Sequence Splitting relies on a particular feature of the PRNG called Jump Ahead, or Skip Ahead. Instead of unrolling the generator, it enables us to compute any given advanced state of the random sequence without processing the previous ones (for example X (i+1)L can be computed directly from X iL ). Here, we discern two categories of algorithms. Some PRNGs contain an algorithm able to realize Jump Ahead, thanks to intrinsic properties of the PRNG [START_REF] Haramoto | Efficient jump ahead for f2-linear random number generators[END_REF]. This allows us to reach any part of the sequence in equal time, regardless of the destination point. Let this technique be considered as Intrinsic Jump Ahead hereafter. The other solution is to emulate Skipping Ahead: to do so, we have to compute an advanced state by processing step by step the previous ones (for example, starting from X iL , X (i+1)L can be computed by running the PRNG L times in a sequential way in order to get X iL+1 , . . . , X (i+1)L-1 and finally X (i+1)L ). In fact, whichever PRNG you use, you can unfold the sequence to the desired point, and store the state vector at this point in order to be able to load it later. Such a vector is named seed status in [START_REF] Passerat-Palmbach | Reliable initialization of gpuenabled parallel stochastic simulations using mersenne twister for graphics processors[END_REF], since it is able to set a generator in a predefined state. Emulated Jump Ahead can become very costly though: indeed, the further you need to go in the sequence, the more time it takes to compute the seed status.

To conclude on Sequence Splitting, if an intrinsic Jump Ahead algorithm is available for the involved PRNG, Sequence Splitting is a very good approach for GP-GPUs. However, as far as we know there are few ports of algorithms from this family to GP-GPU. On the other hand, Emulated Jump Ahead is not GPU-compliant because statuses computation is a purely sequential operation (we need X n to compute X n+1 ). As a consequence, different threads may display different lengths of time for these computations, thus decreasing the SIMD parallelism overall gain. To solve this problem, we propose to pre-compute substreams on the host side, store the seed status at each substream starting point and then transfer all these statuses to the device. In the past, Sequence Splitting was also considered as unsafe by several studies [START_REF] De Matteis | Computation of critical distances within multiplicative congruential pseudorandom number sequences[END_REF][START_REF] De Matteis | Critical distances in pseudorandom sequences generated with composite moduli[END_REF], but this remark applies to any PRNGs that showed potential long-range correlations, such as congruential ones.

We have just seen that emulated Jump Ahead led to expensive calculations in order to compute an initial status for each sub-sequence. Another approach based upon probabilities is known as Random Spacing. The idea is to avoid Emulated Jump Ahead cost by choosing random statuses in the base stream, and set them as initial statuses for the resulting substreams as represented in Figure 6. Although Random Spacing spares us from solving the Jump Ahead problem, it should not be used without care. [START_REF] Wu | Parallel use of multiplicative congruential random number generators[END_REF] have given the minimum distance between N sub-sequences which status were randomly chosen: it is in average equal to 1/N 2 multiplied by the period length. More precisely, the probability of overlapping between N sequences of length L, issued from a PRNG of period P is equal to: 1 -(1 -N L/(P -1)) (N -1) . That is equivalent to N (N -1)L/P when N L/(P -1) is in the neighbourhood of 0.

!"#$%&'()" * +" * ," * -" * ." * /" * 0" * 1" * 2" * 3" * 4" * ,+" * ,," * ,-" * ,." * ,/" * ,0" Overlapping risks become sizeable with short period PRNGs. All the PRNGs tested by Pierre L'Ecuyer and Richard Simard in [START_REF] Ecuyer | Testu01: A c library for empirical testing of random number generators[END_REF] display periods P from 2 24 to 2 131072 . Most of them have log 2 P of the order of a few dozens. Now, given that nowadays longest simulations can consume up to several hundreds of billions of random numbers, (L = 10 11 ), a hundred of such replications (N = 100) leads to a N (N -1)L of approximately 10 15 , i.e. 2 50 . The probability to see an overlapping between two sub-sequences issued from a PRNG of period P far bigger than 2 50 is negligible. Nonetheless, from 20 LCGs PRNGs tested in [START_REF] Ecuyer | Testu01: A c library for empirical testing of random number generators[END_REF], 14 laid out an overlapping probability greater than 99.9% (period less than 2 47 ). Only 3 generators have an overlapping probability less than 0.01% (log 2 P > 59.7). In the same way, widely spread PRNGs such as Comblec88 of period P = 2 61 (Combined LCGs), are on the acceptance borderline for this technique. They are shipped with several renowned software packages though, including: RANLIB, CERNLIB, Boost, Octave and Scilab [START_REF] Ecuyer | Testu01: A c library for empirical testing of random number generators[END_REF].

As a conclusion, this technique fits GP-GPUs well, but the risk of overlapping between subsequences must be evaluated according to the number and the length of the sub-sequences and to the period of the PRNG used. If we select generators with large periods, such as WELLs and Mersenne Twister, this risk is negligible.

Techniques presented so far tried to split a single stream into several substreams. Another approach consists in using several declinations of the same PRNG: each generator has the same structure and generation mechanism with a unique parameter set, called parameterized status hereafter. By allocating a parameterized status per processor, parameterization ensures our first requirement, namely that each thread should dispose of its own sequence. This situation is described in Figure 7. We always have to remember that this cannot ensure our third requirement: strict independence between parallel random streams.

Although no mathematical proof can establish this independence, some implementations of Parameterization are safe according to the current state of the art [START_REF] Mascagni | Parameterizing parallel multiplicative lagged-fibonacci generators[END_REF]. We especially think to the Dynamic Creator (DC) algorithm [START_REF] Matsumoto | Dynamic creation of pseudorandom number generators[END_REF] coming along with the Mersenne Twister for Graphics Processors (MTGP) PRNG [START_REF] Saito | A variant of mersenne twister suitable for graphics processors[END_REF], a quality-proven work benchmarked in [START_REF] Passerat-Palmbach | Reliable initialization of gpuenabled parallel stochastic simulations using mersenne twister for graphics processors[END_REF]. DC integrates a unique identifier, which belongs to the parameterized status of MTGP. This identifier becomes a part of the characteristic polynomial of the matrix that defines the recurrence of the PRNG.

Figure 7: Three parameterized PRNGs producing three highly independent random sequences Two identifiers will consequently lead to two different parameterized statuses. Furthermore, DC ensures that the characteristic polynomials we obtain are mutually prime and the authors assert that the random sequences generated with such distinct parameterized statuses will be highly independent, even if, as mentioned before, this fact cannot be mathematically proven.

This technique displays some constraints making it difficult to port on GPUs. Unfortunately, few PRNGs propose it intrinsically. Some, such as LCGs, are even recognized as bad candidates for Parameterization [START_REF] De Matteis | Parallelization of random number gemerators and longrange correlations[END_REF]. In addition, storing seed statuses being already problematic, we can scarcely imagine spending vast amounts of memory to store a parameterized status per thread.

Every technique introduced in this section, presents advantages and drawbacks. Most of them are related to the chosen PRNG. Depending on the application and environment you own, you might be forced to select a PRNG knowing it has some flaws in particular cases. In this way, Table 1 states PRNG kinds and parallelization techniques that work well together:

Conclusion

We have seen that more and more applications, and especially stochastic simulations, tend to take advantage of recent GP-GPU architectures in order to improve their performance. However GPU computing needs to offer the same tools as other platforms. Good quality PRNGs belong to this category. This paper has shown how difficult it could be to generate good quality pseudo random numbers on GP-GPUs. Indeed, it implies taking into consideration two different domains: GP-GPU programming and PRNG parallelization techniques. Issuing a PRNG that can produce independent stochastic streams when used in parallel is a first hurdle that not all PRNGs can get over. When you have at your disposal one that fulfills this requirement, it has to be ported to GP-GPU. It means that you need to think about a GPU parallelization of your PRNG, if it is not already available. This paper introduced the main problems that you will be faced with when trying to port your stochastic simulations on GP-GPU. To avoid these difficulties, and above all, errors and performance drops that could result from the use of a hazardous GPU-enabled PRNG, we first proposed PRNGs criteria dedicated to GPUs. Then we defined parallelization techniques requirements in order to make these PRNGs fit GPU peculiarities. They sum up the experience accumulated in our research team concerning GPU-enabled stochastic simulations.

We encountered lots of parameters brought up by PRNGs and GPU programming. As long as it can dramatically impact both the overall performance of the simulation and the quality of its results, it might be a good point to propose a straightforward API to use well-defined PRNGs on GPUs. In this way, libraries laid out in Section 2 represent in our mind an elegant manner to do so. They allow every user to easily call PRNG functions without wasting time on how parameters should be set. In the same way, they enable advanced users to adjust parameters in their own fashion. As a matter of fact, further works of ours should target this kind of development, embedding PRNGs following the presently established requirements into a GPU-enabled library with a unified API.

!Figure 1 :

 1 Figure 1: Simple representation of the major elements of a GPU

Figure 2 :

 2 Figure 2: PRNGs implementation scopes and their location in the different memory areas of a GPU

Figure 3 :

 3 Figure 3: Bijective relation between threads and stochastic streams

Figure 4 :

 4 Figure 4: Different threads numbers leading to different random substreams through the Leap Frog method

Figure 5 :

 5 Figure 5: Two random streams parallelizations based upon Sequence Splitting with two different sub-sequences lengths

Figure 6 :

 6 Figure 6: Random Spacing creation of three sub-sequences of equal length but differently spaced from each other

Table 1 :

 1 Summary of the potential PRNG/Parallelization technique associations

	Technique	Preferred	PRNGs
		PRNGs	to avoid
	Leap Frog	None
		(disable reproducibility) Linear generators
	Sequence Splitting	Intrinsic Jump-Ahead	Emulated
		compliant	Jump-Ahead
	Random Spacing	Large period	Short period
	Parameterization	MT family	LCG

Acknowledgements

Special thanks to Susan Arbon, Natacha Vandereruch and Florian Guillochon for their careful reading and useful suggestions on the draft.

funded by the Auvergne Regional Council and the FEDER