1 /*
2 * $Id: diehard_parking_lot.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 *========================================================================
10 * This is the Diehard Parking Lot test, rewritten from the description
11 * in tests.txt on George Marsaglia's diehard site.
12 *
13 * :: THIS IS A PARKING LOT TEST ::
14 * :: In a square of side 100, randomly "park" a car---a circle of ::
15 * :: radius 1. Then try to park a 2nd, a 3rd, and so on, each ::
16 * :: time parking "by ear". That is, if an attempt to park a car ::
17 * :: causes a crash with one already parked, try again at a new ::
18 * :: random location. (To avoid path problems, consider parking ::
19 * :: helicopters rather than cars.) Each attempt leads to either ::
20 * :: a crash or a success, the latter followed by an increment to ::
21 * :: the list of cars already parked. If we plot n: the number of ::
22 * :: attempts, versus k:: the number successfully parked, we get a::
23 * :: curve that should be similar to those provided by a perfect ::
24 * :: random number generator. Theory for the behavior of such a ::
25 * :: random curve seems beyond reach, and as graphics displays are ::
26 * :: not available for this battery of tests, a simple characteriz ::
27 * :: ation of the random experiment is used: k, the number of cars ::
28 * :: successfully parked after n=12,000 attempts. Simulation shows ::
29 * :: that k should average 3523 with sigma 21.9 and is very close ::
30 * :: to normally distributed. Thus (k-3523)/21.9 should be a st- ::
31 * :: andard normal variable, which, converted to a uniform varia- ::
32 * :: ble, provides input to a KSTEST based on a sample of 10. ::
33 *
34 * Comments
35 *
36 * First, the description above is incorrect in two regards.
37 * As seen in the original code, the test measures
38 * overlap of SQUARES of radius one, a thing I only observed after
39 * actually programming circles as described (which are
40 * also easy, although a bit more expensive to evaluate crashes
41 * for). Circles produce an altogether different mean and are
42 * probably a bit more sensitive to 2d striping at arbitrary angles.
43 *
44 * Note that I strongly suspect that this test is basically
45 * equivalent to Knuth's better conceived hyperplane test, which
46 * measures aggregation of N dimensional sets of "coordinates" in
47 * hyperplanes. To put it another way, if something fails a
48 * hyperplane test in 2d, it will certainly fail this test as well.
49 * If something fails this test, I'd bet serious money that it
50 * is because of aggregation of points on hyperplanes although
51 * there MAY be other failure patterns as well.
52 *
53 * Finally, note that the probability that any given k is
54 * obtained for a normal random distribution is just determined
55 * from the erf() -- this is just an Xtest().
56 *
57 * As always, we will increase the number of tsamples and hopefully improve
58 * the resolution of the test. However, it should be carefully noted
59 * that modern random number generators can almost certainly add many
60 * decimal places to the simulation value used in this test. In other
61 * words, test failure at higher resolution can be INVERTED -- it can
62 * indicate the relative failure of the generators used to produce the
63 * earlier result! This is really a subject for future research...
64 *========================================================================
65 */
66
67
68 #include <dieharder/libdieharder.h>
69
70 typedef struct {
71 double x;
72 double y;
73 } Cars;
74
diehard_parking_lot(Test ** test,int irun)75 int diehard_parking_lot(Test **test, int irun)
76 {
77
78 /*
79 * This is the most that could under any circumstances be parked.
80 */
81 Cars parked[12000];
82 uint k,n,i,crashed;
83 double xtry,ytry;
84 Xtest ptest;
85
86 /*
87 * for display only. 0 means "ignored".
88 */
89 test[0]->ntuple = 0;
90 test[0]->tsamples = 12000;
91
92 /*
93 * ptest.x = (double) k
94 * ptest.y = 3523.0
95 * ptest.sigma = 21.9
96 * This will generate ptest->pvalue when Xtest(ptest) is called
97 */
98 ptest.y = 3523.0;
99 ptest.sigma = 21.9;
100
101 /*
102 * Clear the parking lot the fast way.
103 */
104 memset(parked,0,12000*sizeof(Cars));
105
106 /*
107 * Park a single car to have something to avoid and count it.
108 */
109 parked[0].x = 100.0*gsl_rng_uniform(rng);
110 parked[0].y = 100.0*gsl_rng_uniform(rng);
111 k = 1;
112
113
114 /*
115 * This is now a really simple test. Park them cars! We try to park
116 * 12000 times, and increment k (the number successfully parked) on
117 * successes. We brute force the crash test.
118 */
119 for(n=1;n<12000;n++){
120 xtry = 100.0*gsl_rng_uniform(rng);
121 ytry = 100.0*gsl_rng_uniform(rng);
122 crashed = 0;
123 for(i=0;i<k;i++){
124 /*
125 * We do this REASONABLY efficiently. As soon as we know we didn't
126 * crash we move on until we learn that we crashed, trying to skip
127 * arithmetic. Once we've crashed, we break out of the loop.
128 * Uncrashed survivors join the parked list.
129 */
130 if( (fabs(parked[i].x - xtry) <= 1.0) && (fabs(parked[i].y - ytry) <= 1.0)){
131 crashed = 1; /* We crashed! */
132 break; /* So quit the loop here */
133 }
134 }
135 /*
136 * Save uncrashed helicopter coordinates.
137 */
138 if(crashed == 0){
139 parked[k].x = xtry;
140 parked[k].y = ytry;
141 crashed = 0;
142 k++;
143 }
144 }
145
146 ptest.x = (double)k;
147 Xtest_eval(&ptest);
148 test[0]->pvalues[irun] = ptest.pvalue;
149
150 MYDEBUG(D_DIEHARD_PARKING_LOT) {
151 printf("# diehard_parking_lot(): test[0]->pvalues[%u] = %10.5f\n",irun,test[0]->pvalues[irun]);
152 }
153
154 return(0);
155
156 }
157
158