Logo AND Algorithmique Numérique Distribuée

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