Actual source code: ex3.c
1: static char help[] = "Run Birthday Spacing Tests for PetscRandom.\n\n";
3: #include <petscsys.h>
4: #include <petscviewer.h>
6: /* L'Ecuyer & Simard, 2001.
7: * "On the performance of birthday spacings tests with certain families of random number generators"
8: * https://doi.org/10.1016/S0378-4754(00)00253-6
9: */
11: static PetscErrorCode PoissonTailProbability(PetscReal lambda, PetscInt Y, PetscReal *prob)
12: {
13: PetscReal p = 1.;
14: PetscReal logLambda;
15: PetscInt i;
16: PetscReal logFact = 0.;
18: PetscFunctionBegin;
19: logLambda = PetscLogReal(lambda);
20: for (i = 0; i < Y; i++) {
21: PetscReal exponent = -lambda + i * logLambda - logFact;
23: logFact += PetscLogReal(i + 1);
25: p -= PetscExpReal(exponent);
26: }
27: *prob = p;
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: int main(int argc, char **argv)
32: {
33: PetscMPIInt size;
34: PetscInt log2d, log2n, t, N, Y;
35: PetscInt64 d, k, *X;
36: PetscInt n, i;
37: PetscReal lambda, p;
38: PetscRandom random;
40: PetscFunctionBeginUser;
41: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
42: t = 8;
43: log2d = 7;
44: log2n = 20;
45: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Birthday spacing test parameters", "PetscRandom");
46: PetscCall(PetscOptionsInt("-t", "t, the dimension of the sample space", "ex3.c", t, &t, NULL));
47: PetscCall(PetscOptionsInt("-log2d", "The log of d, the number of bins per direction", "ex3.c", log2d, &log2d, NULL));
48: PetscCall(PetscOptionsInt("-log2n", "The log of n, the number of samples per process", "ex3.c", log2n, &log2n, NULL));
49: PetscOptionsEnd();
51: PetscCheck((size_t)log2d * t <= sizeof(PetscInt64) * 8 - 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of bins (2^%" PetscInt_FMT ") is too big for PetscInt64.", log2d * t);
52: d = ((PetscInt64)1) << log2d;
53: k = ((PetscInt64)1) << (log2d * t);
54: PetscCheck((size_t)log2n <= sizeof(size_t) * 8 - 1, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "The number of samples per process (2^%" PetscInt_FMT ") is too big for size_t.", log2n);
55: n = (PetscInt)(((size_t)1) << log2n);
56: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
57: N = size;
58: lambda = PetscPowRealInt(2., 3 * log2n - (2 + log2d * t));
60: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Generating %" PetscInt_FMT " 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));
61: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected spacing collisions per process %g (%g total).\n", (double)lambda, (double)(N * lambda)));
62: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &random));
63: PetscCall(PetscRandomSetFromOptions(random));
64: PetscCall(PetscRandomSetInterval(random, 0.0, 1.0));
65: PetscCall(PetscRandomView(random, PETSC_VIEWER_STDOUT_WORLD));
66: PetscCall(PetscMalloc1(n, &X));
67: for (i = 0; i < n; i++) {
68: PetscInt j;
69: PetscInt64 bin = 0;
70: PetscInt64 mult = 1;
72: for (j = 0; j < t; j++, mult *= d) {
73: PetscReal x;
74: PetscInt64 slot;
76: PetscCall(PetscRandomGetValueReal(random, &x));
77: /* worried about precision here */
78: slot = (PetscInt64)(x * d);
79: bin += mult * slot;
80: }
81: PetscCheck(bin < k, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Generated point in bin %" PetscInt64_FMT ", but only %" PetscInt64_FMT " bins", bin, k);
82: X[i] = bin;
83: }
85: PetscCall(PetscSortInt64(n, X));
86: for (i = 0; i < n - 1; i++) X[i] = X[i + 1] - X[i];
87: PetscCall(PetscSortInt64(n - 1, X));
88: for (i = 0, Y = 0; i < n - 2; i++) Y += (X[i + 1] == X[i]);
90: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &Y, 1, MPIU_INT, MPI_SUM, MPI_COMM_WORLD));
91: PetscCall(PoissonTailProbability(N * lambda, Y, &p));
92: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT " total collisions counted: that many or more should occur with probability %g.\n", Y, (double)p));
94: PetscCall(PetscFree(X));
95: PetscCall(PetscRandomDestroy(&random));
96: PetscCall(PetscFinalize());
97: return 0;
98: }
100: /*TEST
102: test:
103: args: -t 4 -log2d 7 -log2n 10
105: TEST*/