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