1 /*
2 * $Id: diehard_3dsphere.c 231 2006-08-22 16:18:05Z rgb $
3 *
4 * See copyright in copyright.h and the accompanying file COPYING
5 */
6
7 /*
8 *========================================================================
9 * This is the Diehard 3d spheres test, rewritten from the description
10 * in tests.txt on George Marsaglia's diehard site.
11 *
12 * Basically, we choose 4000 points in a cube of side 1000. Compute the
13 * smallest nearest neighbor distance (radius R of the smallest sphere
14 * that doesn't overlap any neighboring point). R^3 is exponentially
15 * distributed with an empirical exponential distribution with mean 30.
16 * Thus p = 1.0 - exp(-R^3/30.0) should be a uniform distribution. Run
17 * a KS test on a vector of independent samples of this entire test to
18 * find out.
19 *========================================================================
20 */
21
22
23 #include <dieharder/libdieharder.h>
24
25 #define POINTS_3D 4000
26 #define DIM_3D 3
27
28 typedef struct {
29 double x[DIM_3D];
30 } C3_3D;
31
diehard_3dsphere(Test ** test,int irun)32 int diehard_3dsphere(Test **test, int irun)
33 {
34
35 int j,k;
36 C3_3D *c3;
37 double r1,r2,r3,rmin,r3min;
38 double xdelta,ydelta,zdelta;
39
40 /*
41 * for display only. Test dimension is 3, of course.
42 */
43 test[0]->ntuple = 3;
44
45 r3min = 0;
46
47 /*
48 * This one should be pretty straightforward. Generate a vector
49 * of three random coordinates in the range 0-1000 (check the
50 * diehard code to see what "in" a 1000^3 cube means, but I'm assuming
51 * real number coordinates greater than 0 and less than 1000). Do
52 * a simple double loop through to float the smallest separation out.
53 * Generate p, save in a sample vector. Apply KS test.
54 */
55 c3 = (C3_3D *)malloc(POINTS_3D*sizeof(C3_3D));
56
57 rmin = 2000.0;
58 for(j=0;j<POINTS_3D;j++){
59 /*
60 * Generate a new point in the cube.
61 */
62 for(k=0;k<DIM_3D;k++) c3[j].x[k] = 1000.0*gsl_rng_uniform_pos(rng);
63 if(verbose == D_DIEHARD_3DSPHERE || verbose == D_ALL){
64 printf("%d: (%8.2f,%8.2f,%8.2f)\n",j,c3[j].x[0],c3[j].x[1],c3[j].x[2]);
65 }
66
67 /*
68 * Now compute the distance between the new point and all previously
69 * picked points.
70 */
71 for(k=j-1;k>=0;k--){
72 xdelta = c3[j].x[0]-c3[k].x[0];
73 ydelta = c3[j].x[1]-c3[k].x[1];
74 zdelta = c3[j].x[2]-c3[k].x[2];
75 r2 = xdelta*xdelta + ydelta*ydelta + zdelta*zdelta;
76 r1 = sqrt(r2);
77 r3 = r2*r1;
78 if(verbose == D_DIEHARD_3DSPHERE || verbose == D_ALL){
79 printf("%d-%d: |(%6.2f,%6.2f,%6.2f)| = r1 = %f rmin = %f, \n",
80 j,k,xdelta,ydelta,zdelta,r1,rmin);
81 }
82 if(r1<rmin) {
83 rmin = r1;
84 r3min = r3;
85 }
86 }
87 }
88
89 MYDEBUG(D_DIEHARD_3DSPHERE) {
90 printf("Found rmin = %f (r^3 = %f)\n",rmin,r3min);
91 }
92 test[0]->pvalues[irun] = 1.0 - exp(-r3min/30.0);
93
94 MYDEBUG(D_DIEHARD_3DSPHERE) {
95 printf("# diehard_3dsphere(): test[0]->pvalues[%u] = %10.5f\n",irun,test[0]->pvalues[irun]);
96 }
97
98 nullfree(c3);
99
100 return(0);
101
102 }
103
104