Force subkeys to be sorted in structured keys.

[u/mdw/catacomb] / maurer.c

/* -*-c-*-
 *
 * $Id: maurer.c,v 1.3 2000/08/16 17:56:59 mdw Exp $
 *
 * Maurer's universal statistical test
 *
 * (c) 2000 Straylight/Edgeware
 */

/*----- Licensing notice --------------------------------------------------* 
 *
 * This file is part of Catacomb.
 *
 * Catacomb is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Library General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 * 
 * Catacomb is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Library General Public License for more details.
 * 
 * You should have received a copy of the GNU Library General Public
 * License along with Catacomb; if not, write to the Free
 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 * MA 02111-1307, USA.
 */

/*----- Revision history --------------------------------------------------* 
 *
 * $Log: maurer.c,v $
 * Revision 1.3  2000/08/16 17:56:59  mdw
 * (more): Remove spurious function.
 *
 * Revision 1.2  2000/08/11 21:34:59  mdw
 * New restartable interface to Maurer testing.
 *
 * Revision 1.1  2000/06/17 11:29:49  mdw
 * Maurer's universal statistical test.
 *
 */

/*----- Header files ------------------------------------------------------*/

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#include <mLib/alloc.h>
#include <mLib/bits.h>

#include "maurer.h"

/*----- Data structures ---------------------------------------------------*/

typedef struct bitscan {
  const octet *p;
  unsigned b;
  uint32 a;
} bitscan;

/*----- Main code ---------------------------------------------------------*/

/* --- @maurer_init@ --- *
 *
 * Arguments:	@maurer_ctx *m@ = pointer to a context to initialize
 *		@unsigned l@ = number of bits to take at a time
 *
 * Returns:	Minimum recommended amount of data to test.
 *
 * Use:		Initializes a Maurer testing context.  Note that @l@ values
 *		larger than 16 are not supported, and values less than 6 can
 *		give bizarre results.
 */

unsigned long maurer_init(maurer_ctx *m, unsigned l)
{
  size_t i;

  assert(((void)"(maurer_init): chunks larger than 16 bits not supported",
	  0 < l && l <= 16));

  m->l = l;
  m->a = 0;
  m->b = 0;
  m->n = 0;
  m->x = 0;

  m->t = xmalloc(sizeof(unsigned long) << l);
  for (i = 0; i < (1 << l); i++)
    m->t[i] = 0;
  return (((1000ul << l) * l + 7)/8);
}

/* --- @bits@ --- *
 *
 * Arguments:	@maurer_ctx *m@ = pointer to testing context
 *		@const octet **p@ = address of a buffer pointer
 *		@const octet *q@ = limit of the buffer pointer
 *		@unsigned *x@ = address to store next chunk
 *
 * Returns:	Nonzero if some more bits arrived.
 *
 * Use:		Reads some bits from a buffer.
 */

static int bits(maurer_ctx *m, const octet **p, const octet *q, unsigned *x)
{
  while (m->b < m->l) {
    if (*p >= q)
      return (0);
    m->a |= U8(*(*p)++) << m->b;
    m->b += 8;
  }
  *x = m->a & ((1 << m->l) - 1);
  m->a >>= m->l;
  m->b -= m->l;
  m->n++;
  return (1);
}

/* --- @maurer_test@ --- *
 *
 * Arguments:	@maurer_ctx *m@ = pointer to a Maurer testing context
 *		@const void *buf@ = pointer to a buffer of data
 *		@size_t sz@ = size of the buffer
 *
 * Returns:	---
 *
 * Use:		Scans a buffer of data and updates the testing context.
 */

void maurer_test(maurer_ctx *m, const void *buf, size_t sz)
{
  const octet *p = buf, *l = p + sz;
  unsigned long q = 10 << m->l;
  unsigned x;

  /* --- Initialize the table --- */

  while (m->n < q) {
    if (!bits(m, &p, l, &x))
      return;
    m->t[x] = m->n;
  }

  /* --- Update the statistic --- */

  while (bits(m, &p, l, &x)) {
    m->x += log(m->n - m->t[x]);
    m->t[x] = m->n;
  }
}

/* --- @maurer_done@ --- *
 *
 * Arguments:	@maurer_ctx *m@ = pointer to a Maurer testing context
 *
 * Returns:	The statistic %$Z_u = (X_u - \mu)/\sigma$%, which should be
 *		normally distributed with a mean of 0 and variance of 1.
 *
 * Use:		Returns the result of Maurer's universal statistical test.
 *		Also frees the memory allocated during initialization.
 *
 *		If insufficient data was received, the value @HUGE_VAL@ is
 *		returned.
 */

double maurer_done(maurer_ctx *m)
{
  double x, mu, c, sigma;
  unsigned long q = 10 << m->l, k;

  /* --- Table for means and variances --- */

  struct { double mu, var_1; } vtab[] = {
    {  0.7326495, 0.690 },
    {  1.5374383, 1.338 },
    {  2.4016068, 1.901 },
    {  3.3112247, 2.358 },
    {  4.2534266, 2.705 },
    {  5.2177052, 2.945 },
    {  6.1962507, 3.125 },
    {  7.1836656, 3.238 },
    {  8.1764248, 3.311 },
    {  9.1723243, 3.356 },
    { 10.170032 , 3.384 },
    { 11.168765 , 3.401 },
    { 12.168070 , 3.410 },
    { 13.167693 , 3.416 },
    { 14.167488 , 3.419 },
    { 15.167379 , 3.421 }
  };

  /* --- Check for bogus requests --- */

  if (m->n < q) {
    x = HUGE_VAL;
    goto done;
  }

  /* --- Do the maths --- *
   *
   * This produces %$X_u$%, which should be normally distributed with a mean
   * and variance we can compute.
   */

  k = m->n - q;
  x = m->x/(k * log(2));

  /* --- Now translate this into %$Z_u$% distributed as %$N(0, 1)$% --- */

  mu = vtab[m->l - 1].mu;
  c = 0.7 - 0.8/m->l + (1.6 + 12.8/m->l) * pow(k, -4.0/m->l);
  sigma = sqrt(c * c * vtab[m->l - 1].var_1 / k);
  x = (x - mu)/sigma;

  /* --- Done (phew!) --- */

done:
  xfree(m->t);
  return (x);
}

/* --- @maurer@ --- *
 *
 * Arguments:	@const octet *p@ = pointer to a buffer of data
 *		@size_t sz@ = size of the buffer
 *		@unsigned l@ = number of bits to take at a time
 *
 * Returns:	The statistic %$Z_u$%.
 *
 * Use:		Simple interface to Maurer's universal statistical test.  For
 *		large %$l$% you should use the more complicated restartable
 *		interface.
 */

double maurer(const octet *p, size_t sz, unsigned l)
{
  maurer_ctx m;

  maurer_init(&m, l);
  maurer_test(&m, p, sz);
  return (maurer_done(&m));
}

/*----- That's all, folks -------------------------------------------------*/