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);
 51:   t     = 8;
 52:   log2d = 7;
 53:   log2n = 20;
 54:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Birthday spacing test parameters", "PetscRandom");
 55:   PetscOptionsInt("-t", "t, the dimension of the sample space", "ex3.c", t, &t, NULL);
 56:   PetscOptionsInt("-log2d", "The log of d, the number of bins per direction", "ex3.c", log2d, &log2d, NULL);
 57:   PetscOptionsInt("-log2n", "The log of n, the number of samples per process", "ex3.c", log2n, &log2n, NULL);
 58:   PetscOptionsEnd();

 61:   d = ((PetscInt64)1) << log2d;
 62:   k = ((PetscInt64)1) << (log2d * t);
 64:   n = ((size_t)1) << log2n;
 65:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 66:   N      = size;
 67:   lambda = PetscPowRealInt(2., (3 * log2n - (2 + log2d * t)));

 69:   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);
 70:   PetscPrintf(PETSC_COMM_WORLD, "Expected spacing collisions per process %g (%g total).\n", (double)lambda, (double)(N * lambda));
 71:   PetscRandomCreate(PETSC_COMM_WORLD, &random);
 72:   PetscRandomSetFromOptions(random);
 73:   PetscRandomSetInterval(random, 0.0, 1.0);
 74:   PetscRandomView(random, PETSC_VIEWER_STDOUT_WORLD);
 75:   PetscMalloc1(n, &X);
 76:   for (i = 0; i < n; i++) {
 77:     PetscInt   j;
 78:     PetscInt64 bin  = 0;
 79:     PetscInt64 mult = 1;

 81:     for (j = 0; j < t; j++, mult *= d) {
 82:       PetscReal  x;
 83:       PetscInt64 slot;

 85:       PetscRandomGetValueReal(random, &x);
 86:       /* worried about precision here */
 87:       slot = (PetscInt64)(x * d);
 88:       bin += mult * slot;
 89:     }
 91:     X[i] = bin;
 92:   }

 94:   qsort(X, n, sizeof(PetscInt64), PetscInt64Compare);
 95:   for (i = 0; i < n - 1; i++) X[i] = X[i + 1] - X[i];
 96:   qsort(X, n - 1, sizeof(PetscInt64), PetscInt64Compare);
 97:   for (i = 0, Y = 0; i < n - 2; i++) Y += (X[i + 1] == X[i]);

 99:   MPI_Allreduce(MPI_IN_PLACE, &Y, 1, MPIU_INT, MPI_SUM, MPI_COMM_WORLD);
100:   PoissonTailProbability(N * lambda, Y, &p);
101:   PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT " total collisions counted: that many or more should occur with probability %g.\n", Y, (double)p);

103:   PetscFree(X);
104:   PetscRandomDestroy(&random);
105:   PetscFinalize();
106:   return 0;
107: }

109: /*TEST

111:   test:
112:     args: -t 4 -log2d 7 -log2n 10

114: TEST*/