1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
|
/* This code was taken from http://www.fourmilab.ch/random/ where it states that:
This software is in the public domain. Permission to use, copy, modify, and distribute
this software and its documentation for any purpose and without fee is hereby granted,
without any conditions or restrictions. This software is provided “as is” without
express or implied warranty. */
/*
Apply various randomness tests to a stream of bytes
by John Walker -- September 1996
http://www.fourmilab.ch/
*/
#include <math.h>
#define FALSE 0
#define TRUE 1
#define log2of10 3.32192809488736234787
static int binary = FALSE; /* Treat input as a bitstream */
static long ccount[256], /* Bins to count occurrences of values */
totalc = 0; /* Total bytes counted */
static double prob[256]; /* Probabilities per bin for entropy */
/* RT_LOG2 -- Calculate log to the base 2 */
static double rt_log2(double x)
{
return log2of10 * log10(x);
}
#define MONTEN 6 /* Bytes used as Monte Carlo
co-ordinates. This should be no more
bits than the mantissa of your
"double" floating point type. */
static int mp, sccfirst;
static unsigned int monte[MONTEN];
static long inmont, mcount;
static double cexp, incirc, montex, montey, montepi,
scc, sccun, sccu0, scclast, scct1, scct2, scct3,
ent, chisq, datasum;
/* RT_INIT -- Initialise random test counters. */
void rt_init(int binmode)
{
int i;
binary = binmode; /* Set binary / byte mode */
/* Initialise for calculations */
ent = 0.0; /* Clear entropy accumulator */
chisq = 0.0; /* Clear Chi-Square */
datasum = 0.0; /* Clear sum of bytes for arithmetic mean */
mp = 0; /* Reset Monte Carlo accumulator pointer */
mcount = 0; /* Clear Monte Carlo tries */
inmont = 0; /* Clear Monte Carlo inside count */
incirc = 65535.0 * 65535.0;/* In-circle distance for Monte Carlo */
sccfirst = TRUE; /* Mark first time for serial correlation */
scct1 = scct2 = scct3 = 0.0; /* Clear serial correlation terms */
incirc = pow(pow(256.0, (double) (MONTEN / 2)) - 1, 2.0);
for (i = 0; i < 256; i++) {
ccount[i] = 0;
}
totalc = 0;
}
/* RT_ADD -- Add one or more bytes to accumulation. */
void rt_add(void *buf, int bufl)
{
unsigned char *bp = (unsigned char *)buf;
int oc, c, bean;
while (bean = 0, (bufl-- > 0)) {
oc = *bp++;
do {
if (binary) {
c = !!(oc & 0x80);
} else {
c = oc;
}
ccount[c]++; /* Update counter for this bin */
totalc++;
/* Update inside / outside circle counts for Monte Carlo
computation of PI */
if (bean == 0) {
monte[mp++] = oc; /* Save character for Monte Carlo */
if (mp >= MONTEN) { /* Calculate every MONTEN character */
int mj;
mp = 0;
mcount++;
montex = montey = 0;
for (mj = 0; mj < MONTEN / 2; mj++) {
montex = (montex * 256.0) + monte[mj];
montey = (montey * 256.0) + monte[(MONTEN / 2) + mj];
}
if ((montex * montex + montey * montey) <= incirc) {
inmont++;
}
}
}
/* Update calculation of serial correlation coefficient */
sccun = c;
if (sccfirst) {
sccfirst = FALSE;
scclast = 0;
sccu0 = sccun;
} else {
scct1 = scct1 + scclast * sccun;
}
scct2 = scct2 + sccun;
scct3 = scct3 + (sccun * sccun);
scclast = sccun;
oc <<= 1;
} while (binary && (++bean < 8));
}
}
/* RT_END -- Complete calculation and return results. */
void rt_end(double *r_ent, double *r_chisq, double *r_mean,
double *r_montepicalc, double *r_scc)
{
int i;
/* Complete calculation of serial correlation coefficient */
scct1 = scct1 + scclast * sccu0;
scct2 = scct2 * scct2;
scc = totalc * scct3 - scct2;
if (scc == 0.0) {
scc = -100000;
} else {
scc = (totalc * scct1 - scct2) / scc;
}
/* Scan bins and calculate probability for each bin and
Chi-Square distribution. The probability will be reused
in the entropy calculation below. While we're at it,
we sum of all the data which will be used to compute the
mean. */
cexp = totalc / (binary ? 2.0 : 256.0); /* Expected count per bin */
for (i = 0; i < (binary ? 2 : 256); i++) {
double a = ccount[i] - cexp;;
prob[i] = ((double) ccount[i]) / totalc;
chisq += (a * a) / cexp;
datasum += ((double) i) * ccount[i];
}
/* Calculate entropy */
for (i = 0; i < (binary ? 2 : 256); i++) {
if (prob[i] > 0.0) {
ent += prob[i] * rt_log2(1 / prob[i]);
}
}
/* Calculate Monte Carlo value for PI from percentage of hits
within the circle */
montepi = 4.0 * (((double) inmont) / mcount);
/* Return results through arguments */
*r_ent = ent;
*r_chisq = chisq;
*r_mean = datasum / totalc;
*r_montepicalc = montepi;
*r_scc = scc;
}
|
