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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
|
/*
* MLRAND.C - random number generation routines
* - most of the routines were written by Jim Rathkopf, LLNL
* - with modifications for portability and efficiency by SAB
*
* Source Version: 2.0
* Software Release #92-0043
*
*/
#include "cpyright.h"
#include "pml.h"
#ifndef RAND_MAX
#define RAND_MAX 2147483647.0
#endif
#define BASE24 16777216.0 /* 2^24 */
#define TWO8 256.0 /* 2^8 */
#define TWO16 65536.0 /* 2^16 */
#define TWO48 281474976710656.0 /* 2^48 */
#define UPPERPART(r) floor(r/BASE24)
#define LOWERPART(r, rupper) (r - rupper*BASE24)
typedef union u_rand_seed rand_seed;
union u_rand_seed
{double sd;
unsigned short sa[4];};
static double
dseed[2] = {16555217.0, 9732691.0},
dmult[2] = {15184245.0, 2651554.0},
dadder = 0.0;
static unsigned short
multplr[] = {0xb175, 0xa2e7, 0x2875};
static void
SC_DECLARE(_PM_seed_48, (unsigned short rs[3], unsigned short aso[3])),
SC_DECLARE(_PM_rand48_16to24, (unsigned short x16[3], double x24[2])),
SC_DECLARE(_PM_rand48_24to16, (double x24[2], unsigned short x16[3])),
SC_DECLARE(_PM_rnset, (unsigned short s[3], unsigned short mlt[3]));
static double
SC_DECLARE(_PM_random_48, (void));
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/* _PM_RANDOM_48 - generate pseudo-random numbers through a 48-bit linear
* - congruential algorithm
* - this emulates the drand48 library
* - of random number generators available on many, but not all,
* - UNIX machines
* - return a double x such that 0.0 <= x < 1.0
*
* algorithm:
* x(n+1) = (a*x(n) + c)mod m
*
* where
*
* defaults, for standard unix rand48, Cray defaults used in coding
*
* double name hexdecimal decimal
* x - seed-dseed[1], [2] 1234abcd330e 20017429951246
* a - multiplier-dmult[1], [2] 5deece66d 25214903917
* c - adder-dadder b 11
* m - 2**481000000000000 281474976710656
*
* 24-bit defaults (decimal) (lower bits listed first)
* x - dseed[1], [2]- 13447950.0, 1193131.0
* a - dmult[1], [2]- 15525485.0, 1502.0
* c - dadder- 11.0
*
*/
static double _PM_random_48()
{double t1_48, t2_48, t1_24[2], t2_24[2];
double d;
/* perform 48-bit arithmetic using two part data */
t1_48 = dmult[0]*dseed[0] + dadder;
t2_48 = dmult[1]*dseed[0] + dmult[0]*dseed[1];
t1_24[1] = UPPERPART(t1_48);
t1_24[0] = LOWERPART(t1_48, t1_24[1]);
t2_24[1] = UPPERPART(t2_48);
t2_24[0] = LOWERPART(t2_48, t2_24[1]);
t2_48 = t2_24[0] + t1_24[1];
t2_24[1] = UPPERPART(t2_48);
t2_24[0] = LOWERPART(t2_48, t2_24[1]);
d = (dseed[0] + dseed[1]*BASE24)/TWO48;
dseed[0] = t1_24[0];
dseed[1] = t2_24[0];
return(d);}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/* _PM_SEED_48 - convert a 3 by 16 representation of a seed
* - return a pointer to old seed
* - RS, new seed
*
*/
static void _PM_seed_48(rs, aso)
unsigned short rs[3], aso[3];
{double sdo[2];
sdo[0] = dseed[0];
sdo[1] = dseed[1];
_PM_rand48_24to16(sdo, aso);
_PM_rand48_16to24(rs, dseed);
return;}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/* _PM_RAND48_16TO24 - take a number stored in 3 16-bit shorts and move to
* - 2 doubles each containing 24 bits of data
* - X16, the 3 16-bit shorts
* - X24, the 2 doubles of 24 bits returned
*
*/
static void _PM_rand48_16to24(x16, x24)
unsigned short x16[3];
double x24[2];
{double t0, t1, t2, t1u, t1l;
t0 = (double) x16[0];
t1 = (double) x16[1];
t2 = (double) x16[2];
t1u = floor(t1/TWO8);
t1l = t1 - t1u*TWO8;
x24[0] = t0 + t1l*TWO16;
x24[1] = t1u + t2*TWO8;
return;}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/* _PM_RAND48_24TO16 - take a number stored in 2 doubles each containing
* - 24 bits of data move to 3 16-bit shorts
* - X24, the 2 doubles of 24 bits
* - X16, the 3 16-bit shorts returned
*/
static void _PM_rand48_24to16(x24, x16)
double x24[2];
unsigned short x16[3];
{double t0u, t0l, t1u, t1l;
t0u = floor(x24[0]/TWO16);
t0l = x24[0] - t0u*TWO16;
t1u = floor(x24[1]/TWO8);
t1l = x24[1] - t1u*TWO8;
x16[0] = t0l;
x16[1] = t0u + t1l*TWO8;
x16[2] = t1u;
return;}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/* _PM_RNSET - set random number generator seed, multiplier
* - S, the random number seed
* - MLT, the random number multiplier
*/
static void _PM_rnset(s, mlt)
unsigned short s[3];
unsigned short mlt[3];
{
/* construct 24-bit versions of 'seed' and 'mlt'
* and store in 64-bit doubles
*/
_PM_rand48_16to24(s, dseed);
_PM_rand48_16to24(mlt, dmult);
return;}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/* PM_SRAND_48 - change the value of the random number seed
* - ranset must be passed a 64 bit double
*/
void PM_srand_48(x)
double x;
{unsigned short wrk_seed[3];
int j;
rand_seed seed;
if (x == 0.0)
{wrk_seed[0] = 0x9cd1;
wrk_seed[1] = 0x53fc;
wrk_seed[2] = 0x9482;}
/* store x in wrk_seed */
else
{seed.sd = x;
for (j = 0; j < 3; j++)
wrk_seed[j] = seed.sa[j];};
_PM_rnset(wrk_seed, multplr);
return;}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/* PM_GRAND_48 - return the value of the current seed */
double PM_grand_48()
{int j;
unsigned short wrk_seed[3];
rand_seed seed;
static unsigned short SeedDum[3] = {0, 0, 0};
/* get current seed */
_PM_seed_48(SeedDum, wrk_seed);
/* set seed and multiplier */
_PM_rnset(wrk_seed, multplr);
/* store wrk_seed in seed */
for (j = 0; j < 3; j++)
seed.sa[j] = wrk_seed[j];
seed.sa[3] = 0x0;
return(seed.sd);}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/* PM_RANDOM_48 - generate a pseudo-random number
* - return a pseudo-random number between -1 and 1
* - use the linear congruential sequence function
*/
double PM_random_48(x)
double x;
{double v;
v = 2.0*_PM_random_48() - 1.0;
return(v);}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/* PM_RANDOM_S - return a random number between -1 and 1
* - use the system rand function
*/
double PM_random_s(x)
double x;
{double v;
static double rn = 0.0;
if (rn == 0.0)
rn = 2.0/RAND_MAX;
v = rn*rand() - 1.0;
return(v);}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/* PMRNDF - F77 interface to PM_random_48 */
REAL F77_ID(pmrndf_, pmrndf, PMRNDF)(px)
REAL *px;
{
return((REAL) PM_random_48((double) *px));}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/* PMRNDF - F77 interface to PM_random_s */
REAL F77_ID(pmrnds_, pmrnds, PMRNDS)(px)
REAL *px;
{
return((REAL) PM_random_s((double) *px));}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
|
