chiark - git - mdw - catacomb/blob - maurer.c

   1 /* -*-c-*-
   2  *
   3  * $Id: maurer.c,v 1.4 2004/04/08 01:36:15 mdw Exp $
   4  *
   5  * Maurer's universal statistical test
   6  *
   7  * (c) 2000 Straylight/Edgeware
   8  */
   9
  10 /*----- Licensing notice --------------------------------------------------*
  11  *
  12  * This file is part of Catacomb.
  13  *
  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.
  18  *
  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.
  23  *
  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,
  27  * MA 02111-1307, USA.
  28  */
  29
  30 /*----- Header files ------------------------------------------------------*/
  31
  32 #include <assert.h>
  33 #include <math.h>
  34 #include <stdio.h>
  35 #include <stdlib.h>
  36
  37 #include <mLib/alloc.h>
  38 #include <mLib/bits.h>
  39
  40 #include "maurer.h"
  41
  42 /*----- Data structures ---------------------------------------------------*/
  43
  44 typedef struct bitscan {
  45   const octet *p;
  46   unsigned b;
  47   uint32 a;
  48 } bitscan;
  49
  50 /*----- Main code ---------------------------------------------------------*/
  51
  52 /* --- @maurer_init@ --- *
  53  *
  54  * Arguments:   @maurer_ctx *m@ = pointer to a context to initialize
  55  *              @unsigned l@ = number of bits to take at a time
  56  *
  57  * Returns:     Minimum recommended amount of data to test.
  58  *
  59  * Use:         Initializes a Maurer testing context.  Note that @l@ values
  60  *              larger than 16 are not supported, and values less than 6 can
  61  *              give bizarre results.
  62  */
  63
  64 unsigned long maurer_init(maurer_ctx *m, unsigned l)
  65 {
  66   size_t i;
  67
  68   assert(((void)"(maurer_init): chunks larger than 16 bits not supported",
  69           0 < l && l <= 16));
  70
  71   m->l = l;
  72   m->a = 0;
  73   m->b = 0;
  74   m->n = 0;
  75   m->x = 0;
  76
  77   m->t = xmalloc(sizeof(unsigned long) << l);
  78   for (i = 0; i < (1 << l); i++)
  79     m->t[i] = 0;
  80   return (((1000ul << l) * l + 7)/8);
  81 }
  82
  83 /* --- @bits@ --- *
  84  *
  85  * Arguments:   @maurer_ctx *m@ = pointer to testing context
  86  *              @const octet **p@ = address of a buffer pointer
  87  *              @const octet *q@ = limit of the buffer pointer
  88  *              @unsigned *x@ = address to store next chunk
  89  *
  90  * Returns:     Nonzero if some more bits arrived.
  91  *
  92  * Use:         Reads some bits from a buffer.
  93  */
  94
  95 static int bits(maurer_ctx *m, const octet **p, const octet *q, unsigned *x)
  96 {
  97   while (m->b < m->l) {
  98     if (*p >= q)
  99       return (0);
 100     m->a |= U8(*(*p)++) << m->b;
 101     m->b += 8;
 102   }
 103   *x = m->a & ((1 << m->l) - 1);
 104   m->a >>= m->l;
 105   m->b -= m->l;
 106   m->n++;
 107   return (1);
 108 }
 109
 110 /* --- @maurer_test@ --- *
 111  *
 112  * Arguments:   @maurer_ctx *m@ = pointer to a Maurer testing context
 113  *              @const void *buf@ = pointer to a buffer of data
 114  *              @size_t sz@ = size of the buffer
 115  *
 116  * Returns:     ---
 117  *
 118  * Use:         Scans a buffer of data and updates the testing context.
 119  */
 120
 121 void maurer_test(maurer_ctx *m, const void *buf, size_t sz)
 122 {
 123   const octet *p = buf, *l = p + sz;
 124   unsigned long q = 10 << m->l;
 125   unsigned x;
 126
 127   /* --- Initialize the table --- */
 128
 129   while (m->n < q) {
 130     if (!bits(m, &p, l, &x))
 131       return;
 132     m->t[x] = m->n;
 133   }
 134
 135   /* --- Update the statistic --- */
 136
 137   while (bits(m, &p, l, &x)) {
 138     m->x += log(m->n - m->t[x]);
 139     m->t[x] = m->n;
 140   }
 141 }
 142
 143 /* --- @maurer_done@ --- *
 144  *
 145  * Arguments:   @maurer_ctx *m@ = pointer to a Maurer testing context
 146  *
 147  * Returns:     The statistic %$Z_u = (X_u - \mu)/\sigma$%, which should be
 148  *              normally distributed with a mean of 0 and variance of 1.
 149  *
 150  * Use:         Returns the result of Maurer's universal statistical test.
 151  *              Also frees the memory allocated during initialization.
 152  *
 153  *              If insufficient data was received, the value @HUGE_VAL@ is
 154  *              returned.
 155  */
 156
 157 double maurer_done(maurer_ctx *m)
 158 {
 159   double x, mu, c, sigma;
 160   unsigned long q = 10 << m->l, k;
 161
 162   /* --- Table for means and variances --- */
 163
 164   struct { double mu, var_1; } vtab[] = {
 165     {  0.7326495, 0.690 },
 166     {  1.5374383, 1.338 },
 167     {  2.4016068, 1.901 },
 168     {  3.3112247, 2.358 },
 169     {  4.2534266, 2.705 },
 170     {  5.2177052, 2.945 },
 171     {  6.1962507, 3.125 },
 172     {  7.1836656, 3.238 },
 173     {  8.1764248, 3.311 },
 174     {  9.1723243, 3.356 },
 175     { 10.170032 , 3.384 },
 176     { 11.168765 , 3.401 },
 177     { 12.168070 , 3.410 },
 178     { 13.167693 , 3.416 },
 179     { 14.167488 , 3.419 },
 180     { 15.167379 , 3.421 }
 181   };
 182
 183   /* --- Check for bogus requests --- */
 184
 185   if (m->n < q) {
 186     x = HUGE_VAL;
 187     goto done;
 188   }
 189
 190   /* --- Do the maths --- *
 191    *
 192    * This produces %$X_u$%, which should be normally distributed with a mean
 193    * and variance we can compute.
 194    */
 195
 196   k = m->n - q;
 197   x = m->x/(k * log(2));
 198
 199   /* --- Now translate this into %$Z_u$% distributed as %$N(0, 1)$% --- */
 200
 201   mu = vtab[m->l - 1].mu;
 202   c = 0.7 - 0.8/m->l + (1.6 + 12.8/m->l) * pow(k, -4.0/m->l);
 203   sigma = sqrt(c * c * vtab[m->l - 1].var_1 / k);
 204   x = (x - mu)/sigma;
 205
 206   /* --- Done (phew!) --- */
 207
 208 done:
 209   xfree(m->t);
 210   return (x);
 211 }
 212
 213 /* --- @maurer@ --- *
 214  *
 215  * Arguments:   @const octet *p@ = pointer to a buffer of data
 216  *              @size_t sz@ = size of the buffer
 217  *              @unsigned l@ = number of bits to take at a time
 218  *
 219  * Returns:     The statistic %$Z_u$%.
 220  *
 221  * Use:         Simple interface to Maurer's universal statistical test.  For
 222  *              large %$l$% you should use the more complicated restartable
 223  *              interface.
 224  */
 225
 226 double maurer(const octet *p, size_t sz, unsigned l)
 227 {
 228   maurer_ctx m;
 229
 230   maurer_init(&m, l);
 231   maurer_test(&m, p, sz);
 232   return (maurer_done(&m));
 233 }
 234
 235 /*----- That's all, folks -------------------------------------------------*/