Renamed from `rsa-decrypt', since the name was no longer appropriate.

[u/mdw/catacomb] / maurer.c

/* -*-c-*-
 *
 * $Id: maurer.c,v 1.1 2000/06/17 11:29:49 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.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 ---------------------------------------------------------*/

/* --- @more@ --- *
 *
 * Arguments:	@bitscan *b@ = pointer to a bitscanner structure
 *		@unsigned l@ = number of bits wanted
 *
 * Returns:	An integer representing the next @l@ bits of the stream.
 *
 * Use:		Fetches the next @l@ bits of a data stream.  The bitscanner
 *		has no way of knowing whether it's run out of data, so make
 *		sure that doesn't happen.
 */

static unsigned more(bitscan *b, unsigned l)
{
  unsigned x, m;

  while (b->b < l) {
    b->a |= U8(*b->p++) << b->b;
    b->b += 8;
  }
  m = (1 << l) - 1;
  x = b->a & m;
  b->a >>= l;
  b->b -= l;
  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 = (X_u - \mu)/\sigma$%, which should be
 *		normally distributed with a mean of zero and standard
 *		deviation of 1.
 *
 * Use:		Performs Maurer's universal statistical test on a buffer of
 *		data.  You must ensure that the buffer size @sz@ is larger
 *		than %$11 l \cdot 2^l$% bits long, and ideally much larger. 
 *
 *		Note that, under Unix systems, this function requires the
 *		maths library.
 */

double maurer(const octet *p, size_t sz, unsigned l)
{
  size_t *t;				/* Table of occurrences */
  size_t q, k;				/* Block counts */
  size_t i;				/* Loop counter */
  bitscan b;				/* Bitscanning context */
  double x;				/* Statistic accumulator */
  double mu, c, sigma;			/* Other useful things */

  /* --- 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 }
  };

  assert(((void)"bit chunk size out of range", 0 < l && l <= 16));

  /* --- Initialize everything --- */

  t = xmalloc(sizeof(size_t) << l);
  for (i = 0; i < (1 << l); i++)
    t[i] = 0;

  b.p = p;
  b.a = 0;
  b.b = 0;

  k = (sz * 8)/l;
  q = 10 << l;
  assert(((void)"buffer too small for test", k > q && k - q > (1 << l)));

  /* --- Set up the table --- */

  for (i = 0; i < q; i++) {
    unsigned a = more(&b, l);
    t[a] = i;
  }

  /* --- Now compute the statistic --- */

  x = 0;
  for (; i < k; i++) {
    unsigned a = more(&b, l);
    x += log(i - t[a]);
    t[a] = i;
  }

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

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

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

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

  x = (x - mu)/sigma;

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

  xfree(t);
  return (x);
}

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