math/mpreduce.h: Missing include files.

[u/mdw/catacomb] / rand / fibrand.c

/* -*-c-*-
 *
 * Fibonacci generator
 *
 * (c) 1999 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.
 */

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

#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

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

#include "fibrand.h"
#include "grand.h"
#include "lcrand.h"

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

/* --- @fibrand_step@ --- *
 *
 * Arguments:	@fibrand *f@ = pointer to Fibonacci generator context
 *
 * Returns:	Next output from generator.
 *
 * Use:		Steps the generator.  Returns
 *		%$x_{i - 24} + x_{i - 55} \bmod 2^{32}$%.
 */

uint32 fibrand_step(fibrand *f)
{
  unsigned i = f->i;
  unsigned j = i + (FIB_SZ - FIB_TAP);
  uint32 x;
  if (j >= FIB_SZ)
    j -= FIB_SZ;
  x = f->x[i] = U32(f->x[i] + f->x[j]);
  i++;
  if (i >= FIB_SZ)
    i = 0;
  f->i = i;
  return (x);
}

/* --- @fibrand_seed@ --- *
 *
 * Arguments:	@fibrand *f@ = pointer to Fibonacci generator context
 *		@grand *r@ = random number generator to extract words from
 *
 * Returns:	---
 *
 * Use:		Initializes a Fibonacci generator using word outputs from the
 *		given random number source @r@.
 */

void fibrand_seed(fibrand *f, grand *r)
{
  int i;
  unsigned p = 0;

  for (i = 0; i < FIB_SZ; i++)
    p |= f->x[i] = r->ops->word(r);
  if (!(p & 1)) {
    i = r->ops->range(r, FIB_SZ);
    f->x[i] |= 1;
  }
  f->i = 0;
}

/* --- @fibrand_lcseed@ --- *
 *
 * Arguments:	@fibrand *f@ = pointer to Fibonacci generator context
 *		@uint32 seed@ = seed value
 *
 * Returns:	---
 *
 * Use:		Initializes a Fibonacci generator using outputs from the
 *		@lcrand@ generator seeded from @seed@.  This is faster than
 *		using a generic @lcrand@-based generator and @fibrand_rseed@
 *		because it uses raw outputs rather than uniformly distributed
 *		32-bit words.
 */

void fibrand_lcseed(fibrand *f, uint32 seed)
{
  int i;
  unsigned p = 0;

  for (i = 0; i < FIB_SZ; i++)
    p |= f->x[i] = seed = lcrand(seed);
  if (!(p & 1)) {
    i = lcrand_range(&seed, FIB_SZ);
    f->x[i] |= 1;
  }
  f->i = 0;
}

/* --- @fibrand_range@ --- *
 *
 * Arguments:	@fibrand *f@ = pointer to Fibonacci generator context
 *		@uint32 m@ = limit
 *
 * Returns:	A uniformly distributed pseudorandom integer in the interval
 *		%$[0, m)$%.
 */

uint32 fibrand_range(fibrand *f, uint32 m)
{
  uint32 r = 0xffffffff - (0xffffffff % m);
  uint32 x;

  /* --- Now generate numbers until a good one comes along --- */

  do x = fibrand_step(f); while (x >= r);
  return (x % m);
}

/*----- Generic interface -------------------------------------------------*/

typedef struct gctx {
  grand r;
  fibrand f;
} gctx;

static void gdestroy(grand *r)
{
  gctx *g = (gctx *)r;
  DESTROY(g);
}

static int gmisc(grand *r, unsigned op, ...)
{
  gctx *g = (gctx *)r;
  va_list ap;
  int rc = 0;
  va_start(ap, op);

  switch (op) {
    case GRAND_CHECK:
      switch (va_arg(ap, unsigned)) {
	case GRAND_CHECK:
	case GRAND_SEEDINT:
	case GRAND_SEEDUINT32:
	case GRAND_SEEDRAND:
	  rc = 1;
	  break;
	default:
	  rc = 0;
	  break;
      }
      break;
    case GRAND_SEEDINT:
      fibrand_lcseed(&g->f, va_arg(ap, unsigned));
      break;
    case GRAND_SEEDUINT32:
      fibrand_lcseed(&g->f, va_arg(ap, uint32));
      break;
    case GRAND_SEEDRAND:
      fibrand_seed(&g->f, va_arg(ap, grand *));
      break;
    default:
      GRAND_BADOP;
      break;
  }

  va_end(ap);
  return (rc);
}

static octet gbyte(grand *r)
{
  gctx *g = (gctx *)r;
  return (U8(fibrand_step(&g->f)));
}

static uint32 gword(grand *r)
{
  gctx *g = (gctx *)r;
  return (fibrand_step(&g->f));
}

static uint32 grange(grand *r, uint32 l)
{
  gctx *g = (gctx *)r;
  return (fibrand_range(&g->f, l));
}

static void gfill(grand *r, void *p, size_t sz)
{
  gctx *g = (gctx *)r;
  octet *q = p;
  while (sz) {
    *q++ = U8(fibrand_step(&g->f));
    sz--;
  }
}

static const grand_ops gops = {
  "fibrand",
  0, 0,
  gmisc, gdestroy,
  gword, gbyte, gword, grange, gfill
};

/* --- @fibrand_create@ --- *
 *
 * Arguments:	@uint32 seed@ = initial seed
 *
 * Returns:	Pointer to a generic generator.
 *
 * Use:		Constructs a generic generator interface over a Fibonacci
 *		generator.  The generator is seeded using @fibrand_lcseed@.
 */

grand *fibrand_create(uint32 seed)
{
  gctx *g = CREATE(gctx);
  g->r.ops = &gops;
  fibrand_lcseed(&g->f, seed);
  return (&g->r);
}

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