1c4762a1bSJed Brown 2c4762a1bSJed Brown static char help[] = "Run Birthday Spacing Tests for PetscRandom.\n\n"; 3c4762a1bSJed Brown 4c4762a1bSJed Brown #include <petscsys.h> 5c4762a1bSJed Brown #include <petscviewer.h> 6c4762a1bSJed Brown 7c4762a1bSJed Brown /* L'Ecuyer & Simard, 2001. 8c4762a1bSJed Brown * "On the performance of birthday spacings tests with certain families of random number generators" 9c4762a1bSJed Brown * https://doi.org/10.1016/S0378-4754(00)00253-6 10c4762a1bSJed Brown */ 11c4762a1bSJed Brown 12c4762a1bSJed Brown static int PetscInt64Compare(const void *a, const void *b) 13c4762a1bSJed Brown { 14c4762a1bSJed Brown PetscInt64 A = *((const PetscInt64 *)a); 15c4762a1bSJed Brown PetscInt64 B = *((const PetscInt64 *)b); 16c4762a1bSJed Brown if (A < B) return -1; 17c4762a1bSJed Brown if (A == B) return 0; 18c4762a1bSJed Brown return 1; 19c4762a1bSJed Brown } 20c4762a1bSJed Brown 21c4762a1bSJed Brown static PetscErrorCode PoissonTailProbability(PetscReal lambda, PetscInt Y, PetscReal *prob) 22c4762a1bSJed Brown { 23c4762a1bSJed Brown PetscReal p = 1.; 24c4762a1bSJed Brown PetscReal logLambda; 25c4762a1bSJed Brown PetscInt i; 26c4762a1bSJed Brown PetscReal logFact = 0.; 27c4762a1bSJed Brown 28c4762a1bSJed Brown PetscFunctionBegin; 29c4762a1bSJed Brown logLambda = PetscLogReal(lambda); 30c4762a1bSJed Brown for (i = 0; i < Y; i++) { 31c4762a1bSJed Brown PetscReal exponent = -lambda + i * logLambda - logFact; 32c4762a1bSJed Brown 33c4762a1bSJed Brown logFact += PetscLogReal(i+1); 34c4762a1bSJed Brown 35c4762a1bSJed Brown p -= PetscExpReal(exponent); 36c4762a1bSJed Brown } 37c4762a1bSJed Brown *prob = p; 38c4762a1bSJed Brown PetscFunctionReturn(0); 39c4762a1bSJed Brown } 40c4762a1bSJed Brown 41c4762a1bSJed Brown int main(int argc, char **argv) 42c4762a1bSJed Brown { 43c4762a1bSJed Brown PetscMPIInt size; 44c4762a1bSJed Brown PetscInt log2d, log2n, t, N, Y; 45c4762a1bSJed Brown PetscInt64 d, k, *X; 46c4762a1bSJed Brown size_t n, i; 47c4762a1bSJed Brown PetscReal lambda, p; 48c4762a1bSJed Brown PetscRandom random; 49c4762a1bSJed Brown PetscErrorCode ierr; 50c4762a1bSJed Brown 51*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscInitialize(&argc,&argv,NULL,help)); 52c4762a1bSJed Brown 53c4762a1bSJed Brown t = 8; 54c4762a1bSJed Brown log2d = 7; 55c4762a1bSJed Brown log2n = 20; 56c4762a1bSJed Brown ierr = PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Birthday spacing test parameters", "PetscRandom");CHKERRQ(ierr); 575f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsInt("-t", "t, the dimension of the sample space", "ex3.c", t, &t, NULL)); 585f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsInt("-log2d", "The log of d, the number of bins per direction", "ex3.c", log2d, &log2d, NULL)); 595f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsInt("-log2n", "The log of n, the number of samples per process", "ex3.c", log2n, &log2n, NULL)); 60c4762a1bSJed Brown ierr = PetscOptionsEnd();CHKERRQ(ierr); 61c4762a1bSJed Brown 622c71b3e2SJacob Faibussowitsch PetscCheckFalse((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); 63c4762a1bSJed Brown d = ((PetscInt64) 1) << log2d; 64c4762a1bSJed Brown k = ((PetscInt64) 1) << (log2d * t); 652c71b3e2SJacob Faibussowitsch PetscCheckFalse((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); 66c4762a1bSJed Brown n = ((size_t) 1) << log2n; 675f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 68c4762a1bSJed Brown N = size; 69c4762a1bSJed Brown lambda = PetscPowRealInt(2.,(3 * log2n - (2 + log2d * t))); 70c4762a1bSJed Brown 715f80ce2aSJacob Faibussowitsch CHKERRQ(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)); 725f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Expected spacing collisions per process %g (%g total).\n", (double) lambda, (double)(N * lambda))); 735f80ce2aSJacob Faibussowitsch CHKERRQ(PetscRandomCreate(PETSC_COMM_WORLD,&random)); 745f80ce2aSJacob Faibussowitsch CHKERRQ(PetscRandomSetFromOptions(random)); 755f80ce2aSJacob Faibussowitsch CHKERRQ(PetscRandomSetInterval(random,0.0,1.0)); 765f80ce2aSJacob Faibussowitsch CHKERRQ(PetscRandomView(random,PETSC_VIEWER_STDOUT_WORLD)); 775f80ce2aSJacob Faibussowitsch CHKERRQ(PetscMalloc1(n,&X)); 78c4762a1bSJed Brown for (i = 0; i < n; i++) { 79c4762a1bSJed Brown PetscInt j; 80c4762a1bSJed Brown PetscInt64 bin = 0; 81c4762a1bSJed Brown PetscInt64 mult = 1; 82c4762a1bSJed Brown 83c4762a1bSJed Brown for (j = 0; j < t; j++, mult *= d) { 84c4762a1bSJed Brown PetscReal x; 85c4762a1bSJed Brown PetscInt64 slot; 86c4762a1bSJed Brown 875f80ce2aSJacob Faibussowitsch CHKERRQ(PetscRandomGetValueReal(random,&x)); 88c4762a1bSJed Brown /* worried about precision here */ 89c4762a1bSJed Brown slot = (PetscInt64) (x * d); 90c4762a1bSJed Brown bin += mult * slot; 91c4762a1bSJed Brown } 922c71b3e2SJacob Faibussowitsch PetscCheckFalse(bin >= k,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Generated point in bin %" PetscInt64_FMT ", but only %" PetscInt64_FMT " bins",bin,k); 93c4762a1bSJed Brown X[i] = bin; 94c4762a1bSJed Brown } 95c4762a1bSJed Brown 96c4762a1bSJed Brown qsort(X, n, sizeof(PetscInt64), PetscInt64Compare); 97c4762a1bSJed Brown for (i = 0; i < n - 1; i++) { 98c4762a1bSJed Brown X[i] = X[i + 1] - X[i]; 99c4762a1bSJed Brown } 100c4762a1bSJed Brown qsort(X, n - 1, sizeof(PetscInt64), PetscInt64Compare); 101c4762a1bSJed Brown for (i = 0, Y = 0; i < n - 2; i++) { 102c4762a1bSJed Brown Y += (X[i + 1] == X[i]); 103c4762a1bSJed Brown } 104c4762a1bSJed Brown 1055f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Allreduce(MPI_IN_PLACE, &Y, 1, MPIU_INT, MPI_SUM, MPI_COMM_WORLD)); 1065f80ce2aSJacob Faibussowitsch CHKERRQ(PoissonTailProbability(N*lambda,Y,&p)); 1075f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT " total collisions counted: that many or more should occur with probabilty %g.\n",Y,(double)p)); 108c4762a1bSJed Brown 1095f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFree(X)); 1105f80ce2aSJacob Faibussowitsch CHKERRQ(PetscRandomDestroy(&random)); 111*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscFinalize()); 112*b122ec5aSJacob Faibussowitsch return 0; 113c4762a1bSJed Brown } 114c4762a1bSJed Brown 115c4762a1bSJed Brown /*TEST 116c4762a1bSJed Brown 117c4762a1bSJed Brown test: 118c4762a1bSJed Brown args: -t 4 -log2d 7 -log2n 10 119c4762a1bSJed Brown 120c4762a1bSJed Brown TEST*/ 121