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

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