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.;

 29:   logLambda = PetscLogReal(lambda);
 30:   for (i = 0; i < Y; i++) {
 31:     PetscReal exponent = -lambda + i * logLambda - logFact;

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

 35:     p -= PetscExpReal(exponent);
 36:   }
 37:   *prob = p;
 38:   return(0);
 39: }

 41: int main(int argc, char **argv)
 42: {
 43:   PetscMPIInt    size;
 44:   PetscInt       log2d, log2n, t, N, Y;
 45:   PetscInt64     d, k, *X;
 46:   size_t         n, i;
 47:   PetscReal      lambda, p;
 48:   PetscRandom    random;

 51:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;

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

 62:   if ((size_t)log2d * t > sizeof(PetscInt64) * 8 - 2) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE, "The number of bins (2^%D) is too big for PetscInt64.", log2d * t);
 63:   d = ((PetscInt64) 1) << log2d;
 64:   k = ((PetscInt64) 1) << (log2d * t);
 65:   if ((size_t)log2n > sizeof(size_t) * 8 - 1) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE, "The number of samples per process (2^%D) is too big for size_t.", log2n);
 66:   n = ((size_t) 1) << log2n;
 67:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 68:   N    = size;
 69:   lambda = PetscPowRealInt(2.,(3 * log2n - (2 + log2d * t)));

 71:   PetscPrintf(PETSC_COMM_WORLD,"Generating %zu samples (%g GB) per process in a %D dimensional space with %" PetscInt64_FMT " bins per dimension.\n", n, (n*1.e-9)*sizeof(PetscInt64), t, d);
 72:   PetscPrintf(PETSC_COMM_WORLD,"Expected spacing collisions per process %f (%f total).\n", (double) lambda, N * lambda);
 73:   PetscRandomCreate(PETSC_COMM_WORLD,&random);
 74:   PetscRandomSetFromOptions(random);
 75:   PetscRandomSetInterval(random,0.0,1.0);
 76:   PetscRandomView(random,PETSC_VIEWER_STDOUT_WORLD);
 77:   PetscMalloc1(n,&X);
 78:   for (i = 0; i < n; i++) {
 79:     PetscInt   j;
 80:     PetscInt64 bin  = 0;
 81:     PetscInt64 mult = 1;

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

 87:       PetscRandomGetValueReal(random,&x);
 88:       /* worried about precision here */
 89:       slot = (PetscInt64) (x * d);
 90:       bin += mult * slot;
 91:     }
 92:     if (bin >= k) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Generated point in bin %" PetscInt64_FMT ", but only %" PetscInt64_FMT " bins\n",bin,k);
 93:     X[i] = bin;
 94:   }

 96:   qsort(X, n, sizeof(PetscInt64), PetscInt64Compare);
 97:   for (i = 0; i < n - 1; i++) {
 98:     X[i] = X[i + 1] - X[i];
 99:   }
100:   qsort(X, n - 1, sizeof(PetscInt64), PetscInt64Compare);
101:   for (i = 0, Y = 0; i < n - 2; i++) {
102:     Y += (X[i + 1] == X[i]);
103:   }

105:   MPI_Allreduce(MPI_IN_PLACE, &Y, 1, MPIU_INT, MPI_SUM, MPI_COMM_WORLD);
106:   PoissonTailProbability(N*lambda,Y,&p);
107:   PetscPrintf(PETSC_COMM_WORLD,"%D total collisions counted: that many or more should occur with probabilty %g.\n",Y,(double)p);

109:   PetscFree(X);
110:   PetscRandomDestroy(&random);
111:   PetscFinalize();
112:   return ierr;
113: }

115: /*TEST

117:   test:
118:     args: -t 4 -log2d 7 -log2n 10

120: TEST*/