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

   1 /* -*-c-*-
   2  *
   3  * $Id: maurer.c,v 1.1 2000/06/17 11:29:49 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 /*----- Revision history --------------------------------------------------*
  31  *
  32  * $Log: maurer.c,v $
  33  * Revision 1.1  2000/06/17 11:29:49  mdw
  34  * Maurer's universal statistical test.
  35  *
  36  */
  37
  38 /*----- Header files ------------------------------------------------------*/
  39
  40 #include <assert.h>
  41 #include <math.h>
  42 #include <stdio.h>
  43 #include <stdlib.h>
  44
  45 #include <mLib/alloc.h>
  46 #include <mLib/bits.h>
  47
  48 #include "maurer.h"
  49
  50 /*----- Data structures ---------------------------------------------------*/
  51
  52 typedef struct bitscan {
  53   const octet *p;
  54   unsigned b;
  55   uint32 a;
  56 } bitscan;
  57
  58 /*----- Main code ---------------------------------------------------------*/
  59
  60 /* --- @more@ --- *
  61  *
  62  * Arguments:   @bitscan *b@ = pointer to a bitscanner structure
  63  *              @unsigned l@ = number of bits wanted
  64  *
  65  * Returns:     An integer representing the next @l@ bits of the stream.
  66  *
  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.
  70  */
  71
  72 static unsigned more(bitscan *b, unsigned l)
  73 {
  74   unsigned x, m;
  75
  76   while (b->b < l) {
  77     b->a |= U8(*b->p++) << b->b;
  78     b->b += 8;
  79   }
  80   m = (1 << l) - 1;
  81   x = b->a & m;
  82   b->a >>= l;
  83   b->b -= l;
  84   return (x);
  85 }
  86
  87 /* --- @maurer@ --- *
  88  *
  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
  92  *
  93  * Returns:     The statistic %$Z_u = (X_u - \mu)/\sigma$%, which should be
  94  *              normally distributed with a mean of zero and standard
  95  *              deviation of 1.
  96  *
  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.
 100  *
 101  *              Note that, under Unix systems, this function requires the
 102  *              maths library.
 103  */
 104
 105 double maurer(const octet *p, size_t sz, unsigned l)
 106 {
 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 */
 113
 114   /* --- Table for means and variances --- */
 115
 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 }
 133   };
 134
 135   assert(((void)"bit chunk size out of range", 0 < l && l <= 16));
 136
 137   /* --- Initialize everything --- */
 138
 139   t = xmalloc(sizeof(size_t) << l);
 140   for (i = 0; i < (1 << l); i++)
 141     t[i] = 0;
 142
 143   b.p = p;
 144   b.a = 0;
 145   b.b = 0;
 146
 147   k = (sz * 8)/l;
 148   q = 10 << l;
 149   assert(((void)"buffer too small for test", k > q && k - q > (1 << l)));
 150
 151   /* --- Set up the table --- */
 152
 153   for (i = 0; i < q; i++) {
 154     unsigned a = more(&b, l);
 155     t[a] = i;
 156   }
 157
 158   /* --- Now compute the statistic --- */
 159
 160   x = 0;
 161   for (; i < k; i++) {
 162     unsigned a = more(&b, l);
 163     x += log(i - t[a]);
 164     t[a] = i;
 165   }
 166
 167   /* --- Twiddle at the end --- *
 168    *
 169    * This produces %$X_u$%, which should be normally distributed with a mean
 170    * and variance we can compute.
 171    */
 172
 173   k -= q;
 174   x /= k * log(2);
 175
 176   /* --- Now translate this into %$Z_u$% distributed as %$N(0, 1)$% --- */
 177
 178   mu = vtab[l - 1].mu;
 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);
 181
 182   x = (x - mu)/sigma;
 183
 184   /* --- Done (phew!) --- */
 185
 186   xfree(t);
 187   return (x);
 188 }
 189
 190 /*----- That's all, folks -------------------------------------------------*/