1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*---------------------------------------------------------------------------*/
8 /*
9 Some utility routines for generating random numbers.
10
11 File: randgen.c
12 Author: B. Douglas Ward
13 Date: 15 February 1999
14
15 Mod: Correction to rand_normal
16 Date: 15 September 1999
17
18 Mod: Added routine rand_initialize.
19 Date: 28 January 2000
20 */
21
22
23 /*---------------------------------------------------------------------------*/
24 /*
25 Initialize the random number generator.
26 */
27
rand_initialize(long int seedval)28 void rand_initialize (long int seedval)
29 {
30 srand48 (seedval);
31 }
32
33
34 /*---------------------------------------------------------------------------*/
35 /*
36 Routine to generate a uniform U(a,b) random variate.
37 */
38
rand_uniform(float a,float b)39 float rand_uniform (float a, float b)
40 {
41 return (a + (float)drand48() * (b-a) );
42 }
43
44
45 /*---------------------------------------------------------------------------*/
46 /*
47 Routine to generate a normal N(mu,var) random variate.
48 */
49
rand_normal(float mu,float var)50 float rand_normal (float mu, float var)
51 {
52 float u1, u2;
53 float r, n;
54
55
56 u1 = 0.0;
57 while (u1 <= 0.0)
58 {
59 u1 = rand_uniform (0.0, 1.0);
60 }
61 u2 = rand_uniform (0.0, 1.0);
62
63 r = sqrt(-2.0*log(u1));
64 n = r * cos(2.0*PI*u2);
65
66 return (mu + n * sqrt(var));
67 }
68
69
70 /*---------------------------------------------------------------------------*/
71 /*
72 Routine to generate two independent normal N(mu,var) random variates.
73 */
74
rand_binormal(float mu,float var,float * n1,float * n2)75 void rand_binormal (float mu, float var, float * n1, float * n2)
76 {
77 float u1, u2;
78 float r, sigma;
79
80
81 u1 = 0.0;
82 while (u1 <= 0.0)
83 {
84 u1 = rand_uniform (0.0, 1.0);
85 }
86 u2 = rand_uniform (0.0, 1.0);
87
88 r = sqrt(-2.0*log(u1));
89 sigma = sqrt (var);
90
91 *n1 = mu + r * cos(2.0*PI*u2) * sigma;
92 *n2 = mu + r * sin(2.0*PI*u2) * sigma;
93 }
94
95
96 /*---------------------------------------------------------------------------*/
97
98
99
100