How to generate a massive amount of high quality Random Numbers?

Question

I'm working on a random walk simulation of particles moving in a lattice. For that reason I must create a massive amount of random numbers, about 10^12 and above. Currently I'm using the possibilities C++11 provides with <random>. When profiling my program, I see that a major amount of time is spent in <random>. The vast majority of those numbers are between 0 and 1, evenly distributed. Here a then I need a number from a binomial distribution. But the focus lies on the 0..1 numbers.

The question is: What can I do to reduce the CPU time needed to generate these numbers and what would the impact be on their quality?

As you can see, I tried different engines, but that had no big effect on CPU time. Further, what is the difference between my uniform01(gen) and generate_canonical<double,numeric_limits<double>::digits>(gen) anyhow?

Edit: Reading through the answers I conclude that there is not THE ideal solution for my problem. Thus I decided to first make my program multi threading capable and run multiple RNG in different threads (seeded with one random_device number + an thread individual increment). For the time being this seams to be the most unavoidable step (multi threading would be required anyhow). As a further step, pending on exact requirements I consider switching to the suggested Intel RNG or to Thrust. Meaning that my RNG implementation should not be to complex, which, currently is is not. But for now I like to focus on the physical correctness of my model and not on programming stuff, this comes as soon as the output of my program is physically correct. Thrust Concerning Intel RNG

Here is what I do currently:

class Generator {
public:
    Generator();
    virtual ~Generator();
    double rand01(); //random number [0,1)
    int binomial(int n, double p); //binomial distribution with n samples with probability p
private:
    std::random_device randev; //seed
    /*Engines*/
    std::mt19937_64 gen;
    //std::mt19937 gen;
    //std::default_random_engine gen;
    /*Distributions*/
    std::uniform_real_distribution<double> uniform01;
    std::binomial_distribution<> binomialdist;
};

Generator::Generator() : randev(), gen(randev()), uniform01(0.,1.), binomial(1,1.) {
}

Generator::~Generator() { }

double Generator::rand01() {
    //return uniform01(gen);
    return generate_canonical<double,numeric_limits<double>::digits>(gen);
}

int Generator::binomialdist(int n, double p) {
    binomial.param(binomial_distribution<>::param_type(n,p));
    return binomial(gen);
}

Algo

Although it is not in C++, you should definitely look into picking up a copy of the very well reviewed A Million Random Digits with 100,000 Normal Deviates. I was so happy with my first copy that I bought a second one. – Chargeable

@HostileFork: I'm fighting with myself to not give you +1, as your comment is in fact not helpful :) – Scholem

It could be compiler, operating system and hardware specific. On Linux you might use /dev/random at least to seed regularily your PRNG. You could then choose a simpler PRNG (so faster, but worse) like drand48 – Cribriform

I indeed thought the same... because 1'000'000 samples are not enough ;) – Algo

You could have an auxiliary file (with real random numbers) where you use a random number to index which block of the file you are going to use. Maybe you could have more than one index and combine the blocks between them self , like xoring. I am not a specialist, but if I am not wrong, xoring two real random numbers give you another real random number. – Hereabout

@AndréPuel: given typical floating point binary representations, a naive XOR wouldn't work very well - particularly against the criteria to constrain output to the 0..1 interval - but a combination of XORing the bits in the mantissa with a weighted selection of new exponent has potential. – Miscreance

@TonyD I was thinking in fixed point representation actually, you get the 4 bytes unsigned integer and divide by double(0xFFFFFFFFu) – Hereabout

@AndréPuel ok - that'd work if 32 bits of precision is adequate (double supports 52) and a good way to permute the random numbers adequately before reuse. In your comment above I took "real random numbers" to mean floating point given the mathematical meaning of "real" - sorry for the confusion. – Miscreance

I do not guaranty the speed improvment, but for a flat random distribution in [0,1], you could use exponential inside a sine or cosine function. The exponential slope soon overrides the sine period and you start to have random-like numbers for each x increment by one. I used it once long long ago, but I don't remember if it had speeded calculations. (It was one between several recipes for generating pseudo random numbers and chaotic values). – Mumble

Thinking laterally, why not use an extremely fast hasher as your random source? Start with any document as you like as the seed and use hash(N) as the input to get hash(N+1). xxhash promises 5Gb/sec. Is that fast enough for you? – Surefooted

Intel's new CPU has builtin random number generation instruction. Don't know if you can use it though. – Meldameldoh

@BasileStarynkevitch: urandom please. random has no legitimate use that I'm aware of. – Sendal

FYI, you don't need to keep around the random_device; just initialize gen(std::random_device()()) – Steeplechase

oblig xkcd.com/221 – Berners

Answer 1

You can pre-process random numbers and use them when you need.

If you need true random numbers I suggest you to use a service like http://www.random.org/ that ensures random numbers calculated by environment ambient instead that some algorithm.

And, speaking about random numbers, you must also check this:

Spry

Answer 2

If you need a massive amount of random numbers, and I mean MASSIVE, do a careful search on the internet for IBM's floating point random number generator, published maybe ten years ago. You'll have to buy either a PowerPC machine, or a newer Intel machine with fused multiply-add. They achieved random numbers at a rate of one per cycle per core. So if you bought a new Mac Pro, you could achieve probably 50 billion random numbers per second.

Glengarry

Answer 3

Perhaps instead of using a CPU you could use a GPU to generate many numbers concurrently?

Efficient Random Number Generation and Application Using CUDA

Gothicize

Answer 4

On my i3, the following program runs in about five seconds:

#include <random>
std::mt19937_64 foo;

double drand() {
  union {
    double d;
    long long l;
  } x;
  x.d = 1.0;
  x.l |= foo() & (1LL<<53)-1;
  return x.d-1;
}

int main() {
  double d;
  for (int i = 0; i < 1e9; i++)
    d += drand();
  printf("%g\n", d);
}

whereas replacing the drand() call with the following results in a program that runs in about ten seconds:

double drand2() {
  return std::generate_canonical<double,
      std::numeric_limits<double>::digits>(foo);
}

Using the following instead of drand() also results in a program that runs in about ten seconds:

std::uniform_real_distribution<double> uni;
double drand3() {
  return uni(foo);
}

Perhaps the hacky drand() above suits your purposes better than the standard solutions..

Sendal

Answer 5

Task Definition

OP asks to get answer for both the

1. Speed of generation -- assuming a set of 10E+012 random numbers to be "massive"

and

2. Quality of generator -- with a weak assumption that numbers just evenly distributed over some range of values are also random

However, there are more cardinal aspects to be addressed and successfully solved for the real system:

A. Define, whether your system simulation needs to be provided with a guarantee of a repeatability of the sequence of the random numbers for future re-runs of an experiment.

If this is not the case, the re-runs of the simulated experiment will yield principally different results then the randomizer process ( or pre-randomizer and randomized-selector ) need not worry about their re-entrant, state-full mode of operation and will get much simpler implementation.

B. Define, to what level do you need to proof a quality of randomness of the generated random numbers ( or does the generated sets of random numbers have to belong to some specific law of statistic theory ( some known synthetic distributions or truly random with an utmost Kolmogorov complexity of the resulting set of random numbers )). One need not be NSA expert to state that numerical generators of true-random sequences is a very hard issue and has it's computational costs associated with production of high-randomness products.

Hyper-chaotic and true-random sequences are computationally extemely expensive. Using low- or poor-randomness generators is not an option for randomness-quality sensitive applications ( whatever the marketing papers may say, no MIL-STD- or NSA-graded system will ever try this compromised quality in enviroments, where the results indeed matter, so why to settle for less in scientific simulations? Perhaps not a problem if you do not mind to miss so many "unvisited" states of the simulated phenomena ).

C. Verify, how many random numbers does your simulation system need to "consume per [usec]" and whether this design requirement parameter is constant or may get scaled-up by going into multi-threaded, vectorised, Grid-/Cloud-based distributed computation framework.

D. Does your simulation system require to maintain a global or per-thread- or perGrid/CloudNode- individual access management to the pool-of-randomized numbers in case of vectorized or Grid/Cloud-based computational strategy.

Task Solution Approach

Fastest [1] and best [2] solution with [A] and [B] solved and options for [D] is to pre-generate an utmost randomness quality numbers into an adequate access-pool ( and pay an acceptable cost of [C] and [D] on access-policy and access-management controls to re-read from the pool, rather than to re-generate ).

Pelaga