1: static char help[] = "Run Birthday Spacing Tests for PetscRandom.\n\n";

  3: #include <petscsys.h>

  5: /* L'Ecuyer & Simard, 2001.
  6:  * "On the performance of birthday spacings tests with certain families of random number generators"
  7:  * https://doi.org/10.1016/S0378-4754(00)00253-6
  8:  */

 10: static PetscErrorCode PoissonTailProbability(PetscReal lambda, PetscInt Y, PetscReal *prob)
 11: {
 12:   PetscReal p = 1.;
 13:   PetscReal logLambda;
 14:   PetscInt  i;
 15:   PetscReal logFact = 0.;

 17:   PetscFunctionBegin;
 18:   logLambda = PetscLogReal(lambda);
 19:   for (i = 0; i < Y; i++) {
 20:     PetscReal exponent = -lambda + i * logLambda - logFact;

 22:     logFact += PetscLogReal(i + 1);

 24:     p -= PetscExpReal(exponent);
 25:   }
 26:   *prob = p;
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 28: }

 30: int main(int argc, char **argv)
 31: {
 32:   PetscMPIInt size;
 33:   PetscInt    log2d, log2n, t, N, Y;
 34:   PetscInt64  d, k, *X;
 35:   PetscInt    n, i;
 36:   PetscReal   lambda, p;
 37:   PetscRandom random;

 39:   PetscFunctionBeginUser;
 40:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 41:   t     = 8;
 42:   log2d = 7;
 43:   log2n = 20;
 44:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Birthday spacing test parameters", "PetscRandom");
 45:   PetscCall(PetscOptionsInt("-t", "t, the dimension of the sample space", "ex3.c", t, &t, NULL));
 46:   PetscCall(PetscOptionsInt("-log2d", "The log of d, the number of bins per direction", "ex3.c", log2d, &log2d, NULL));
 47:   PetscCall(PetscOptionsInt("-log2n", "The log of n, the number of samples per process", "ex3.c", log2n, &log2n, NULL));
 48:   PetscOptionsEnd();

 50:   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);
 51:   d = ((PetscInt64)1) << log2d;
 52:   k = ((PetscInt64)1) << (log2d * t);
 53:   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);
 54:   n = (PetscInt)(((size_t)1) << log2n);
 55:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 56:   N      = size;
 57:   lambda = PetscPowRealInt(2., 3 * log2n - (2 + log2d * t));

 59:   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));
 60:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected spacing collisions per process %g (%g total).\n", (double)lambda, (double)(N * lambda)));
 61:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &random));
 62:   PetscCall(PetscRandomSetFromOptions(random));
 63:   PetscCall(PetscRandomSetInterval(random, 0.0, 1.0));
 64:   PetscCall(PetscRandomView(random, PETSC_VIEWER_STDOUT_WORLD));
 65:   PetscCall(PetscMalloc1(n, &X));
 66:   for (i = 0; i < n; i++) {
 67:     PetscInt64 bin  = 0;
 68:     PetscInt64 mult = 1;

 70:     for (PetscInt j = 0; j < t; j++, mult *= d) {
 71:       PetscReal  x;
 72:       PetscInt64 slot;

 74:       PetscCall(PetscRandomGetValueReal(random, &x));
 75:       /* worried about precision here */
 76:       slot = (PetscInt64)(x * d);
 77:       bin += mult * slot;
 78:     }
 79:     PetscCheck(bin < k, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Generated point in bin %" PetscInt64_FMT ", but only %" PetscInt64_FMT " bins", bin, k);
 80:     X[i] = bin;
 81:   }

 83:   PetscCall(PetscSortInt64(n, X));
 84:   for (i = 0; i < n - 1; i++) X[i] = X[i + 1] - X[i];
 85:   PetscCall(PetscSortInt64(n - 1, X));
 86:   for (i = 0, Y = 0; i < n - 2; i++) Y += (X[i + 1] == X[i]);

 88:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &Y, 1, MPIU_INT, MPI_SUM, MPI_COMM_WORLD));
 89:   PetscCall(PoissonTailProbability(N * lambda, Y, &p));
 90:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT " total collisions counted: that many or more should occur with probability %g.\n", Y, (double)p));

 92:   PetscCall(PetscFree(X));
 93:   PetscCall(PetscRandomDestroy(&random));
 94:   PetscCall(PetscFinalize());
 95:   return 0;
 96: }

 98: /*TEST

100:   test:
101:     args: -t 4 -log2d 7 -log2n 10

103: TEST*/