1 /*
2 * Kuklomenos
3 * Copyright (C) 2008-2009 Martin Bays <mbays@sdf.lonestar.org>
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see http://www.gnu.org/licenses/.
17 */
18
19 #include <cmath>
20 #include <cstdlib>
21
ranf()22 float ranf()
23 {
24 return (float)(rand())/RAND_MAX;
25 }
ranf(float m)26 float ranf(float m)
27 {
28 return m*ranf();
29 }
30
31 // rani(n) returns uniformly random integer in [0,n).
rani(int m)32 int rani(int m)
33 {
34 float r;
35 do
36 { r = ranf(); } while (r == 1);
37
38 return (int)floorf(r*m);
39 }
40
41 /* Returns a random number with distribution N(0,1). Uses Box-Muller. Code
42 * taken from http://www.taygeta.com/random/gaussian.html.
43 */
gaussian()44 float gaussian()
45 {
46 float x1, x2, w, y;
47 static float spare_y = 1000;
48
49 // the algorithm below returns two independent normally distributed
50 // numbers, so we store the spare one in a static variable for use on the
51 // subsequent call.
52 if (spare_y != 1000)
53 {
54 float ret = spare_y;
55 spare_y = 1000;
56 return ret;
57 }
58
59 do {
60 x1 = 2.0 * ranf() - 1.0;
61 x2 = 2.0 * ranf() - 1.0;
62 w = x1 * x1 + x2 * x2;
63 } while ( w >= 1.0 );
64
65 w = sqrt( (-2.0 * log( w ) ) / w );
66 y = x1 * w;
67 spare_y = x2 * w;
68
69 return y;
70 }
71