1*c4762a1bSJed Brown 2*c4762a1bSJed Brown static char help[] = "Run Birthday Spacing Tests for PetscRandom.\n\n"; 3*c4762a1bSJed Brown 4*c4762a1bSJed Brown #include <petscsys.h> 5*c4762a1bSJed Brown #include <petscviewer.h> 6*c4762a1bSJed Brown 7*c4762a1bSJed Brown /* L'Ecuyer & Simard, 2001. 8*c4762a1bSJed Brown * "On the performance of birthday spacings tests with certain families of random number generators" 9*c4762a1bSJed Brown * https://doi.org/10.1016/S0378-4754(00)00253-6 10*c4762a1bSJed Brown */ 11*c4762a1bSJed Brown 12*c4762a1bSJed Brown static int PetscInt64Compare(const void *a, const void *b) 13*c4762a1bSJed Brown { 14*c4762a1bSJed Brown PetscInt64 A = *((const PetscInt64 *)a); 15*c4762a1bSJed Brown PetscInt64 B = *((const PetscInt64 *)b); 16*c4762a1bSJed Brown if (A < B) return -1; 17*c4762a1bSJed Brown if (A == B) return 0; 18*c4762a1bSJed Brown return 1; 19*c4762a1bSJed Brown } 20*c4762a1bSJed Brown 21*c4762a1bSJed Brown static PetscErrorCode PoissonTailProbability(PetscReal lambda, PetscInt Y, PetscReal *prob) 22*c4762a1bSJed Brown { 23*c4762a1bSJed Brown PetscReal p = 1.; 24*c4762a1bSJed Brown PetscReal logLambda; 25*c4762a1bSJed Brown PetscInt i; 26*c4762a1bSJed Brown PetscReal logFact = 0.; 27*c4762a1bSJed Brown 28*c4762a1bSJed Brown PetscFunctionBegin; 29*c4762a1bSJed Brown logLambda = PetscLogReal(lambda); 30*c4762a1bSJed Brown for (i = 0; i < Y; i++) { 31*c4762a1bSJed Brown PetscReal exponent = -lambda + i * logLambda - logFact; 32*c4762a1bSJed Brown 33*c4762a1bSJed Brown logFact += PetscLogReal(i+1); 34*c4762a1bSJed Brown 35*c4762a1bSJed Brown p -= PetscExpReal(exponent); 36*c4762a1bSJed Brown } 37*c4762a1bSJed Brown *prob = p; 38*c4762a1bSJed Brown PetscFunctionReturn(0); 39*c4762a1bSJed Brown } 40*c4762a1bSJed Brown 41*c4762a1bSJed Brown int main(int argc, char **argv) 42*c4762a1bSJed Brown { 43*c4762a1bSJed Brown PetscMPIInt size; 44*c4762a1bSJed Brown PetscInt log2d, log2n, t, N, Y; 45*c4762a1bSJed Brown PetscInt64 d, k, *X; 46*c4762a1bSJed Brown size_t n, i; 47*c4762a1bSJed Brown PetscReal lambda, p; 48*c4762a1bSJed Brown PetscRandom random; 49*c4762a1bSJed Brown PetscErrorCode ierr; 50*c4762a1bSJed Brown 51*c4762a1bSJed Brown ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 52*c4762a1bSJed Brown 53*c4762a1bSJed Brown t = 8; 54*c4762a1bSJed Brown log2d = 7; 55*c4762a1bSJed Brown log2n = 20; 56*c4762a1bSJed Brown ierr = PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Birthday spacing test parameters", "PetscRandom");CHKERRQ(ierr); 57*c4762a1bSJed Brown ierr = PetscOptionsInt("-t", "t, the dimension of the sample space", "ex3.c", t, &t, NULL);CHKERRQ(ierr); 58*c4762a1bSJed Brown ierr = PetscOptionsInt("-log2d", "The log of d, the number of bins per direction", "ex3.c", log2d, &log2d, NULL);CHKERRQ(ierr); 59*c4762a1bSJed Brown ierr = PetscOptionsInt("-log2n", "The log of n, the number of samples per process", "ex3.c", log2n, &log2n, NULL);CHKERRQ(ierr); 60*c4762a1bSJed Brown ierr = PetscOptionsEnd();CHKERRQ(ierr); 61*c4762a1bSJed Brown 62*c4762a1bSJed Brown 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*c4762a1bSJed Brown d = ((PetscInt64) 1) << log2d; 64*c4762a1bSJed Brown k = ((PetscInt64) 1) << (log2d * t); 65*c4762a1bSJed Brown 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*c4762a1bSJed Brown n = ((size_t) 1) << log2n; 67*c4762a1bSJed Brown ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 68*c4762a1bSJed Brown N = size; 69*c4762a1bSJed Brown lambda = PetscPowRealInt(2.,(3 * log2n - (2 + log2d * t))); 70*c4762a1bSJed Brown 71*c4762a1bSJed Brown 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*c4762a1bSJed Brown ierr = PetscPrintf(PETSC_COMM_WORLD,"Expected spacing collisions per process %f (%f total).\n", (double) lambda, N * lambda);CHKERRQ(ierr); 73*c4762a1bSJed Brown ierr = PetscRandomCreate(PETSC_COMM_WORLD,&random);CHKERRQ(ierr); 74*c4762a1bSJed Brown ierr = PetscRandomSetFromOptions(random);CHKERRQ(ierr); 75*c4762a1bSJed Brown ierr = PetscRandomSetInterval(random,0.0,1.0);CHKERRQ(ierr); 76*c4762a1bSJed Brown ierr = PetscRandomView(random,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 77*c4762a1bSJed Brown ierr = PetscMalloc1(n,&X);CHKERRQ(ierr); 78*c4762a1bSJed Brown for (i = 0; i < n; i++) { 79*c4762a1bSJed Brown PetscInt j; 80*c4762a1bSJed Brown PetscInt64 bin = 0; 81*c4762a1bSJed Brown PetscInt64 mult = 1; 82*c4762a1bSJed Brown 83*c4762a1bSJed Brown for (j = 0; j < t; j++, mult *= d) { 84*c4762a1bSJed Brown PetscReal x; 85*c4762a1bSJed Brown PetscInt64 slot; 86*c4762a1bSJed Brown 87*c4762a1bSJed Brown ierr = PetscRandomGetValueReal(random,&x);CHKERRQ(ierr); 88*c4762a1bSJed Brown /* worried about precision here */ 89*c4762a1bSJed Brown slot = (PetscInt64) (x * d); 90*c4762a1bSJed Brown bin += mult * slot; 91*c4762a1bSJed Brown } 92*c4762a1bSJed Brown 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*c4762a1bSJed Brown X[i] = bin; 94*c4762a1bSJed Brown } 95*c4762a1bSJed Brown 96*c4762a1bSJed Brown qsort(X, n, sizeof(PetscInt64), PetscInt64Compare); 97*c4762a1bSJed Brown for (i = 0; i < n - 1; i++) { 98*c4762a1bSJed Brown X[i] = X[i + 1] - X[i]; 99*c4762a1bSJed Brown } 100*c4762a1bSJed Brown qsort(X, n - 1, sizeof(PetscInt64), PetscInt64Compare); 101*c4762a1bSJed Brown for (i = 0, Y = 0; i < n - 2; i++) { 102*c4762a1bSJed Brown Y += (X[i + 1] == X[i]); 103*c4762a1bSJed Brown } 104*c4762a1bSJed Brown 105*c4762a1bSJed Brown ierr = MPI_Allreduce(MPI_IN_PLACE, &Y, 1, MPIU_INT, MPI_SUM, MPI_COMM_WORLD);CHKERRQ(ierr); 106*c4762a1bSJed Brown ierr = PoissonTailProbability(N*lambda,Y,&p);CHKERRQ(ierr); 107*c4762a1bSJed Brown 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*c4762a1bSJed Brown 109*c4762a1bSJed Brown ierr = PetscFree(X);CHKERRQ(ierr); 110*c4762a1bSJed Brown ierr = PetscRandomDestroy(&random);CHKERRQ(ierr); 111*c4762a1bSJed Brown ierr = PetscFinalize(); 112*c4762a1bSJed Brown return ierr; 113*c4762a1bSJed Brown } 114*c4762a1bSJed Brown 115*c4762a1bSJed Brown /*TEST 116*c4762a1bSJed Brown 117*c4762a1bSJed Brown test: 118*c4762a1bSJed Brown args: -t 4 -log2d 7 -log2n 10 119*c4762a1bSJed Brown 120*c4762a1bSJed Brown TEST*/ 121