1 2 static char help[] = "Run Birthday Spacing Tests for PetscRandom.\n\n"; 3 4 #include <petscsys.h> 5 #include <petscviewer.h> 6 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 */ 11 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 } 20 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.; 27 28 PetscFunctionBegin; 29 logLambda = PetscLogReal(lambda); 30 for (i = 0; i < Y; i++) { 31 PetscReal exponent = -lambda + i * logLambda - logFact; 32 33 logFact += PetscLogReal(i+1); 34 35 p -= PetscExpReal(exponent); 36 } 37 *prob = p; 38 PetscFunctionReturn(0); 39 } 40 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; 49 PetscErrorCode ierr; 50 51 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 52 53 t = 8; 54 log2d = 7; 55 log2n = 20; 56 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Birthday spacing test parameters", "PetscRandom");CHKERRQ(ierr); 57 ierr = PetscOptionsInt("-t", "t, the dimension of the sample space", "ex3.c", t, &t, NULL);CHKERRQ(ierr); 58 ierr = PetscOptionsInt("-log2d", "The log of d, the number of bins per direction", "ex3.c", log2d, &log2d, NULL);CHKERRQ(ierr); 59 ierr = PetscOptionsInt("-log2n", "The log of n, the number of samples per process", "ex3.c", log2n, &log2n, NULL);CHKERRQ(ierr); 60 ierr = PetscOptionsEnd();CHKERRQ(ierr); 61 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 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 68 N = size; 69 lambda = PetscPowRealInt(2.,(3 * log2n - (2 + log2d * t))); 70 71 ierr = 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);CHKERRQ(ierr); 72 ierr = PetscPrintf(PETSC_COMM_WORLD,"Expected spacing collisions per process %f (%f total).\n", (double) lambda, N * lambda);CHKERRQ(ierr); 73 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&random);CHKERRQ(ierr); 74 ierr = PetscRandomSetFromOptions(random);CHKERRQ(ierr); 75 ierr = PetscRandomSetInterval(random,0.0,1.0);CHKERRQ(ierr); 76 ierr = PetscRandomView(random,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 77 ierr = PetscMalloc1(n,&X);CHKERRQ(ierr); 78 for (i = 0; i < n; i++) { 79 PetscInt j; 80 PetscInt64 bin = 0; 81 PetscInt64 mult = 1; 82 83 for (j = 0; j < t; j++, mult *= d) { 84 PetscReal x; 85 PetscInt64 slot; 86 87 ierr = PetscRandomGetValueReal(random,&x);CHKERRQ(ierr); 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 } 95 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 } 104 105 ierr = MPI_Allreduce(MPI_IN_PLACE, &Y, 1, MPIU_INT, MPI_SUM, MPI_COMM_WORLD);CHKERRQ(ierr); 106 ierr = PoissonTailProbability(N*lambda,Y,&p);CHKERRQ(ierr); 107 ierr = PetscPrintf(PETSC_COMM_WORLD,"%D total collisions counted: that many or more should occur with probabilty %g.\n",Y,(double)p);CHKERRQ(ierr); 108 109 ierr = PetscFree(X);CHKERRQ(ierr); 110 ierr = PetscRandomDestroy(&random);CHKERRQ(ierr); 111 ierr = PetscFinalize(); 112 return ierr; 113 } 114 115 /*TEST 116 117 test: 118 args: -t 4 -log2d 7 -log2n 10 119 120 TEST*/ 121