1 static char help[] = "Run Birthday Spacing Tests for PetscRandom.\n\n"; 2 3 #include <petscsys.h> 4 #include <petscviewer.h> 5 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 */ 10 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.; 17 18 PetscFunctionBegin; 19 logLambda = PetscLogReal(lambda); 20 for (i = 0; i < Y; i++) { 21 PetscReal exponent = -lambda + i * logLambda - logFact; 22 23 logFact += PetscLogReal(i + 1); 24 25 p -= PetscExpReal(exponent); 26 } 27 *prob = p; 28 PetscFunctionReturn(PETSC_SUCCESS); 29 } 30 31 int main(int argc, char **argv) 32 { 33 PetscMPIInt size; 34 PetscInt log2d, log2n, t, N, Y; 35 PetscInt64 d, k, *X; 36 size_t n, i; 37 PetscReal lambda, p; 38 PetscRandom random; 39 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(); 50 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 = ((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))); 59 60 PetscCall(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)); 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; 71 72 for (j = 0; j < t; j++, mult *= d) { 73 PetscReal x; 74 PetscInt64 slot; 75 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 } 84 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]); 89 90 PetscCall(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)); 93 94 PetscCall(PetscFree(X)); 95 PetscCall(PetscRandomDestroy(&random)); 96 PetscCall(PetscFinalize()); 97 return 0; 98 } 99 100 /*TEST 101 102 test: 103 args: -t 4 -log2d 7 -log2n 10 104 105 TEST*/ 106