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
|
/* Integer and floating point random number generator routines */
/*
* Copyright 2000 Graeme W. Gill
* All rights reserved.
*
* This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
* see the License.txt file for licencing details.
*/
#include "numsup.h"
#include "rand.h"
/* 32 bit pseudo random sequencer based on XOR feedback */
/* generates number between 1 and 4294967295 */
#define PSRAND32(S) (((S) & 0x80000000) ? (((S) << 1) ^ 0xa398655d) : ((S) << 1))
/* 32 bit linear congruent generator */
/* generates number between 0 and 4294967295 */
/* (From Knuth & H.W.Lewis) */
#define PSRAND32L(S) ((S) * 1664525L + 1013904223L)
/* Return a 32 bit number between 0 and 4294967295 */
/* Use Knuth shuffle to improve PSRAND32 sequence */
unsigned int
rand32( /* Return 32 bit random number */
unsigned int seed /* Optional seed. Non-zero re-initialized with that seed */
) {
#define TSIZE 2843 /* Prime */
static unsigned int last, ran = 0x12345678; /* Default seed */
static unsigned int pvs[TSIZE];
static int pvs_inited = 0;
int i;
if (seed != 0)
{
pvs_inited = 0;
ran = seed;;
}
/* Init random storage locations */
if (pvs_inited == 0)
{
for (i = 0; i < TSIZE; i++)
pvs[i] = ran = PSRAND32(ran);
last = ran;
pvs_inited = 1;
}
i = last % TSIZE; /* New location */
last = pvs[i]; /* Value generated */
pvs[i] = ran = PSRAND32(ran); /* New value */
return last-1;
}
/* return a random number between 0.0 and 1.0 */
/* based on rand32 */
double ranno(void) {
return rand32(0) / 4294967295.0;
}
/* Return a random double in the range min to max */
double
d_rand(double min, double max) {
double tt;
tt = ranno();
return min + (max - min) * tt;
}
/* Return a random integer in the range min to max inclusive */
int
i_rand(int min, int max) {
double tt;
tt = ranno();
return min + (int)floor(0.5 + ((double)(max - min)) * tt);
}
/* Return a random floating point number with a gausian/normal */
/* distribution, centered about 0.0, with standard deviation 1.0 */
/* This uses the Box-Muller transformation */
double norm_rand(void) {
static int r2 = 0; /* Can use 2nd number generated */
static double nr2;
if (r2 == 0) {
double v1, v2, t1, t2, r1;
do {
v1 = d_rand(-1.0, 1.0);
v2 = d_rand(-1.0, 1.0);
t1 = v1 * v1 + v2 * v2;
} while (t1 == 0.0 || t1 >= 1.0);
t2 = sqrt(-2.0 * log(t1)/t1);
nr2 = v2 * t2; /* One for next time */
r2 = 1;
r1 = v1 * t2;
return r1;
} else {
r2 = 0;
return nr2;
}
}
|
