blob: d4ef9950056a3cbf6b7747380c401cc430c0ec5c [file] [log] [blame]
jardineb5d44e2003-12-23 08:09:43 +00001/*********************************************************************/
2/* */
3/* current processor time in seconds */
4/* difference between two calls is processor time spent by your code */
5/* needs: <sys/types.h>, <sys/times.h> */
6/* depends on compiler and OS */
7/* */
8/*********************************************************************/
9
10#include <sys/types.h>
11#include <sys/times.h>
12
13float timer()
14 { struct tms hold;
15
16 times(&hold);
17 return (float)(hold.tms_utime) / 60.0;
18 }
19
20
21/*********************************************************************/
22/* */
23/* Family of random number generators */
24/* */
25/* Initialisation: */
26/* void init_rand ( seed ); */
27/* long seed - any positive number */
28/* if seed<=0 init_rand takes time */
29/* from timer instead of seed */
30/* */
31/* Whole number uniformly distributed on [0,n): */
32/* long nrand (n); */
33/* long n */
34/* */
35/* Real number uniformly distributed on [0,1] */
36/* double rand01(); */
37/* */
38/* Real number with Gauss(0,1) disitribution: */
39/* double randg01(); */
40/* */
41/* Algorithm: */
42/* x(n+1) = (x(n) * 5^13) mod 2^31 */
43/* */
44/*********************************************************************/
45
46unsigned long internal_seed;
47
48void init_rand ( init_seed )
49
50long init_seed;
51
52{ internal_seed = ( init_seed > 0 )
53 ? (unsigned long) init_seed
54 : (unsigned long) timer();
55
56
57 /* only odd numbers are acceptable */
58 if ( internal_seed % 2 == 0 ) internal_seed --;
59}
60
61/*********************************************************************/
62/* */
63/* Internal function irand may depend on OS and compiler */
64/* */
65/* irand assumption: */
66/* unsigned long i,j; */
67/* if i*j > max(unsigned long) */
68/* 1. No overflow interruption */
69/* 2. i*j = i*j mod max(unsigned long) */
70/* */
71/* This assumption is true for a lot of computers. */
72/* If your computer fails: */
73/* rename: irand <---> xrand */
74/* */
75/*********************************************************************/
76
77#define A 1220703125
78#define B 2147483647
79#define BF 2147483647.
80
81static long irand ()
82
83{ internal_seed = ( internal_seed * A ) & B;
84 return (long) internal_seed ;
85}
86
87/*********************************************************************/
88/* */
89/* computer independent variant of irand */
90/* */
91/*********************************************************************/
92
93
94#define T15 32768
95#define T16 65536
96#define A1 37252
97#define A2 29589
98
99static long xrand()
100
101{ unsigned long is1, is2;
102
103 is1 = internal_seed / T15;
104 is2 = internal_seed % T15;
105
106 internal_seed = ( (((is2 * A1) + (is1 * A2))% T16 )* T15 + (is2 * A2) ) & B;
107 return (long) ( internal_seed ) ;
108}
109
110
111/*********************************************************************/
112
113
114double rand01()
115
116{ return (double) irand() / BF ;
117}
118
119/*********************************************************************/
120
121#define NK 12
122
123double randg01()
124
125{ int i;
126 double sum = 0;
127
128 for ( i = 0; i < NK; i++ ) sum += rand01();
129 return sum - 6.;
130
131 /* if NK != 12 then you must return (12/NK)*sum - (NK/2) */
132}
133
134#undef NK
135
136
137/*********************************************************************/
138
139long nrand ( n )
140
141long n;
142
143{ return (long) ( rand01() * (double) n );
144}
145
146/*********************************************************************/
147
148#undef A
149#undef A1
150#undef A2
151#undef B
152#undef BF
153#undef T15
154#undef T16