There is a lot of literature that indicates using rand() for simulations might be harmful. A couple of them are as follows:

A common recommendation has been to use the ran1/ran2/ran3 functions as defined in Numerical Recipes.

But these seem pretty old criticisms of the rand function. So my questions is, how do its current implementations (GCC/clang) measure up for the purpose of being a good RNG in, say Monte Carlo simulations?

This is the current implementation of the function used by rand() in glibc: (rand calls random, which calls random_r)

  • 1
    In order to answer this question you should at least provide some objective metric by which the implementations can be measured. What about rand is important to you ? The speed ? the distribution of the random numbers ? – Bruce Becker Nov 5 at 9:25
  • @BruceBecker The most important thing for me would be the non-existence of "correlations" in the random number sequence. The reason I did not add this is because it is very ill-defined: there are an infinite number of correlation metrics definable. – physkets Nov 5 at 10:46

But these seem pretty old criticisms of the rand function

This may be a nitpick, but I want to point out what I think is a flaw in this logic. Compilers are often extremely conservative about changing program behaviour, even when that behaviour (foolishly) depends on implementation details.

This may or may not be true for the big compilers you’re familiar with, but it is conceivable that a compiler would keep a broken legacy RNG for compatibility reasons.

The following program on my macbook:

#include <cstdlib>
#include <iostream>
using namespace std;

int main() {
  srand(1);
  for (int i = 0; i < 20; ++i) {
    cout << (i > 0 ? ", " : "") << rand();
  }
  cout << endl;
  return 0;
}

produces the output

16807, 282475249, 1622650073, 984943658, 1144108930, 470211272, 101027544, 1457850878, 1458777923, 2007237709, 823564440, 1115438165, 1784484492, 74243042, 114807987, 1137522503, 1441282327, 16531729, 823378840, 143542612

which you can look up in OEIS: https://oeis.org/A096550, it is the sequence $$16807^n\bmod(2^{31}-1).$$

How this does under the SmallCrush RNG test (in the TestU01 library):

========= Summary results of SmallCrush =========

 Version:          TestU01 1.2.3
 Generator:        ulcg_CreateLCG
 Number of statistics:  15
 Total CPU time:   00:00:08.05
 The following tests gave p-values outside [0.001, 0.9990]:
 (eps  means a value < 1.0e-300):
 (eps1 means a value < 1.0e-15):

       Test                          p-value
 ----------------------------------------------
  1  BirthdaySpacings                 eps
  2  Collision                        eps
  6  MaxOft                           eps
 ----------------------------------------------
 All other tests were passed

As far as I know there is simply no good reason to assume that rand() might be okay. If you assume otherwise, you'll usually end up writing some extremely non-portable code that will break under different compilers on different systems.

  • Okay; even though I'll have control over what compiler I use, I suppose I should just use a standard rng implementation like ran1 or Ranq1 from numerical recipes. – physkets Nov 9 at 5:32
  • 1
    @physkets I would check its results on BigCrush before using it. – Kirill Nov 9 at 8:32
  • I ran Ranq1 through BigCrush, and it failed one test: ClosePairs mNP2S, t = 16, with a p-value of 8.7e-4. – physkets yesterday
  • @physkets that’s not terrible, but I believe there are much more standard/common PRNGs that might be both faster and also passing all the tests. – Kirill yesterday
  • 1
    True, but I think this is fast enough, and the main operational part of it is just 3 statements! – physkets yesterday

Though the question is about C, I still consider the following information regarding C++11 relevant. In C++11, <random> was added that is both designed to fix many of C rand() implementation flaws, as well as to provide additional capabilities. So, if you have an option of using C++11, I would point your attention to <random> and its benefits.

In most cases, <random> function family would provide sufficient quality random numbers (as opposed to rand(), which is mentioned in the materials you linked in your question. The critique to rand() is still relevant (as ANSI standard did not change); however, certainly exaggerated.

My advice would be to first assess the need in higher quality random numbers and go step-by-step. There is a high chance, that for particular purposes <random> or even rand() will provide decent enough results. Otherwise (unlikely), you might have to go to special libraries or even hardware generators.

  • I am aware of <random> in C++, but I'm restricted to C. rand() supposedly has some issues that makes it unsuitable for simulations. I'm currently checking out the Ranq1() function from Numerical Recipes. – physkets Nov 6 at 8:01

There is a finite number of correlation metrics relevant to your problem which you ideally should know before trying to implement MC!

Regarding rand(): I have written (parallel) MC methods which could easily draw $>10^9$ PRNs at which point you cycle the sequence of PRNs provided by rand(), whose period is $2^{32}-1\approx 4.3\cdot 10^9$ if I understand the comments in the code you link to correctly. Now you either have to analyse whether potential correlations will impact your problem, or you can choose a PRNG with a longer period, thereby avoiding the problem (or at least pushing it off to larger problem sizes).

Especially in parallel environments (the common way to run a MC simulation) the usage of rand() harbors the additional problem that multiple processes may end up initializing the PRNG to either start at the same position, or at positions at which their sequences will eventually overlap. The latter will result in two processes using the same PRNs for (supposedly) independent simulations. Other PRNGs, like the Mersenne Twister engine used in C++11 suffer the same problem, but to a lesser degree since their (limited) parameter sets produce different sequences, so as long as you do not initialize them with the same seed they will yield non-overlapping sequences.

A C code I have worked with used rand() for a parallel MC method used to approximate an inverse of a (very large) matrix. I have empirically determined that replacing it by a TRNG ( https://www.numbercrunch.de/trng/ ) - had to write a C-wrapper for it, though - measurably reduced the errors of the final result. Especially for large matrices (e.g., those that needed a lot of PRNs). Note here: Only because you are limited to work with C does not mean you can not use C++ libraries, you just have to write wrappers for the functions you need.

  • I do know what correlations I'm looking at, but a higher order correlation can easily affect behavior at criticality. So technically, an infinite number of correlations can affect my result (of-course, the higher the order, the more diminished the effect). – physkets Nov 9 at 5:34
  • If you are looking at something akin to a phase transition I would recommend to use /dev/urandom as a source of (bit)entropy instead of looking for pseudorandom number generators. There's a paper (can't find it atm, though) on the quality of urandom available. – Nox Nov 11 at 23:09

Your Answer

 

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.