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
|
/*
Apply various randomness tests to a stream of bytes
by John Walker -- September 1996
http://www.fourmilab.ch/
*/
#define _ISOC99_SOURCE
#include <math.h>
#define FALSE 0
#define TRUE 1
static int binary = FALSE; /* Treat input as a bitstream */
static unsigned long 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 */
#define rt_log2 log2
#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 = 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;
}
|