File: randtest.c

package info (click to toggle)
ent 1.2debian-2
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 260 kB
  • sloc: ansic: 410; makefile: 41
file content (179 lines) | stat: -rw-r--r-- 4,646 bytes parent folder | download | duplicates (3)
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;
}