I was recently starting to write a monte carlo simulation involving part survival, and got some very strange results. For one single time-step, 30 time-steps after the start of the simulation, parts broke less often. I made the probability of failure the same for every time-step and got the following plot. Note the large downward spike at 30 in the otherwise smooth decaying exponential curve.

Since I was still writing code during the quick-and-dirty phase, just to get something started, I was using the POSIX system call random() to generate my random numbers. Now I have always known that the system-provided random numbers are awful, but I had never actually encountered such a striking failure. Here is the code after it was hacked down to eliminate errors on my part:

```
#include <cstdio>
#include <cstdlib>
int
main()
{
const int seed = 1 ;
const double p_m = 0.05 ;
const int NDIST = 200 ;
const int NRUN = 1000000 ;
srandom(seed) ;
static int dist_m[NDIST] ;
for( int d = 0 ; d < NDIST ; ++d )
dist_m[d] = 0 ;
for( int run = 0 ; run < NRUN ; ++run )
{
int M = 0 ;
for(;;)
{
double r_m = double( random() ) / double( RAND_MAX ) ;
if( r_m < p_m )
break ;
else
M++ ;
}
dist_m[ M < NDIST ? M : NDIST-1 ]++ ;
}
for( int d = 0 ; d < NDIST ; ++d )
{
if( dist_m[d] )
printf( "%3d %f\n", d, dist_m[d] / double(NRUN) ) ;
}
}
```

Replacing the POSIX random() call with something better fixed the problem nicely.

## Comments

There are currently no comments on this article.

## Please leave a comment