src/surf/random_mgr.c

   1 #include "surf/random_mgr.h"
   2 #include "xbt/sysdep.h"
   3
   4 #ifdef WIN32
   5
   6 static unsigned int _seed = 2147483647;
   7
   8 typedef unsigned __int64 uint64_t;
   9 typedef unsigned int uint32_t;
  10
  11 struct drand48_data
  12 {
  13         unsigned short int __x[3];  /* Current state.  */
  14         unsigned short int __old_x[3]; /* Old state.  */
  15         unsigned short int __c;     /* Additive const. in congruential formula.  */
  16         unsigned short int __init;  /* Flag for initializing.  */
  17         unsigned long long int __a; /* Factor in congruential formula.  */
  18 };
  19
  20 static struct drand48_data __libc_drand48_data = {0};
  21
  22 union ieee754_double
  23   {
  24         double d;
  25
  26         /* This is the IEEE 754 double-precision format.  */
  27         struct
  28         {
  29                 /* Together these comprise the mantissa.  */
  30                 unsigned int mantissa1:32;
  31                 unsigned int mantissa0:20;
  32                 unsigned int exponent:11;
  33                 unsigned int negative:1;
  34                 /* Little endian.  */
  35         } ieee;
  36
  37         /* This format makes it easier to see if a NaN is a signalling NaN.  */
  38         struct
  39         {
  40                 /* Together these comprise the mantissa.  */
  41                 unsigned int mantissa1:32;
  42                 unsigned int mantissa0:19;
  43                 unsigned int quiet_nan:1;
  44                 unsigned int exponent:11;
  45                 unsigned int negative:1;
  46
  47         } ieee_nan;
  48 };
  49
  50 #define IEEE754_DOUBLE_BIAS     0x3ff                                   /* Added to exponent.  */
  51
  52 double
  53 drand48 (void);
  54
  55 int
  56 _drand48_iterate (unsigned short int xsubi[3], struct drand48_data *buffer);
  57
  58 int
  59 _erand48_r (unsigned short int xsubi[3], struct drand48_data *buffer, double *result);
  60
  61
  62 int
  63 _erand48_r (unsigned short int xsubi[3], struct drand48_data *buffer, double *result)
  64 {
  65         union ieee754_double temp;
  66
  67         /* Compute next state.  */
  68         if (_drand48_iterate(xsubi, buffer) < 0)
  69                 return -1;
  70
  71         /* Construct a positive double with the 48 random bits distributed over
  72         its fractional part so the resulting FP number is [0.0,1.0).  */
  73
  74         temp.ieee.negative = 0;
  75         temp.ieee.exponent = IEEE754_DOUBLE_BIAS;
  76         temp.ieee.mantissa0 = (xsubi[2] << 4) | (xsubi[1] >> 12);
  77         temp.ieee.mantissa1 = ((xsubi[1] & 0xfff) << 20) | (xsubi[0] << 4);
  78
  79         /* Please note the lower 4 bits of mantissa1 are always 0.  */
  80         *result = temp.d - 1.0;
  81
  82         return 0;
  83 }
  84
  85 int
  86 _drand48_iterate (unsigned short int xsubi[3], struct drand48_data *buffer)
  87 {
  88         uint64_t X;
  89         uint64_t result;
  90
  91         /* Initialize buffer, if not yet done.  */
  92
  93         if(buffer->__init == 0)
  94         {
  95                 buffer->__a = 0x5deece66dull;
  96                 buffer->__c = 0xb;
  97                 buffer->__init = 1;
  98         }
  99
 100         /* Do the real work.  We choose a data type which contains at least
 101         48 bits.  Because we compute the modulus it does not care how
 102         many bits really are computed.  */
 103
 104         X = (uint64_t) xsubi[2] << 32 | (uint32_t) xsubi[1] << 16 | xsubi[0];
 105
 106         result = X * buffer->__a + buffer->__c;
 107
 108
 109         xsubi[0] = result & 0xffff;
 110         xsubi[1] = (result >> 16) & 0xffff;
 111         xsubi[2] = (result >> 32) & 0xffff;
 112
 113         return 0;
 114 }
 115
 116
 117 double
 118 _drand48 (void)
 119 {
 120         double result;
 121
 122         (void) _erand48_r (__libc_drand48_data.__x, &__libc_drand48_data, &result);
 123
 124          return result;
 125  }
 126
 127 void
 128 _srand(unsigned int seed)
 129 {
 130         _seed = seed;
 131 }
 132
 133 int
 134 _rand(void)
 135 {
 136         const long a = 16807;
 137         const long m = 2147483647;
 138         const long q = 127773; /* (m/a) */
 139         const long r = 2836; /* (m%a) */
 140
 141         long lo, k, s;
 142
 143         s = (long)_seed;
 144
 145         k = (long)(s/q);
 146
 147         lo = (s - q * k);
 148
 149         s = a * lo -r * k;
 150
 151         if(s <= 0)
 152                 s += m;
 153
 154         _seed = (int)(s & RAND_MAX);
 155
 156         return _seed;
 157 }
 158
 159 int
 160 _rand_r(unsigned int* pseed)
 161 {
 162         const long a = 16807;
 163         const long m = 2147483647;
 164         const long q = 127773;                  /* (m/a) */
 165         const long r = 2836;                    /* (m%a) */
 166
 167         long lo, k, s;
 168
 169         s = (long)*pseed;
 170
 171         k = (long)(s/q);
 172
 173         lo = (s - q * k);
 174
 175         s = a * lo -r * k;
 176
 177         if(s <= 0)
 178                 s += m;
 179
 180         return (int)(s & RAND_MAX);
 181
 182 }
 183
 184
 185 #define rand_r _rand_r
 186 #define drand48 _drand48
 187
 188 #endif
 189
 190 static double custom_random(Generator generator, long int *seed){
 191    switch(generator) {
 192
 193     case DRAND48:
 194       return drand48();
 195     case RAND:
 196       return (double)rand_r((unsigned int*)seed)/RAND_MAX;
 197     default:
 198       return drand48();
 199    }
 200 }
 201
 202 /* Generate numbers between min and max with a given mean and standard deviation */
 203 double random_generate(random_data_t random) {
 204   double a, b;
 205   double alpha, beta, gamma;
 206   double U1, U2, V, W, X;
 207
 208   if (random == NULL) return 0.0f;
 209
 210   if (random->std == 0)
 211      return random->mean * (random->max - random->min) + random->min;
 212
 213   a = random->mean * ( random->mean * (1 - random->mean) / (random->std*random->std) - 1 );
 214   b = (1 - random->mean) * ( random->mean * (1 - random->mean) / (random->std*random->std) - 1 );
 215
 216   alpha = a + b;
 217   if (a <= 1. || b <= 1.)
 218     beta = ((1./a)>(1./b))?(1./a):(1./b);
 219   else
 220     beta = sqrt ((alpha-2.) / (2.*a*b - alpha));
 221   gamma = a + 1./beta;
 222
 223   do {
 224     /* Random generation for the Beta distribution based on
 225      *   R. C. H. Cheng (1978). Generating beta variates with nonintegral shape parameters. _Communications of the ACM_, *21*, 317-322.
 226      *   It is good for speed because it does not call math functions many times and respect the 4 given constraints
 227      */
 228     U1 = custom_random(random->generator,&(random->seed));
 229     U2 = custom_random(random->generator,&(random->seed));
 230
 231     V = beta * log(U1/(1-U1));
 232     W = a * exp(V);
 233   } while (alpha * log(alpha/(b + W)) + gamma*V - log(4) < log(U1*U1*U2));
 234
 235   X = W / (b + W);
 236
 237   return X * (random->max - random->min) + random->min;
 238 }
 239
 240 random_data_t random_new(Generator generator, long int seed,
 241                          double min, double max,
 242                          double mean, double std){
 243   random_data_t random = xbt_new0(s_random_data_t, 1);
 244
 245   random->generator = generator;
 246   random->seed = seed;
 247   random->min = min;
 248   random->max = max;
 249
 250   /* Check user stupidities */
 251   if (max < min)
 252      THROW2(arg_error,0,"random->max < random->min (%f < %f)",max, min);
 253   if (mean < min)
 254      THROW2(arg_error,0,"random->mean < random->min (%f < %f)",mean, min);
 255   if (mean > max)
 256      THROW2(arg_error,0,"random->mean > random->max (%f > %f)",mean, max);
 257
 258   /* normalize the mean and standard deviation before storing */
 259   random->mean = (mean - min) / (max - min);
 260   random->std = std / (max - min);
 261
 262   if (random->mean * (1-random->mean) < random->std*random->std)
 263      THROW2(arg_error,0,"Invalid mean and standard deviation (%f and %f)",random->mean, random->std);
 264
 265   return random;
 266 }