1 /*
2 *=========================================================
3 * DAB Fill Tree 2 Test
4 * Bit version of Fill Tree test.
5 * This test fills small binary trees of fixed depth with
6 * "visited" markers. When a marker cannot be placed, the
7 * current count of markers in the tree and the position
8 * that the marker would have been inserted, if it hadn't
9 * already been marked.
10 *
11 * For each bit in the RNG input, the test takes a step
12 * right (for a zero) or left (for a one) in the tree.
13 * If the node hasn't been marked, it is marked, and the
14 * path restarts. Otherwise, the test continues with the
15 * next bit.
16 *
17 * The test returns two p-values. The first is a Pearson
18 * chi-sq test against the expected values (which were
19 * estimated empirically. The second is a Pearson chi-sq
20 * test for a uniform distribution of the positions at
21 * which the insert failed.
22 *
23 * Because of the target data for the first p-value,
24 * ntuple must be kept at the default (128).
25 */
26 #include <dieharder/libdieharder.h>
27
28 typedef unsigned char uchar;
29
30 static double targetData1[32] __attribute__((unused)) = { // size=32, generated from 3e9 samples
31 0.00000000000e+00, 0.00000000000e+00, 0.00000000000e+00, 0.00000000000e+00,
32 9.76265666667e-04, 3.90648133333e-03, 9.42791500000e-03, 1.77898240000e-02,
33 2.88606903333e-02, 4.21206876667e-02, 5.67006123333e-02, 7.13000270000e-02,
34 8.43831060000e-02, 9.43500813333e-02, 9.97647353333e-02, 9.97074473333e-02,
35 9.39825660000e-02, 8.32745740000e-02, 6.90631193333e-02, 5.33001223333e-02,
36 3.80133193333e-02, 2.48386790000e-02, 1.47058170000e-02, 7.79203100000e-03,
37 3.63280466667e-03, 1.45833733333e-03, 4.88868000000e-04, 1.31773666667e-04,
38 2.63546666667e-05, 3.51800000000e-06, 2.42333333333e-07, 0.00000000000e+00
39 };
40
41 static double targetData2[64] __attribute__((unused)) = { // size=64, generated from 3e9 samples
42 0.00000000000e+00, 0.00000000000e+00, 0.00000000000e+00, 0.00000000000e+00,
43 0.00000000000e+00, 3.03990000000e-05, 1.52768666667e-04, 4.47074666667e-04,
44 1.00459133333e-03, 1.91267566667e-03, 3.25090066667e-03, 5.08490633333e-03,
45 7.45162400000e-03, 1.03865720000e-02, 1.38770320000e-02, 1.78957393333e-02,
46 2.23788223333e-02, 2.72281453333e-02, 3.23125090000e-02, 3.74760433333e-02,
47 4.25407143333e-02, 4.72809176667e-02, 5.15021953333e-02, 5.49909926667e-02,
48 5.75592136667e-02, 5.90709556667e-02, 5.94040870000e-02, 5.85191160000e-02,
49 5.64374923333e-02, 5.32542393333e-02, 4.91296690000e-02, 4.42833420000e-02,
50 3.89462263333e-02, 3.33966756667e-02, 2.78853806667e-02, 2.26466276667e-02,
51 1.78641310000e-02, 1.36666050000e-02, 1.01212106667e-02, 7.24765200000e-03,
52 5.00267066667e-03, 3.32686800000e-03, 2.12376733333e-03, 1.29971400000e-03,
53 7.59987666667e-04, 4.23040000000e-04, 2.23585333333e-04, 1.11777000000e-04,
54 5.27836666667e-05, 2.33956666667e-05, 9.67033333333e-06, 3.63066666667e-06,
55 1.29500000000e-06, 4.04666666667e-07, 1.24666666667e-07, 2.80000000000e-08,
56 8.33333333333e-09, 1.33333333333e-09, 0.00000000000e+00, 0.00000000000e+00,
57 0.00000000000e+00, 0.00000000000e+00, 0.00000000000e+00, 0.00000000000e+00
58 };
59
60 static double targetData[128] = { // size=128, generated from 6e9 samples
61 0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,
62 0.00000000000e+00,0.00000000000e+00,4.77166666667e-07,2.85516666667e-06,
63 9.85200000000e-06,2.55708333333e-05,5.54055000000e-05,1.06338000000e-04,
64 1.86151333333e-04,3.02864333333e-04,4.66803833333e-04,6.86296166667e-04,
65 9.71489833333e-04,1.33154316667e-03,1.77454416667e-03,2.31186450000e-03,
66 2.94724100000e-03,3.68750350000e-03,4.53733233333e-03,5.50155300000e-03,
67 6.58318550000e-03,7.77896366667e-03,9.08643266667e-03,1.05029766667e-02,
68 1.20238296667e-02,1.36316733333e-02,1.53215870000e-02,1.70715931667e-02,
69 1.88714935000e-02,2.06986750000e-02,2.25274171667e-02,2.43387205000e-02,
70 2.61018481667e-02,2.77880516667e-02,2.93792388333e-02,3.08369918333e-02,
71 3.21362530000e-02,3.32571148333e-02,3.41773491667e-02,3.48649186667e-02,
72 3.53142736667e-02,3.55033806667e-02,3.54345836667e-02,3.50962926667e-02,
73 3.44943555000e-02,3.36465076667e-02,3.25496588333e-02,3.12430565000e-02,
74 2.97450508333e-02,2.80761470000e-02,2.62786038333e-02,2.43732936667e-02,
75 2.24096238333e-02,2.04087746667e-02,1.84149658333e-02,1.64520138333e-02,
76 1.45553360000e-02,1.27460330000e-02,1.10448845000e-02,9.47002933333e-03,
77 8.03035333333e-03,6.73145016667e-03,5.58088483333e-03,4.56875066667e-03,
78 3.69710383333e-03,2.95267200000e-03,2.32891533333e-03,1.81311533333e-03,
79 1.39093383333e-03,1.05322383333e-03,7.86386833333e-04,5.78384666667e-04,
80 4.18726333333e-04,2.98206500000e-04,2.09200000000e-04,1.44533833333e-04,
81 9.79670000000e-05,6.52683333333e-05,4.27531666667e-05,2.73350000000e-05,
82 1.72115000000e-05,1.06390000000e-05,6.38166666667e-06,3.79950000000e-06,
83 2.21116666667e-06,1.25500000000e-06,6.95833333333e-07,3.83500000000e-07,
84 1.98000000000e-07,9.91666666667e-08,4.73333333333e-08,2.66666666667e-08,
85 1.13333333333e-08,5.50000000000e-09,2.66666666667e-09,1.00000000000e-09,
86 6.66666666667e-10,3.33333333333e-10,0.00000000000e+00,0.00000000000e+00,
87 0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,
88 0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,
89 0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,
90 0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,
91 0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,
92 0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,0.00000000000e+00,
93 };
94
95 inline int insertBit(uint x, uchar *array, uint *i, uint *d);
96
dab_filltree2(Test ** test,int irun)97 int dab_filltree2(Test **test, int irun) {
98 int size = (ntuple == 0) ? 128 : ntuple;
99 uint target = sizeof(targetData)/sizeof(double);
100 int startVal = (size / 2) - 1;
101 uchar *array = (uchar *) malloc(sizeof(*array) * size);
102 double *counts, *expected;
103 int i, j;
104 uint x;
105 uint start = 0;
106 uint end = 0;
107 double *positionCounts;
108 uint bitCount;
109
110 test[0]->ntuple = 0;
111 test[1]->ntuple = 1;
112
113 counts = (double *) malloc(sizeof(double) * target);
114 expected = (double *) malloc(sizeof(double) * target);
115
116 memset(counts, 0, sizeof(double) * target);
117
118 positionCounts = (double *) malloc(sizeof(double) * size/2);
119 memset(positionCounts, 0, sizeof(double) * size/2);
120
121 /* Calculate expected counts. */
122 for (i = 0; i < target; i++) {
123 expected[i] = targetData[i] * test[0]->tsamples;
124 if (expected[i] < 4) {
125 if (end == 0) start = i;
126 } else if (expected[i] > 4) end = i;
127 }
128 start++;
129
130
131 x = gsl_rng_get(rng);
132 bitCount = rmax_bits;
133 for (j = 0; j < test[0]->tsamples; j++) {
134 int ret;
135 memset(array, 0, sizeof(*array) * size);
136 i = 0;
137 do { /* While new markers can be aded to this tree.... */
138 uint index = startVal;
139 uint d = (startVal + 1) / 2;
140 if (i > size * 2) {
141 test[0]->pvalues[irun] = 0;
142 return(0);
143 }
144 do { /* While this path has not yet found a blank node to mark.... */
145 ret = insertBit(x & 0x01, array, &index, &d); /* Keep going */
146 x >>= 1;
147 if (--bitCount == 0) {
148 x = gsl_rng_get(rng);
149 bitCount = rmax_bits;
150 }
151 } while (ret == -2); /* End of path. */
152
153 i++;
154 } while (ret == -1); /* Couldn't insert marker; end of this tree. */
155 positionCounts[ret/2]++;
156
157 counts[i-1]++;
158 }
159
160 /* First p-value is calculated based on the targetData array. */
161 test[0]->pvalues[irun] = chisq_pearson(counts + start, expected + start, end - start);
162
163 /* Second p-value is calculated against a uniform distribution. */
164 for (i = 0; i < size/2; i++) expected[i] = test[0]->tsamples/(size/2);
165 test[1]->pvalues[irun] = chisq_pearson(positionCounts, expected, size/2);
166
167
168 nullfree(positionCounts);
169 nullfree(expected);
170 nullfree(counts);
171 nullfree(array);
172
173 return(0);
174 }
175
176 /*
177 * Insert a bit into the tree, represented by an array.
178 * A value of one is marked; zero is unmarked.
179 * The function returns -2 is still on the path.
180 * The function returns -1 if the path ends by marking a node.
181 * The function returns >= 0 if the path went too deep; the
182 * returned value is the last position of the path.
183 */
insertBit(uint x,uchar * array,uint * i,uint * d)184 inline int insertBit(uint x, uchar *array, uint *i, uint *d) {
185 if (x != 0) {
186 *i += *d;
187 } else {
188 *i -= *d;
189 }
190 *d /= 2;
191
192 if (array[*i] == 0) {
193 array[*i] = 1;
194 return -1;
195 }
196 if (*d == 0) return *i;
197 else return -2;
198 }
199
200