1
2
3 package org.catacomb.numeric.math;
4
5
6 public final class Random {
7
8 private static int jran = 12345;
9 private static int im = 714025;
10 private static int ia = 1366;
11 private static int ic = 150889;
12
13
14
15 private static double[] cof = {76.18009173, -86.50532033, 24.01409822,
16 -1.231739516, 0.120858003e-2, -0.536382e-5
17 };
18
19
20
21
22 public static double random() {
23 jran = (jran * ia + ic) % im;
24 double ran = (1. * jran) / im;
25 return ran;
26 }
27
28 public static int getSeed() {
29 return jran;
30 }
31
32
33 public static void setSeed(int jr) {
34 jran = jr;
35 }
36
37
38
39 public static double nextRandom() {
40 return random();
41 }
42
43
44 public static double uniformRV() {
45 return random();
46 }
47
48
49
50
51 public static int weightedSample(double[] rw) {
52 int n = rw.length;
53 double a = random();
54 int inew = 0;
55 while ((a -= rw[inew]) > 0 && inew < n-1) {
56 inew++;
57 }
58 return inew;
59 }
60
61
62
63 public static double gaussianRV() {
64 return grv();
65 }
66
67
68 public static double grv() {
69 double r, ran1, ran2, fac, g1;
70 r = -1;
71 ran1 = 0.0;
72 ran2 = 0.0;
73 while (r <= 0.0 || r >= 1.0) {
74 jran = (jran * ia + ic) % im;
75 ran1 = (2. * jran) / im - 1;
76
77 jran = (jran * ia + ic) % im;
78 ran2 = (2. * jran) / im - 1;
79
80 r = ran1 * ran1 + ran2 * ran2;
81 }
82 fac = Math.sqrt(-2. * Math.log(r) / r);
83 g1 = ran1 * fac;
84
85 return g1;
86 }
87
88
89
90 public final static int poissonInt(double d) {
91 return (int)(poidev(d));
92 }
93
94
95
96 public final static double gammln(double xx) {
97 double x = xx - 1.0;
98 double tmp = x + 5.5;
99 tmp -= (x+0.5) * Math.log(tmp);
100 double ser = 1.0;
101 for (int j = 0; j <= 5; j++) {
102 x += 1.0;
103 ser += cof[j]/x;
104 }
105 return -tmp+Math.log(2.50662827465*ser);
106 }
107
108
109
110 public final static double poidev(double xm) {
111
112 double em = -1.;
113 if (xm < 12.0) {
114 double g = Math.exp(-xm);
115 double t = 1.0;
116 do {
117 em += 1.0;
118 t *= uniformRV();
119 } while (t > g);
120
121
122 } else {
123 double sq = Math.sqrt(2.0 * xm);
124 double alxm = Math.log(xm);
125 double g = xm * alxm - gammln(xm + 1.0);
126 double y = 0.;
127 double t = -1.;
128 do {
129 do {
130 y = Math.tan(Math.PI*uniformRV());
131 em= sq * y + xm;
132 } while (em < 0.0);
133
134 em = Math.floor(em);
135 t = 0.9 * (1.0 + y*y) * Math.exp(em * alxm - gammln(em + 1.0) - g);
136 } while (uniformRV() > t);
137 }
138 return em;
139 }
140 }
141
142
143
144