Visual C++ Random Number: How do I generate random numbers of specified frequency distribution?

[Pravin Kumar]

Q: How do I generate random numbers of specified frequency distribution?

A: Visual C++ has a function to generate random numbers (integers) between 0 and predefined integer 'RAND_MAX'. This can be used to generate random numbers of integers within a range, say 'IMax' and 'IMin':

int I = IMin + rand() % (IMax - IMin);

or of double type between 'XMax' and 'XMin':

double X = XMin + rand() * (XMax - XMin) / RAND_MAX;

When you generate integer random numbers within a range, make sure that 'IMax > IMin' since a modulus operation is involved.

Generating random numbers according to a specified distribution function involves some things more. Let us say that we have to generate random numbers between 'XMin' and 'XMax' by the distribution function 'Fun(x)'. Evaluation of this function of any value between these limits does not generate errors like division by zero, square root of negative number etc.

First thing is to evaluate the minimum and maximum values of this function (say 'YMin' and 'YMax') within the specified limits 'XMin' and 'XMax'.

Generate random numbers in pairs, say 'X' and 'Y' - 'X' between 'XMin' and 'XMax' and 'Y' between 'YMin' and 'YMax'. Discard all random numbers 'X' where the value 'Y' is below 'Fun(X)'. Remaining numbers will have the required frequency distribution.

The function shown below generates a random number each time between 'xmin' and 'xmax', repeated use of which will generate random numbers of fairly accurate distribution function 'fun()'.

double Random(double (*fun)(double), double xmin = 0, double xmax = 1)
{
  static double (*Fun)(double) = NULL, YMin, YMax;
  static bool First = true;

  // Initialises random generator for first call
  if (First)
  {
    First = false;
    srand((unsigned) time(NULL));
  }

  // Evaluates maximum and minimum of function
  if (fun != Fun)
  {
    Fun = fun;
    YMin = YMax = Fun(xmin);
    for (int iX = 1; iX < RAND_MAX; iX++)
    {
      double X = xmin + (xmax - xmin) * iX / RAND_MAX;
      double Y = Fun(X);
      YMin = Y < YMin ? Y : YMin;
      YMax = Y > YMax ? Y : YMax;
    }
  }

  // Gets random values for X & Y
  double X = xmin + (xmax - xmin) * rand() / RAND_MAX;
  double Y = YMin + (YMax - YMin) * rand() / RAND_MAX;

  // Returns if valid and try again if not valid
  return Y < fun(X) ? X : Random(Fun, xmin, xmax);
}

I have added codes for initialising random number generator and to calculate the maximum and minimum values of the function.