Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
simpleRandom.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file simpleRandom.cc
1 #include "simpleRandom.h"
2 #include <math.h>
3 #if defined(SunOS) || defined(Linux) || defined(OSF1)
4 #include <strings.h>
5 #else
6 #include <bstring.h>
7 #endif
8 
9 
10 /*
11  * Shuffle the bytes into little-endian order within words, as per the
12  * MD5 spec. Note: this code works regardless of the byte order.
13  */
14 void
15 simpleRandom::byteSwap(word32 *buf, unsigned words)
16 {
17  xbyte *p = (xbyte *)buf;
18 
19  do {
20  *buf++ = (word32)((unsigned)p[3] << 8 | p[2]) << 16 |
21  ((unsigned)p[1] << 8 | p[0]);
22  p += 4;
23  } while (--words);
24 }
25 
26 /*
27  * Start MD5 accumulation. Set bit count to 0 and buffer to mysterious
28  * initialization constants.
29  */
30 void
31 simpleRandom::xMD5Init(struct xMD5Context *ctx)
32 {
33  ctx->buf[0] = 0x67452301;
34  ctx->buf[1] = 0xefcdab89;
35  ctx->buf[2] = 0x98badcfe;
36  ctx->buf[3] = 0x10325476;
37 
38  ctx->bytes[0] = 0;
39  ctx->bytes[1] = 0;
40 }
41 
42 /*
43  * Update context to reflect the concatenation of another buffer full
44  * of bytes.
45  */
46 void
47 simpleRandom::xMD5Update(struct xMD5Context *ctx, xbyte const *buf, unsigned int len)
48 {
49  word32 t;
50 
51  /* Update byte count */
52 
53  t = ctx->bytes[0];
54  if ((ctx->bytes[0] = t + len) < t)
55  ctx->bytes[1]++; /* Carry from low to high */
56 
57  t = 64 - (t & 0x3f); /* Space available in ctx->in (at least 1) */
58  if ( t > len) {
59  bcopy(buf, (xbyte *)ctx->in + 64 - (unsigned)t, len);
60  return;
61  }
62  /* First chunk is an odd size */
63  bcopy(buf,(xbyte *)ctx->in + 64 - (unsigned)t, (unsigned)t);
64  byteSwap(ctx->in, 16);
65  xMD5Transform(ctx->buf, ctx->in);
66  buf += (unsigned)t;
67  len -= (unsigned)t;
68 
69  /* Process data in 64-byte chunks */
70  while (len >= 64) {
71  bcopy(buf, ctx->in, 64);
72  byteSwap(ctx->in, 16);
73  xMD5Transform(ctx->buf, ctx->in);
74  buf += 64;
75  len -= 64;
76  }
77 
78  /* Handle any remaining bytes of data. */
79  bcopy(buf, ctx->in, len);
80 }
81 
82 /*
83  * Final wrapup - pad to 64-byte boundary with the bit pattern
84  * 1 0* (64-bit count of bits processed, MSB-first)
85  */
86 void
87 simpleRandom::xMD5Final(xbyte bdigest[16], struct xMD5Context *ctx)
88 {
89  int count = (int)(ctx->bytes[0] & 0x3f); /* Bytes in ctx->in */
90  xbyte *p = (xbyte *)ctx->in + count; /* First unused byte */
91 
92  /* Set the first char of padding to 0x80. There is always room. */
93  *p++ = 0x80;
94 
95  /* Bytes of padding needed to make 56 bytes (-8..55) */
96  count = 56 - 1 - count;
97 
98  if (count < 0) { /* Padding forces an extra block */
99  bzero(p, count+8);
100  byteSwap(ctx->in, 16);
101  xMD5Transform(ctx->buf, ctx->in);
102  p = (xbyte *)ctx->in;
103  count = 56;
104  }
105  bzero(p, count+8);
106  byteSwap(ctx->in, 14);
107 
108  /* Append length in bits and transform */
109  ctx->in[14] = ctx->bytes[0] << 3;
110  ctx->in[15] = ctx->bytes[1] << 3 | ctx->bytes[0] >> 29;
111  xMD5Transform(ctx->buf, ctx->in);
112 
113  byteSwap(ctx->buf, 4);
114  bcopy(ctx->buf, bdigest, 16);
115  bzero(ctx,sizeof(*ctx));
116 }
117 
118 
119 /* The four core functions - F1 is optimized somewhat */
120 
121 /* #define F1(x, y, z) (x & y | ~x & z) */
122 #define F1(x, y, z) (z ^ (x & (y ^ z)))
123 #define F2(x, y, z) F1(z, x, y)
124 #define F3(x, y, z) (x ^ y ^ z)
125 #define F4(x, y, z) (y ^ (x | ~z))
126 
127 /* This is the central step in the MD5 algorithm. */
128 #define MD5STEP(f,w,x,y,z,in,s) \
129  (w += f(x,y,z) + in, w = (w<<s | w>>(32-s)) + x)
130 
131 /*
132  * The core of the MD5 algorithm, this alters an existing MD5 hash to
133  * reflect the addition of 16 longwords of new data. MD5Update blocks
134  * the data and converts bytes into longwords for this routine.
135  */
136 void
137 simpleRandom::xMD5Transform(word32 buf[4], word32 const in[16])
138 {
139  word32 a, b, c, d;
140 
141  a = buf[0];
142  b = buf[1];
143  c = buf[2];
144  d = buf[3];
145 
146  MD5STEP(F1, a, b, c, d, in[0] + 0xd76aa478, 7);
147  MD5STEP(F1, d, a, b, c, in[1] + 0xe8c7b756, 12);
148  MD5STEP(F1, c, d, a, b, in[2] + 0x242070db, 17);
149  MD5STEP(F1, b, c, d, a, in[3] + 0xc1bdceee, 22);
150  MD5STEP(F1, a, b, c, d, in[4] + 0xf57c0faf, 7);
151  MD5STEP(F1, d, a, b, c, in[5] + 0x4787c62a, 12);
152  MD5STEP(F1, c, d, a, b, in[6] + 0xa8304613, 17);
153  MD5STEP(F1, b, c, d, a, in[7] + 0xfd469501, 22);
154  MD5STEP(F1, a, b, c, d, in[8] + 0x698098d8, 7);
155  MD5STEP(F1, d, a, b, c, in[9] + 0x8b44f7af, 12);
156  MD5STEP(F1, c, d, a, b, in[10] + 0xffff5bb1, 17);
157  MD5STEP(F1, b, c, d, a, in[11] + 0x895cd7be, 22);
158  MD5STEP(F1, a, b, c, d, in[12] + 0x6b901122, 7);
159  MD5STEP(F1, d, a, b, c, in[13] + 0xfd987193, 12);
160  MD5STEP(F1, c, d, a, b, in[14] + 0xa679438e, 17);
161  MD5STEP(F1, b, c, d, a, in[15] + 0x49b40821, 22);
162 
163  MD5STEP(F2, a, b, c, d, in[1] + 0xf61e2562, 5);
164  MD5STEP(F2, d, a, b, c, in[6] + 0xc040b340, 9);
165  MD5STEP(F2, c, d, a, b, in[11] + 0x265e5a51, 14);
166  MD5STEP(F2, b, c, d, a, in[0] + 0xe9b6c7aa, 20);
167  MD5STEP(F2, a, b, c, d, in[5] + 0xd62f105d, 5);
168  MD5STEP(F2, d, a, b, c, in[10] + 0x02441453, 9);
169  MD5STEP(F2, c, d, a, b, in[15] + 0xd8a1e681, 14);
170  MD5STEP(F2, b, c, d, a, in[4] + 0xe7d3fbc8, 20);
171  MD5STEP(F2, a, b, c, d, in[9] + 0x21e1cde6, 5);
172  MD5STEP(F2, d, a, b, c, in[14] + 0xc33707d6, 9);
173  MD5STEP(F2, c, d, a, b, in[3] + 0xf4d50d87, 14);
174  MD5STEP(F2, b, c, d, a, in[8] + 0x455a14ed, 20);
175  MD5STEP(F2, a, b, c, d, in[13] + 0xa9e3e905, 5);
176  MD5STEP(F2, d, a, b, c, in[2] + 0xfcefa3f8, 9);
177  MD5STEP(F2, c, d, a, b, in[7] + 0x676f02d9, 14);
178  MD5STEP(F2, b, c, d, a, in[12] + 0x8d2a4c8a, 20);
179 
180  MD5STEP(F3, a, b, c, d, in[5] + 0xfffa3942, 4);
181  MD5STEP(F3, d, a, b, c, in[8] + 0x8771f681, 11);
182  MD5STEP(F3, c, d, a, b, in[11] + 0x6d9d6122, 16);
183  MD5STEP(F3, b, c, d, a, in[14] + 0xfde5380c, 23);
184  MD5STEP(F3, a, b, c, d, in[1] + 0xa4beea44, 4);
185  MD5STEP(F3, d, a, b, c, in[4] + 0x4bdecfa9, 11);
186  MD5STEP(F3, c, d, a, b, in[7] + 0xf6bb4b60, 16);
187  MD5STEP(F3, b, c, d, a, in[10] + 0xbebfbc70, 23);
188  MD5STEP(F3, a, b, c, d, in[13] + 0x289b7ec6, 4);
189  MD5STEP(F3, d, a, b, c, in[0] + 0xeaa127fa, 11);
190  MD5STEP(F3, c, d, a, b, in[3] + 0xd4ef3085, 16);
191  MD5STEP(F3, b, c, d, a, in[6] + 0x04881d05, 23);
192  MD5STEP(F3, a, b, c, d, in[9] + 0xd9d4d039, 4);
193  MD5STEP(F3, d, a, b, c, in[12] + 0xe6db99e5, 11);
194  MD5STEP(F3, c, d, a, b, in[15] + 0x1fa27cf8, 16);
195  MD5STEP(F3, b, c, d, a, in[2] + 0xc4ac5665, 23);
196 
197  MD5STEP(F4, a, b, c, d, in[0] + 0xf4292244, 6);
198  MD5STEP(F4, d, a, b, c, in[7] + 0x432aff97, 10);
199  MD5STEP(F4, c, d, a, b, in[14] + 0xab9423a7, 15);
200  MD5STEP(F4, b, c, d, a, in[5] + 0xfc93a039, 21);
201  MD5STEP(F4, a, b, c, d, in[12] + 0x655b59c3, 6);
202  MD5STEP(F4, d, a, b, c, in[3] + 0x8f0ccc92, 10);
203  MD5STEP(F4, c, d, a, b, in[10] + 0xffeff47d, 15);
204  MD5STEP(F4, b, c, d, a, in[1] + 0x85845dd1, 21);
205  MD5STEP(F4, a, b, c, d, in[8] + 0x6fa87e4f, 6);
206  MD5STEP(F4, d, a, b, c, in[15] + 0xfe2ce6e0, 10);
207  MD5STEP(F4, c, d, a, b, in[6] + 0xa3014314, 15);
208  MD5STEP(F4, b, c, d, a, in[13] + 0x4e0811a1, 21);
209  MD5STEP(F4, a, b, c, d, in[4] + 0xf7537e82, 6);
210  MD5STEP(F4, d, a, b, c, in[11] + 0xbd3af235, 10);
211  MD5STEP(F4, c, d, a, b, in[2] + 0x2ad7d2bb, 15);
212  MD5STEP(F4, b, c, d, a, in[9] + 0xeb86d391, 21);
213 
214  buf[0] += a;
215  buf[1] += b;
216  buf[2] += c;
217  buf[3] += d;
218 }
219 
220 
221 void simpleRandom::MD5(xbyte *dest, const xbyte *orig, unsigned int len)
222 {
223  struct xMD5Context context;
224 
225  xMD5Init(&context);
226  xMD5Update(&context, orig, len);
227  xMD5Final(dest, &context);
228 }
229 
230 /*****************************************************************/
231 
232 simpleRandom::simpleRandom()
233 {
234  digest[0] = 7657635;
235  digest[1] = 5565649;
236  digest[2] = 9827729;
237  digest[3] = 9892898;
238 
239 
240 }
241 
242 simpleRandom::simpleRandom( const int iseed)
243 {
244  digest[0] = iseed;
245  digest[1] = 565649;
246  digest[2] = 6827729;
247  digest[3] = 2892898;
248 
249 }
250 
251 float simpleRandom::gauss(const float mean, const float sigma)
252 {
253 #define BIGNUMBER 100000
254 
255  float y = 0.;
256  while(y == 0.)
257  {
258  y = rnd(0,BIGNUMBER);
259  }
260  float z = rnd(0,BIGNUMBER);
261 
262  y /= BIGNUMBER;
263  z /= BIGNUMBER;
264 
265  float x = z * 6.283185;
266 
267 #if defined(SunOS) || defined(Linux)
268  float result = mean + sigma* sin(x) * sqrt(-2 * log(y) );
269 #else
270  float result = mean + sigma* sinf(x) * sqrtf(-2 * logf(y) );
271 #endif
272 
273  return result;
274 }
275 
276 float simpleRandom::rnd(int low, int high)
277 {
278  unsigned int range = high - low;
279  unsigned int mask = 0;
280  unsigned int num;
281  int r;
282 
283  for (r = range; r; r >>= 1)
284  mask |= r;
285 
286  do {
287  MD5((xbyte*)digest,(const xbyte *) digest, sizeof(digest));
288  num = digest[0] & mask;
289  } while (num > range);
290 
291  return num + low;
292 }