1 package nl.tudelft.simulation.jstats.math;
2
3 import org.djutils.exceptions.Throw;
4
5
6
7
8
9
10
11
12
13
14
15
16 public final class ProbMath
17 {
18
19 public static final double[] FACTORIAL_DOUBLE;
20
21
22 public static final long[] FACTORIAL_LONG;
23
24
25 public static final double[] POW2_DOUBLE;
26
27
28 public static final long[] POW2_LONG;
29
30
31 static
32 {
33 FACTORIAL_DOUBLE = new double[171];
34 FACTORIAL_DOUBLE[0] = 1.0;
35 double d = 1.0;
36 for (int i = 1; i <= 170; i++)
37 {
38 d = d * i;
39 FACTORIAL_DOUBLE[i] = d;
40 }
41
42 FACTORIAL_LONG = new long[21];
43 FACTORIAL_LONG[0] = 1L;
44 long n = 1L;
45 for (int i = 1; i <= 20; i++)
46 {
47 n = n * i;
48 FACTORIAL_LONG[i] = n;
49 }
50
51 POW2_DOUBLE = new double[1024];
52 POW2_DOUBLE[0] = 1.0;
53 for (int i = 1; i < 1024; i++)
54 {
55 POW2_DOUBLE[i] = 2.0 * POW2_DOUBLE[i - 1];
56 }
57
58 POW2_LONG = new long[63];
59 POW2_LONG[0] = 1L;
60 for (int i = 1; i < 63; i++)
61 {
62 POW2_LONG[i] = 2L * POW2_LONG[i - 1];
63 }
64 }
65
66
67
68
69 private ProbMath()
70 {
71
72 }
73
74
75
76
77
78
79 public static double factorial(final int n)
80 {
81 Throw.when(n < 0, IllegalArgumentException.class, "n! with n<0 is invalid");
82 Throw.when(n > 170, IllegalArgumentException.class, "n! with n>170 is larger than Double.MAX_VALUE");
83 return FACTORIAL_DOUBLE[n];
84 }
85
86
87
88
89
90
91 public static double factorial(final long n)
92 {
93 Throw.when(n < 0, IllegalArgumentException.class, "n! with n<0 is invalid");
94 Throw.when(n > 170, IllegalArgumentException.class, "n! with n>170 is larger than Double.MAX_VALUE");
95 return FACTORIAL_DOUBLE[(int) n];
96 }
97
98
99
100
101
102
103 public static long fac(final int n)
104 {
105 Throw.when(n < 0, IllegalArgumentException.class, "n! with n<0 is invalid");
106 Throw.when(n > 20, IllegalArgumentException.class, "n! with n>20 is more than 64 bits long");
107 return FACTORIAL_LONG[n];
108 }
109
110
111
112
113
114
115 public static long fac(final long n)
116 {
117 Throw.when(n < 0, IllegalArgumentException.class, "n! with n<0 is invalid");
118 Throw.when(n > 20, IllegalArgumentException.class, "n! with n>20 is more than 64 bits long");
119 return FACTORIAL_LONG[(int) n];
120 }
121
122
123
124
125
126
127
128 public static double permutations(final int n, final int k)
129 {
130 if (k > n)
131 {
132 throw new IllegalArgumentException("permutations of (n,k) with k>n");
133 }
134 return factorial(n) / factorial(n - k);
135 }
136
137
138
139
140
141
142
143 public static double permutations(final long n, final long k)
144 {
145 if (k > n)
146 {
147 throw new IllegalArgumentException("permutations of (n,k) with k>n");
148 }
149 return factorial(n) / factorial(n - k);
150 }
151
152
153
154
155
156
157
158 public static long perm(final int n, final int k)
159 {
160 if (k > n)
161 {
162 throw new IllegalArgumentException("permutations of (n,k) with k>n");
163 }
164 return fac(n) / fac(n - k);
165 }
166
167
168
169
170
171
172
173 public static long perm(final long n, final long k)
174 {
175 if (k > n)
176 {
177 throw new IllegalArgumentException("permutations of (n,k) with k>n");
178 }
179 return fac(n) / fac(n - k);
180 }
181
182
183
184
185
186
187
188 public static double combinations(final int n, final int k)
189 {
190 if (k > n)
191 {
192 throw new IllegalArgumentException("combinations of (n,k) with k>n");
193 }
194 return factorial(n) / (factorial(k) * factorial(n - k));
195 }
196
197
198
199
200
201
202
203 public static double combinations(final long n, final long k)
204 {
205 if (k > n)
206 {
207 throw new IllegalArgumentException("combinations of (n,k) with k>n");
208 }
209 return factorial(n) / (factorial(k) * factorial(n - k));
210 }
211
212
213
214
215
216
217
218 public static long comb(final int n, final int k)
219 {
220 if (k > n)
221 {
222 throw new IllegalArgumentException("combinations of (n,k) with k>n");
223 }
224 return fac(n) / (fac(k) * fac(n - k));
225 }
226
227
228
229
230
231
232
233 public static long comb(final long n, final long k)
234 {
235 if (k > n)
236 {
237 throw new IllegalArgumentException("combinations of (n,k) with k>n");
238 }
239 return fac(n) / (fac(k) * fac(n - k));
240 }
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258 public static double erf(final double z)
259 {
260 double zpos = Math.abs(z);
261 if (zpos < 0.5)
262 {
263 return erfSmall(z);
264 }
265 if (zpos > 3.8)
266 {
267 return erfBig(z);
268 }
269 return erfTaylor(z);
270 }
271
272
273
274
275
276
277
278
279
280
281 private static double erfSmall(final double z)
282 {
283 double zpos = Math.abs(z);
284
285 double sum = zpos
286 + 2.0 * Math.pow(zpos, 3) / 3.0
287 + 4.0 * Math.pow(zpos, 5) / 15.0
288 + 8.0 * Math.pow(zpos, 7) / 105.0
289 + 16.0 * Math.pow(zpos, 9) / 945.0
290 + 32.0 * Math.pow(zpos, 11) / 10395.0
291 + 64.0 * Math.pow(zpos, 13) / 135135.0
292 + 128.0 * Math.pow(zpos, 15) / 2027025.0
293 + 256.0 * Math.pow(zpos, 17) / 34459425.0
294 + 512.0 * Math.pow(zpos, 19) / 654729075.0
295 + 1024.0 * Math.pow(zpos, 21) / 13749310575.0;
296
297 return Math.signum(z) * sum * 2.0 * Math.exp(-zpos * zpos) / Math.sqrt(Math.PI);
298 }
299
300
301
302
303
304
305
306
307
308
309
310 private static double erfBig(final double z)
311 {
312 double zpos = Math.abs(z);
313
314 double sum = 1.0 / zpos
315 - (1.0 / 2.0) * Math.pow(zpos, -3)
316 + (3.0 / 4.0) * Math.pow(zpos, -5)
317 - (15.0 / 8.0) * Math.pow(zpos, -7)
318 + (105.0 / 16.0) * Math.pow(zpos, -9)
319 - (945.0 / 32.0) * Math.pow(zpos, -11)
320 + (10395.0 / 64.0) * Math.pow(zpos, -13)
321 - (135135.0 / 128.0) * Math.pow(zpos, -15)
322 + (2027025.0 / 256.0) * Math.pow(zpos, -17);
323
324 return Math.signum(z) * (1.0 - sum * Math.exp(-zpos * zpos) / Math.sqrt(Math.PI));
325 }
326
327
328
329
330
331
332
333
334
335
336
337 private static double erfTaylor(final double z)
338 {
339 double zpos = Math.abs(z);
340 double d = zpos;
341 double zpow = zpos;
342 for (int i = 1; i < 64; i++)
343 {
344
345 zpow *= zpos * zpos;
346 double term = zpow / ((2.0 * i + 1.0) * ProbMath.factorial(i));
347 d += term * ((i & 1) == 0 ? 1 : -1);
348 if (term < 1E-16)
349 {
350 break;
351 }
352 }
353 return Math.signum(z) * d * 2.0 / Math.sqrt(Math.PI);
354 }
355
356
357
358
359
360
361
362 public static double erfInv(final double y)
363 {
364 double ax, t, ret;
365 ax = Math.abs(y);
366
367
368
369
370
371 if (ax <= 0.75)
372 {
373
374 double[] p = new double[] {-13.0959967422, 26.785225760, -9.289057635};
375 double[] q = new double[] {-12.0749426297, 30.960614529, -17.149977991, 1.00000000};
376
377 t = ax * ax - 0.75 * 0.75;
378 ret = ax * (p[0] + t * (p[1] + t * p[2])) / (q[0] + t * (q[1] + t * (q[2] + t * q[3])));
379
380 }
381 else if (ax >= 0.75 && ax <= 0.9375)
382 {
383 double[] p = new double[] {-.12402565221, 1.0688059574, -1.9594556078, .4230581357};
384 double[] q = new double[] {-.08827697997, .8900743359, -2.1757031196, 1.0000000000};
385
386
387
388
389
390 t = ax * ax - 0.9375 * 0.9375;
391 ret = ax * (p[0] + t * (p[1] + t * (p[2] + t * p[3]))) / (q[0] + t * (q[1] + t * (q[2] + t * q[3])));
392
393 }
394 else if (ax >= 0.9375 && ax <= (1.0 - 1.0e-9))
395 {
396 double[] p =
397 new double[] {.1550470003116, 1.382719649631, .690969348887, -1.128081391617, .680544246825, -.16444156791};
398 double[] q = new double[] {.155024849822, 1.385228141995, 1.000000000000};
399
400
401
402
403
404 t = 1.0 / Math.sqrt(-Math.log(1.0 - ax));
405 ret = (p[0] / t + p[1] + t * (p[2] + t * (p[3] + t * (p[4] + t * p[5])))) / (q[0] + t * (q[1] + t * (q[2])));
406 }
407 else
408 {
409 ret = Double.POSITIVE_INFINITY;
410 }
411
412 return Math.signum(y) * ret;
413 }
414
415
416 private static final double[] GAMMALN_COF = {76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155,
417 0.1208650973866179e-2, -0.5395239384953e-5};
418
419
420
421
422
423
424
425 public static double gammaln(final double xx)
426 {
427 Throw.when(xx < 0, IllegalArgumentException.class, "gamma function not defined for real values < 0");
428 double x, y, tmp, ser;
429 x = xx;
430 y = x;
431 tmp = x + 5.5;
432 tmp -= (x + 0.5) * Math.log(tmp);
433 ser = 1.000000000190015;
434 for (int j = 0; j <= 5; j++)
435 {
436 ser += GAMMALN_COF[j] / ++y;
437 }
438 return -tmp + Math.log(2.5066282746310005 * ser / x);
439 }
440
441
442
443
444
445
446
447 public static double gamma(final double x)
448 {
449 return Math.exp(gammaln(x));
450 }
451
452
453
454
455
456
457
458
459 public static double beta(final double z, final double w)
460 {
461 Throw.when(z < 0 || w < 0, IllegalArgumentException.class, "beta function not defined for negative arguments");
462 return Math.exp(gammaln(z) + gammaln(w) - gammaln(z + w));
463 }
464
465 }