/* libpbm6.c - pbm utility library part 6
**
** Simple, portable, reasonably robust random number generator.
**
** Copyright (C) 1994 by Jef Poskanzer.
**
** Permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notice appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation.  This software is provided "as is" without express or
** implied warranty.
*/

#include <time.h>
#include "pbm.h"


/* This is a combination of a linear congruential generator and a feedback
** shift register.  Values from the LCG are used to keep a circular buffer
** filled; results are produced by xoring three values from the table.
** The modulus of the LCG must be a power of two for this to produce
** equidistributed results.  This LCG actually uses a modulus that's
** a power of two minus one, but that's close enough.
**
** DO NOT MODIFY, IMPROVE, EXPAND, ENHANCE, OR IN ANY WAY CHANGE this
** generator.  It is used for cryptographic storage of data - if the
** sequence is changed, the data will become unrecoverable.
**
** The linear congruential generator is:
**     Minimal Standard Pseudo-Random Number Generator
**     Author: Fuat C. Baran, Columbia University, 1988
**     Based on code in "Random Number Generators: Good Ones are Hard to Find",
**     by Stephen K. Park and Keith W. Miller in Communications of the ACM,
**     31, 10 (Oct. 1988) pp. 1192-1201.
**
** The feedback shift register is similar to the one described in "Algorithms",
** Robert Sedgewick, 1983, page 38.
*/

#define A 16807L
#define M 2147483647L		/* Mersenne prime 2^31 -1 */
#define Q 127773L		/* M div A (M / A) */
#define R 2836L			/* M mod A (M % A) */
static long value = 1;

#define TABLESIZE 55
#define TAP1 0
#define TAP2 23
#define TAP3 (TABLESIZE-1)
static long table[TABLESIZE];
static int offset;

static long
lcg( void )
    {
    long hi, lo;

    hi = value / Q;
    lo = value % Q;
    value = A * lo - R * hi;
    if ( value <= 0 )
	value += M;
    return value;
    }

void
pm_srandom( unsigned long seed )
    {
    if ( seed == 0 )
	/* Zero doesn't work in this RNG anyway, so we use it as a flag. */
	value = time( (time_t*) 0 ) ^ getpid();
    else
	value = seed;

    for ( offset = 0; offset < TABLESIZE; ++offset )
	table[offset] = lcg();
    }

long
pm_random( void )
    {
    offset = ( offset + 1 ) % TABLESIZE;
    table[offset] = lcg();
    return table[offset] ^			/* TAP1 is zero, optimize */
	   table[( offset + TAP2 ) % TABLESIZE] ^
	   table[( offset + TAP3 ) % TABLESIZE];
    }