Actual source code: ex3.c
2: static char help[] = "Run Birthday Spacing Tests for PetscRandom.\n\n";
4: #include <petscsys.h>
5: #include <petscviewer.h>
7: /* L'Ecuyer & Simard, 2001.
8: * "On the performance of birthday spacings tests with certain families of random number generators"
9: * https://doi.org/10.1016/S0378-4754(00)00253-6
10: */
12: static int PetscInt64Compare(const void *a, const void *b)
13: {
14: PetscInt64 A = *((const PetscInt64 *)a);
15: PetscInt64 B = *((const PetscInt64 *)b);
16: if (A < B) return -1;
17: if (A == B) return 0;
18: return 1;
19: }
21: static PetscErrorCode PoissonTailProbability(PetscReal lambda, PetscInt Y, PetscReal *prob)
22: {
23: PetscReal p = 1.;
24: PetscReal logLambda;
25: PetscInt i;
26: PetscReal logFact = 0.;
28: logLambda = PetscLogReal(lambda);
29: for (i = 0; i < Y; i++) {
30: PetscReal exponent = -lambda + i * logLambda - logFact;
32: logFact += PetscLogReal(i+1);
34: p -= PetscExpReal(exponent);
35: }
36: *prob = p;
37: return 0;
38: }
40: int main(int argc, char **argv)
41: {
42: PetscMPIInt size;
43: PetscInt log2d, log2n, t, N, Y;
44: PetscInt64 d, k, *X;
45: size_t n, i;
46: PetscReal lambda, p;
47: PetscRandom random;
50: PetscInitialize(&argc,&argv,NULL,help);
52: t = 8;
53: log2d = 7;
54: log2n = 20;
55: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Birthday spacing test parameters", "PetscRandom");
56: PetscOptionsInt("-t", "t, the dimension of the sample space", "ex3.c", t, &t, NULL);
57: PetscOptionsInt("-log2d", "The log of d, the number of bins per direction", "ex3.c", log2d, &log2d, NULL);
58: PetscOptionsInt("-log2n", "The log of n, the number of samples per process", "ex3.c", log2n, &log2n, NULL);
59: PetscOptionsEnd();
62: d = ((PetscInt64) 1) << log2d;
63: k = ((PetscInt64) 1) << (log2d * t);
65: n = ((size_t) 1) << log2n;
66: MPI_Comm_size(PETSC_COMM_WORLD,&size);
67: N = size;
68: lambda = PetscPowRealInt(2.,(3 * log2n - (2 + log2d * t)));
70: PetscPrintf(PETSC_COMM_WORLD,"Generating %zu samples (%g GB) per process in a %" PetscInt_FMT " dimensional space with %" PetscInt64_FMT " bins per dimension.\n", n, (n*1.e-9)*sizeof(PetscInt64), t, d);
71: PetscPrintf(PETSC_COMM_WORLD,"Expected spacing collisions per process %g (%g total).\n", (double) lambda, (double)(N * lambda));
72: PetscRandomCreate(PETSC_COMM_WORLD,&random);
73: PetscRandomSetFromOptions(random);
74: PetscRandomSetInterval(random,0.0,1.0);
75: PetscRandomView(random,PETSC_VIEWER_STDOUT_WORLD);
76: PetscMalloc1(n,&X);
77: for (i = 0; i < n; i++) {
78: PetscInt j;
79: PetscInt64 bin = 0;
80: PetscInt64 mult = 1;
82: for (j = 0; j < t; j++, mult *= d) {
83: PetscReal x;
84: PetscInt64 slot;
86: PetscRandomGetValueReal(random,&x);
87: /* worried about precision here */
88: slot = (PetscInt64) (x * d);
89: bin += mult * slot;
90: }
92: X[i] = bin;
93: }
95: qsort(X, n, sizeof(PetscInt64), PetscInt64Compare);
96: for (i = 0; i < n - 1; i++) {
97: X[i] = X[i + 1] - X[i];
98: }
99: qsort(X, n - 1, sizeof(PetscInt64), PetscInt64Compare);
100: for (i = 0, Y = 0; i < n - 2; i++) {
101: Y += (X[i + 1] == X[i]);
102: }
104: MPI_Allreduce(MPI_IN_PLACE, &Y, 1, MPIU_INT, MPI_SUM, MPI_COMM_WORLD);
105: PoissonTailProbability(N*lambda,Y,&p);
106: PetscPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT " total collisions counted: that many or more should occur with probabilty %g.\n",Y,(double)p);
108: PetscFree(X);
109: PetscRandomDestroy(&random);
110: PetscFinalize();
111: return 0;
112: }
114: /*TEST
116: test:
117: args: -t 4 -log2d 7 -log2n 10
119: TEST*/