3 * $Id: maurer.c,v 1.1 2000/06/17 11:29:49 mdw Exp $
5 * Maurer's universal statistical test
7 * (c) 2000 Straylight/Edgeware
10 /*----- Licensing notice --------------------------------------------------*
12 * This file is part of Catacomb.
14 * Catacomb is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU Library General Public License as
16 * published by the Free Software Foundation; either version 2 of the
17 * License, or (at your option) any later version.
19 * Catacomb is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU Library General Public License for more details.
24 * You should have received a copy of the GNU Library General Public
25 * License along with Catacomb; if not, write to the Free
26 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
30 /*----- Revision history --------------------------------------------------*
33 * Revision 1.1 2000/06/17 11:29:49 mdw
34 * Maurer's universal statistical test.
38 /*----- Header files ------------------------------------------------------*/
45 #include <mLib/alloc.h>
46 #include <mLib/bits.h>
50 /*----- Data structures ---------------------------------------------------*/
52 typedef struct bitscan
{
58 /*----- Main code ---------------------------------------------------------*/
62 * Arguments: @bitscan *b@ = pointer to a bitscanner structure
63 * @unsigned l@ = number of bits wanted
65 * Returns: An integer representing the next @l@ bits of the stream.
67 * Use: Fetches the next @l@ bits of a data stream. The bitscanner
68 * has no way of knowing whether it's run out of data, so make
69 * sure that doesn't happen.
72 static unsigned more(bitscan
*b
, unsigned l
)
77 b
->a
|= U8(*b
->p
++) << b
->b
;
89 * Arguments: @const octet *p@ = pointer to a buffer of data
90 * @size_t sz@ = size of the buffer
91 * @unsigned l@ = number of bits to take at a time
93 * Returns: The statistic %$Z_u = (X_u - \mu)/\sigma$%, which should be
94 * normally distributed with a mean of zero and standard
97 * Use: Performs Maurer's universal statistical test on a buffer of
98 * data. You must ensure that the buffer size @sz@ is larger
99 * than %$11 l \cdot 2^l$% bits long, and ideally much larger.
101 * Note that, under Unix systems, this function requires the
105 double maurer(const octet
*p
, size_t sz
, unsigned l
)
107 size_t *t
; /* Table of occurrences */
108 size_t q
, k
; /* Block counts */
109 size_t i
; /* Loop counter */
110 bitscan b
; /* Bitscanning context */
111 double x
; /* Statistic accumulator */
112 double mu
, c
, sigma
; /* Other useful things */
114 /* --- Table for means and variances --- */
116 struct { double mu
, var_1
; } vtab
[] = {
117 { 0.7326495, 0.690 },
118 { 1.5374383, 1.338 },
119 { 2.4016068, 1.901 },
120 { 3.3112247, 2.358 },
121 { 4.2534266, 2.705 },
122 { 5.2177052, 2.945 },
123 { 6.1962507, 3.125 },
124 { 7.1836656, 3.238 },
125 { 8.1764248, 3.311 },
126 { 9.1723243, 3.356 },
127 { 10.170032 , 3.384 },
128 { 11.168765 , 3.401 },
129 { 12.168070 , 3.410 },
130 { 13.167693 , 3.416 },
131 { 14.167488 , 3.419 },
132 { 15.167379 , 3.421 }
135 assert(((void)"bit chunk size out of range", 0 < l
&& l
<= 16));
137 /* --- Initialize everything --- */
139 t
= xmalloc(sizeof(size_t) << l
);
140 for (i
= 0; i
< (1 << l
); i
++)
149 assert(((void)"buffer too small for test", k
> q
&& k
- q
> (1 << l
)));
151 /* --- Set up the table --- */
153 for (i
= 0; i
< q
; i
++) {
154 unsigned a
= more(&b
, l
);
158 /* --- Now compute the statistic --- */
162 unsigned a
= more(&b
, l
);
167 /* --- Twiddle at the end --- *
169 * This produces %$X_u$%, which should be normally distributed with a mean
170 * and variance we can compute.
176 /* --- Now translate this into %$Z_u$% distributed as %$N(0, 1)$% --- */
179 c
= 0.7 - 0.8/l
+ (1.6 + 12.8/l
) * pow(k
, -4.0/l
);
180 sigma
= sqrt(c
* c
* vtab
[l
- 1].var_1
/ k
);
184 /* --- Done (phew!) --- */
190 /*----- That's all, folks -------------------------------------------------*/