1 /*************************************************************************/
2 /*************************************************************************/
3 /* Parallel 64-bit Linear Congruential Generator */
5 /* Author: Ashok Srinivasan, */
6 /* NCSA, University of Illinois, Urbana-Champaign */
7 /* E-Mail: ashoks@ncsa.uiuc.edu */
9 /* Note: The modulus is 2^64 */
11 /* Disclaimer: NCSA expressly disclaims any and all warranties, expressed*/
12 /* or implied, concerning the enclosed software. The intent in sharing */
13 /* this software is to promote the productive interchange of ideas */
14 /* throughout the research community. All software is furnished on an */
15 /* "as is" basis. No further updates to this software should be */
16 /* expected. Although this may occur, no commitment exists. The authors */
17 /* certainly invite your comments as well as the reporting of any bugs. */
18 /* NCSA cannot commit that any or all bugs will be fixed. */
19 /*************************************************************************/
20 /*************************************************************************/
22 /* This is version 1.0, created 20 May 1998 */
28 #define MAX_STREAMS 15613
30 int MAX_STREAMS = (146138719);
33 extern unsigned int _prime_list[];
35 unsigned int _PARAMLIST[3][2] = {{0x87b0b0fdU, 0x27bb2ee6U},
36 {0xe78b6955U,0x2c6fe96eU},
37 {0x31a53f85U,0x369dea0fU}};
39 /*************************************************************************/
40 /* You should not need to look at the next few lines! */
42 #define INIT_SEED1 0x2bc6ffffU
43 #define INIT_SEED0 0x8cfe166dU
45 #define TWO_M22 2.384185791015625e-07 /* 2^(-22) */
46 #define TWO_P22 4194304 /* 2^(22) */
47 #define TWO_M20 9.5367431640625e-07 /* 2^(-20) */
48 #define TWO_P20 1048576 /* 2^(20) */
49 #define TWO_M42 2.273736754432321e-13 /* 2^(-42) */
50 #define TWO_M64 5.4210108624275222e-20 /* 2^(-64) */
53 /************************************************************************/
55 CpvStaticDeclare(int, nstreams);
56 CpvStaticDeclare(CrnStream, _defaultStream);
60 CpvInitialize(int, nstreams);
61 CpvAccess(nstreams) = 0;
62 CpvInitialize(CrnStream, _defaultStream);
63 CrnInitStream(&CpvAccess(_defaultStream), 0, 0);
66 /* Initialize random number stream */
68 void CrnInitStream(CrnStream *genptr, unsigned int seed, int type)
70 unsigned int gennum = seed+CpvAccess(nstreams)*CmiNumPes();
73 genptr->prime = _prime_list[gennum%MAX_STREAMS];
74 genptr->multiplier[0] = (double) (_PARAMLIST[type][0]&0x3fffff);
75 genptr->multiplier[1] = (double)
76 (_PARAMLIST[type][0]>>22 | (_PARAMLIST[type][1]&0xfff)<<10);
77 genptr->multiplier[2] = (double) (_PARAMLIST[type][1]>>12);
78 genptr->state[0] = (double) ((INIT_SEED0^gennum)&0x3fffff);
79 genptr->state[1] = (double)
80 ((INIT_SEED0^gennum)>>22 | ((INIT_SEED1 ^ (unsigned)seed<<1)&0xfff)<<10);
81 genptr->state[2] = (double) ((INIT_SEED1 ^ (unsigned)seed<<1)>>12);
85 CpvAccess(nstreams)++;
89 #define advance_state(genptr) {double t0, t1, t2, t3, st0, st1, st2;\
90 t0 = genptr->state[0]*genptr->multiplier[0] + genptr->prime;\
91 t1 = (double) (int) (t0*TWO_M22); \
92 st0 = t0 - TWO_P22*t1; \
93 t1 += genptr->state[1]*genptr->multiplier[0] + \
94 genptr->state[0]*genptr->multiplier[1]; \
95 t2 = (double) (int) (t1*TWO_M22); \
96 st1 = t1 - TWO_P22*t2; \
97 t2 += genptr->state[2]*genptr->multiplier[0] + \
98 genptr->state[1]*genptr->multiplier[1] + \
99 genptr->state[0]*genptr->multiplier[2];\
100 t3 = (double) (int) (t2*TWO_M20); \
101 st2 = t2 - TWO_P20*t3; \
102 genptr->state[0] = st0; \
103 genptr->state[1] = st1; \
104 genptr->state[2] = st2;}
106 double CrnDouble(CrnStream *genptr)
108 advance_state(genptr); /* next state in sequence */
109 return genptr->state[2]*TWO_M20 + genptr->state[1]*TWO_M42 +
110 genptr->state[0]*TWO_M64;
113 int CrnInt(CrnStream *genptr)
115 return (int) (CrnDouble(genptr)*0x80000000U);
118 float CrnFloat(CrnStream *genptr)
120 return (float) CrnDouble(genptr);
123 void CrnSrand(unsigned int seed)
125 CrnInitStream(&CpvAccess(_defaultStream), seed, 0);
130 return CrnInt(&CpvAccess(_defaultStream));
133 double CrnDrand(void)
135 return CrnDouble(&CpvAccess(_defaultStream));
138 unsigned int _prime_list[MAX_STREAMS] =