blob: d4ef9950056a3cbf6b7747380c401cc430c0ec5c (
plain)
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
|
/*********************************************************************/
/* */
/* current processor time in seconds */
/* difference between two calls is processor time spent by your code */
/* needs: <sys/types.h>, <sys/times.h> */
/* depends on compiler and OS */
/* */
/*********************************************************************/
#include <sys/types.h>
#include <sys/times.h>
float timer()
{ struct tms hold;
times(&hold);
return (float)(hold.tms_utime) / 60.0;
}
/*********************************************************************/
/* */
/* Family of random number generators */
/* */
/* Initialisation: */
/* void init_rand ( seed ); */
/* long seed - any positive number */
/* if seed<=0 init_rand takes time */
/* from timer instead of seed */
/* */
/* Whole number uniformly distributed on [0,n): */
/* long nrand (n); */
/* long n */
/* */
/* Real number uniformly distributed on [0,1] */
/* double rand01(); */
/* */
/* Real number with Gauss(0,1) disitribution: */
/* double randg01(); */
/* */
/* Algorithm: */
/* x(n+1) = (x(n) * 5^13) mod 2^31 */
/* */
/*********************************************************************/
unsigned long internal_seed;
void init_rand ( init_seed )
long init_seed;
{ internal_seed = ( init_seed > 0 )
? (unsigned long) init_seed
: (unsigned long) timer();
/* only odd numbers are acceptable */
if ( internal_seed % 2 == 0 ) internal_seed --;
}
/*********************************************************************/
/* */
/* Internal function irand may depend on OS and compiler */
/* */
/* irand assumption: */
/* unsigned long i,j; */
/* if i*j > max(unsigned long) */
/* 1. No overflow interruption */
/* 2. i*j = i*j mod max(unsigned long) */
/* */
/* This assumption is true for a lot of computers. */
/* If your computer fails: */
/* rename: irand <---> xrand */
/* */
/*********************************************************************/
#define A 1220703125
#define B 2147483647
#define BF 2147483647.
static long irand ()
{ internal_seed = ( internal_seed * A ) & B;
return (long) internal_seed ;
}
/*********************************************************************/
/* */
/* computer independent variant of irand */
/* */
/*********************************************************************/
#define T15 32768
#define T16 65536
#define A1 37252
#define A2 29589
static long xrand()
{ unsigned long is1, is2;
is1 = internal_seed / T15;
is2 = internal_seed % T15;
internal_seed = ( (((is2 * A1) + (is1 * A2))% T16 )* T15 + (is2 * A2) ) & B;
return (long) ( internal_seed ) ;
}
/*********************************************************************/
double rand01()
{ return (double) irand() / BF ;
}
/*********************************************************************/
#define NK 12
double randg01()
{ int i;
double sum = 0;
for ( i = 0; i < NK; i++ ) sum += rand01();
return sum - 6.;
/* if NK != 12 then you must return (12/NK)*sum - (NK/2) */
}
#undef NK
/*********************************************************************/
long nrand ( n )
long n;
{ return (long) ( rand01() * (double) n );
}
/*********************************************************************/
#undef A
#undef A1
#undef A2
#undef B
#undef BF
#undef T15
#undef T16
|
