src/surf/random_mgr.c

   1 #include "surf/random_mgr.h"
   2 #include "xbt/sysdep.h"
   3
   4 #ifdef WIN32
   5 static double drand48(void)
   6 {
   7    THROW_UNIMPLEMENTED;
   8    return -1;
   9 }
  10
  11 static double rand_r(unsigned int* seed)
  12 {
  13    THROW_UNIMPLEMENTED;
  14    return -1;
  15 }
  16 #endif
  17
  18 static double custom_random(Generator generator, long int *seed){
  19    switch(generator) {
  20
  21     case DRAND48:
  22       return drand48();
  23     case RAND:
  24       return (double)rand_r((unsigned int*)seed)/RAND_MAX;
  25     default:
  26       return drand48();
  27    }
  28 }
  29
  30 /* Generate numbers between min and max with a given mean and standard deviation */
  31 double random_generate(random_data_t random) {
  32   double a, b;
  33   double alpha, beta, gamma;
  34   double U1, U2, V, W, X;
  35
  36   if (random == NULL) return 0.0f;
  37
  38   a = random->mean * ( random->mean * (1 - random->mean) / (random->std*random->std) - 1 );
  39   b = (1 - random->mean) * ( random->mean * (1 - random->mean) / (random->std*random->std) - 1 );
  40
  41   alpha = a + b;
  42   if (a <= 1. || b <= 1.)
  43     beta = ((1./a)>(1./b))?(1./a):(1./b);
  44   else
  45     beta = sqrt ((alpha-2.) / (2.*a*b - alpha));
  46   gamma = a + 1./beta;
  47
  48   do {
  49     /* Random generation for the Beta distribution based on
  50      *   R. C. H. Cheng (1978). Generating beta variates with nonintegral shape parameters. _Communications of the ACM_, *21*, 317-322.
  51      *   It is good for speed because it does not call math functions many times and respect the 4 given constraints
  52      */
  53     U1 = custom_random(random->generator,&(random->seed));
  54     U2 = custom_random(random->generator,&(random->seed));
  55
  56     V = beta * log(U1/(1-U1));
  57     W = a * exp(V);
  58   } while (alpha * log(alpha/(b + W)) + gamma*V - log(4) < log(U1*U1*U2));
  59
  60   X = W / (b + W);
  61
  62   return X * (random->max - random->min) + random->min;
  63 }
  64
  65 random_data_t random_new(Generator generator, long int seed,
  66                          double min, double max,
  67                          double mean, double std){
  68   random_data_t random = xbt_new0(s_random_data_t, 1);
  69
  70   random->generator = generator;
  71   random->seed = seed;
  72   random->min = min;
  73   random->max = max;
  74
  75   /* Check user stupidities */
  76   if (max < min)
  77      THROW2(arg_error,0,"random->max < random->min (%f < %f)",max, min);
  78   if (mean < min)
  79      THROW2(arg_error,0,"random->mean < random->min (%f < %f)",mean, min);
  80   if (mean > max)
  81      THROW2(arg_error,0,"random->mean > random->max (%f > %f)",mean, max);
  82
  83   /* normalize the mean and standard deviation before storing */
  84   random->mean = (mean - min) / (max - min);
  85   random->std = std / (max - min);
  86
  87   if (random->mean * (1-random->mean) < random->std*random->std)
  88      THROW2(arg_error,0,"Invalid mean and standard deviation (%f and %f)",random->mean, random->std);
  89
  90   return random;
  91 }