Line | |
---|
1 | #include "madx.h" |
---|
2 | |
---|
3 | // private functions |
---|
4 | |
---|
5 | static void |
---|
6 | irngen(void) |
---|
7 | /* creates random number for frndm() */ |
---|
8 | { |
---|
9 | int i, j; |
---|
10 | for (i = 0; i < NJ_RAND; i++) |
---|
11 | { |
---|
12 | if ((j = irn_rand[i] - irn_rand[i+NR_RAND-NJ_RAND]) < 0) j += MAX_RAND; |
---|
13 | irn_rand[i] = j; |
---|
14 | } |
---|
15 | for (i = NJ_RAND; i < NR_RAND; i++) |
---|
16 | { |
---|
17 | if ((j = irn_rand[i] - irn_rand[i-NJ_RAND]) < 0) j += MAX_RAND; |
---|
18 | irn_rand[i] = j; |
---|
19 | } |
---|
20 | next_rand = 0; |
---|
21 | } |
---|
22 | |
---|
23 | // public interface |
---|
24 | |
---|
25 | void |
---|
26 | init55(int seed) |
---|
27 | /* initializes random number algorithm */ |
---|
28 | { |
---|
29 | int i, ii, k = 1, j = abs(seed)%MAX_RAND; |
---|
30 | irn_rand[NR_RAND-1] = j; |
---|
31 | for (i = 0; i < NR_RAND-1; i++) |
---|
32 | { |
---|
33 | ii = (ND_RAND*(i+1))%NR_RAND; |
---|
34 | irn_rand[ii-1] = k; |
---|
35 | if ((k = j - k) < 0) k += MAX_RAND; |
---|
36 | j = irn_rand[ii-1]; |
---|
37 | } |
---|
38 | /* warm up */ |
---|
39 | for (i = 0; i < 3; i++) irngen(); |
---|
40 | } |
---|
41 | |
---|
42 | double |
---|
43 | frndm(void) |
---|
44 | /* returns random number r with 0 <= r < 1 from flat distribution */ |
---|
45 | { |
---|
46 | const double one = 1; |
---|
47 | double scale = one / MAX_RAND; |
---|
48 | if (next_rand == NR_RAND) irngen(); |
---|
49 | return scale*irn_rand[next_rand++]; |
---|
50 | } |
---|
51 | |
---|
52 | double |
---|
53 | grndm(void) |
---|
54 | /* returns random number x from normal distribution */ |
---|
55 | { |
---|
56 | double xi1 = 2*frndm()-one, xi2=2*frndm()-one, zzr; |
---|
57 | while ((zzr = xi1*xi1+xi2*xi2) > one) |
---|
58 | { |
---|
59 | xi1 = 2*frndm()-one; xi2=2*frndm()-one; |
---|
60 | } |
---|
61 | zzr = sqrt(-2*log(zzr)/zzr); |
---|
62 | return xi1*zzr; |
---|
63 | } |
---|
64 | |
---|
65 | double |
---|
66 | tgrndm(double cut) |
---|
67 | /* returns random variable from normal distribution cat at 'cut' sigmas */ |
---|
68 | { |
---|
69 | double ret = zero; |
---|
70 | if (cut > zero) |
---|
71 | { |
---|
72 | ret = grndm(); |
---|
73 | while (fabs(ret) > fabs(cut)) ret = grndm(); |
---|
74 | } |
---|
75 | return ret; |
---|
76 | } |
---|
77 | |
---|
78 | |
---|
Note: See
TracBrowser
for help on using the repository browser.