blob: c49c08202191e3349015d229dbe527af6c2ec24c [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
Paul Jakma41b36e92006-12-08 01:09:50 +000013/*
14 * Prototypes.
15 */
16unsigned long timer(void);
17void init_rand(long);
18double rand01(void);
19double randg01(void);
20long nrand(long);
21void free_arc(void *);
22
23unsigned long timer ()
jardineb5d44e2003-12-23 08:09:43 +000024 { struct tms hold;
25
26 times(&hold);
hasso29e50b22005-09-01 18:18:47 +000027 return (unsigned long) ((float) (hold.tms_utime) / 60.0);
jardineb5d44e2003-12-23 08:09:43 +000028 }
29
30
31/*********************************************************************/
32/* */
33/* Family of random number generators */
34/* */
35/* Initialisation: */
36/* void init_rand ( seed ); */
37/* long seed - any positive number */
38/* if seed<=0 init_rand takes time */
39/* from timer instead of seed */
40/* */
41/* Whole number uniformly distributed on [0,n): */
42/* long nrand (n); */
43/* long n */
44/* */
45/* Real number uniformly distributed on [0,1] */
46/* double rand01(); */
47/* */
48/* Real number with Gauss(0,1) disitribution: */
49/* double randg01(); */
50/* */
51/* Algorithm: */
52/* x(n+1) = (x(n) * 5^13) mod 2^31 */
53/* */
54/*********************************************************************/
55
56unsigned long internal_seed;
57
58void init_rand ( init_seed )
59
60long init_seed;
61
62{ internal_seed = ( init_seed > 0 )
63 ? (unsigned long) init_seed
64 : (unsigned long) timer();
65
66
67 /* only odd numbers are acceptable */
68 if ( internal_seed % 2 == 0 ) internal_seed --;
69}
70
71/*********************************************************************/
72/* */
73/* Internal function irand may depend on OS and compiler */
74/* */
75/* irand assumption: */
76/* unsigned long i,j; */
77/* if i*j > max(unsigned long) */
78/* 1. No overflow interruption */
79/* 2. i*j = i*j mod max(unsigned long) */
80/* */
81/* This assumption is true for a lot of computers. */
82/* If your computer fails: */
83/* rename: irand <---> xrand */
84/* */
85/*********************************************************************/
86
87#define A 1220703125
88#define B 2147483647
89#define BF 2147483647.
90
91static long irand ()
92
93{ internal_seed = ( internal_seed * A ) & B;
94 return (long) internal_seed ;
95}
96
hasso29e50b22005-09-01 18:18:47 +000097#if 0 /* Not used. */
jardineb5d44e2003-12-23 08:09:43 +000098/*********************************************************************/
99/* */
100/* computer independent variant of irand */
101/* */
102/*********************************************************************/
103
104
105#define T15 32768
106#define T16 65536
107#define A1 37252
108#define A2 29589
109
110static long xrand()
111
112{ unsigned long is1, is2;
113
114 is1 = internal_seed / T15;
115 is2 = internal_seed % T15;
116
117 internal_seed = ( (((is2 * A1) + (is1 * A2))% T16 )* T15 + (is2 * A2) ) & B;
118 return (long) ( internal_seed ) ;
119}
hasso29e50b22005-09-01 18:18:47 +0000120#endif
jardineb5d44e2003-12-23 08:09:43 +0000121
122/*********************************************************************/
123
124
125double rand01()
126
hasso29e50b22005-09-01 18:18:47 +0000127{ return (double) (irand() / BF) ;
jardineb5d44e2003-12-23 08:09:43 +0000128}
129
130/*********************************************************************/
131
132#define NK 12
133
134double randg01()
135
136{ int i;
137 double sum = 0;
138
139 for ( i = 0; i < NK; i++ ) sum += rand01();
140 return sum - 6.;
141
142 /* if NK != 12 then you must return (12/NK)*sum - (NK/2) */
143}
144
145#undef NK
146
147
148/*********************************************************************/
149
150long nrand ( n )
151
152long n;
153
154{ return (long) ( rand01() * (double) n );
155}
156
157/*********************************************************************/
158
159#undef A
160#undef A1
161#undef A2
162#undef B
163#undef BF
164#undef T15
165#undef T16