1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
|
/* r250.c the r250 uniform random number algorithm
Kirkpatrick, S., and E. Stoll, 1981; "A Very Fast
Shift-Register Sequence Random Number Generator",
Journal of Computational Physics, V.40
also:
see W.L. Maier, DDJ May 1991
*/
#include <limits.h>
#include "r250.h"
// static functions
static unsigned int randlcg(int);
#define MSB 0x80000000L
#define ALL_BITS 0xffffffffL
#define HALF_RANGE 0x40000000L
#define STEP 7
#define BITS 32
static unsigned int r250_buffer[ 250 ];
static int r250_index;
static unsigned int randlcg(int sd) /* returns a random unsigned integer */
{
static int quotient1 = LONG_MAX / 16807L;
static int remainder1 = LONG_MAX % 16807L;
if ( sd <= quotient1 )
sd = (sd * 16807L) % LONG_MAX;
else
{
int high_part = sd / quotient1;
int low_part = sd % quotient1;
int test = 16807L * low_part - remainder1 * high_part;
if ( test > 0 )
sd = test;
else
sd = test + LONG_MAX;
}
return sd;
}
void r250_init(int sd)
{
int j, k;
unsigned int mask, msb;
r250_index = 0;
for (j = 0; j < 250; j++) /* fill r250 buffer with BITS-1 bit values */
sd = r250_buffer[j] = randlcg(sd);
for (j = 0; j < 250; j++) /* set some MSBs to 1 */
if ( (sd=randlcg(sd)) > HALF_RANGE )
r250_buffer[j] |= MSB;
msb = MSB; /* turn on diagonal bit */
mask = ALL_BITS; /* turn off the leftmost bits */
for (j = 0; j < BITS; j++)
{
k = STEP * j + 3; /* select a word to operate on */
r250_buffer[k] &= mask; /* turn off bits left of the diagonal */
r250_buffer[k] |= msb; /* turn on the diagonal bit */
mask >>= 1;
msb >>= 1;
}
}
unsigned int r250() /* returns a random unsigned integer */
{
register int j;
register unsigned int new_rand;
if ( r250_index >= 147 )
j = r250_index - 147; /* wrap pointer around */
else
j = r250_index + 103;
new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
r250_buffer[ r250_index ] = new_rand;
if ( r250_index >= 249 ) /* increment pointer for next time */
r250_index = 0;
else
r250_index++;
return new_rand;
}
double dr250() /* returns a random double in range 0..1 */
{
register int j;
register unsigned int new_rand;
if ( r250_index >= 147 )
j = r250_index - 147; /* wrap pointer around */
else
j = r250_index + 103;
new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
r250_buffer[ r250_index ] = new_rand;
if ( r250_index >= 249 ) /* increment pointer for next time */
r250_index = 0;
else
r250_index++;
return (double)new_rand / ALL_BITS;
}
|
