/**CFile***********************************************************************

  FileName    [random.c]

  PackageName [util]

  Synopsis    [Our own portable random number generator.]

  Author      [Fabio Somenzi <Fabio@Colorado.EDU>]

  Copyright   [Copyright (c) 1994-1996 The Regents of the Univ. of California.
  All rights reserved.

  Permission is hereby granted, without written agreement and without license
  or royalty fees, to use, copy, modify, and distribute this software and its
  documentation for any purpose, provided that the above copyright notice and
  the following two paragraphs appear in all copies of this software.

  IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY FOR
  DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT
  OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF
  CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

  THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES,
  INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
  FITNESS FOR A PARTICULAR PURPOSE.  THE SOFTWARE PROVIDED HEREUNDER IS ON AN
  \"AS IS\" BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATION TO PROVIDE
  MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.]

******************************************************************************/

#include "util.h" 

static char rcsid[] UNUSED = "$Id: random.c,v 1.4 2002/09/10 00:05:33 fabio Exp $";

/*---------------------------------------------------------------------------*/
/* Constant declarations                                                     */
/*---------------------------------------------------------------------------*/

/* Random generator constants. */
#define MODULUS1 2147483563
#define LEQA1 40014
#define LEQQ1 53668
#define LEQR1 12211
#define MODULUS2 2147483399
#define LEQA2 40692
#define LEQQ2 52774
#define LEQR2 3791
#define STAB_SIZE 64
#define STAB_DIV (1 + (MODULUS1 - 1) / STAB_SIZE)

/*---------------------------------------------------------------------------*/
/* Type declarations                                                         */
/*---------------------------------------------------------------------------*/


/*---------------------------------------------------------------------------*/
/* Structure declarations                                                    */
/*---------------------------------------------------------------------------*/


/*---------------------------------------------------------------------------*/
/* Variable declarations                                                     */
/*---------------------------------------------------------------------------*/
static long utilRand = 0;
static long utilRand2;
static long shuffleSelect;
static long shuffleTable[STAB_SIZE];

/*---------------------------------------------------------------------------*/
/* Macro declarations                                                        */
/*---------------------------------------------------------------------------*/


/**AutomaticStart*************************************************************/

/*---------------------------------------------------------------------------*/
/* Static function prototypes                                                */
/*---------------------------------------------------------------------------*/

/**AutomaticEnd***************************************************************/


/*---------------------------------------------------------------------------*/
/* Definition of exported functions                                          */
/*---------------------------------------------------------------------------*/


/**Function********************************************************************

  Synopsis    [Initializer for the portable random number generator.]

  Description [Initializer for the portable number generator based on
  ran2 in "Numerical Recipes in C." The input is the seed for the
  generator. If it is negative, its absolute value is taken as seed.
  If it is 0, then 1 is taken as seed. The initialized sets up the two
  recurrences used to generate a long-period stream, and sets up the
  shuffle table.]

  SideEffects [None]

  SeeAlso     [util_random]

******************************************************************************/
void
util_srandom(long seed)
{
    int i;

    if (seed < 0)       utilRand = -seed;
    else if (seed == 0) utilRand = 1;
    else                utilRand = seed;
    utilRand2 = utilRand;
    /* Load the shuffle table (after 11 warm-ups). */
    for (i = 0; i < STAB_SIZE + 11; i++) {
	long int w;
	w = utilRand / LEQQ1;
	utilRand = LEQA1 * (utilRand - w * LEQQ1) - w * LEQR1;
	utilRand += (utilRand < 0) * MODULUS1;
	shuffleTable[i % STAB_SIZE] = utilRand;
    }
    shuffleSelect = shuffleTable[1 % STAB_SIZE];
} /* end of util_srandom */

/**Function********************************************************************

  Synopsis    [Portable random number generator.]

  Description [Portable number generator based on ran2 from "Numerical
  Recipes in C." It is a long period (> 2 * 10^18) random number generator
  of L'Ecuyer with Bays-Durham shuffle. Returns a long integer uniformly
  distributed between 0 and 2147483561 (inclusive of the endpoint values).
  The random generator can be explicitly initialized by calling
  util_srandom. If no explicit initialization is performed, then the
  seed 1 is assumed.]

  SideEffects []

  SeeAlso     [util_srandom]

******************************************************************************/
long
util_random(void)
{
    int i;	/* index in the shuffle table */
    long int w; /* work variable */

    /* utilRand == 0 if the geneartor has not been initialized yet. */
    if (utilRand == 0) util_srandom(1);

    /* Compute utilRand = (utilRand * LEQA1) % MODULUS1 avoiding
    ** overflows by Schrage's method.
    */
    w          = utilRand / LEQQ1;
    utilRand   = LEQA1 * (utilRand - w * LEQQ1) - w * LEQR1;
    utilRand  += (utilRand < 0) * MODULUS1;

    /* Compute utilRand2 = (utilRand2 * LEQA2) % MODULUS2 avoiding
    ** overflows by Schrage's method.
    */
    w          = utilRand2 / LEQQ2;
    utilRand2  = LEQA2 * (utilRand2 - w * LEQQ2) - w * LEQR2;
    utilRand2 += (utilRand2 < 0) * MODULUS2;

    /* utilRand is shuffled with the Bays-Durham algorithm.
    ** shuffleSelect and utilRand2 are combined to generate the output.
    */

    /* Pick one element from the shuffle table; "i" will be in the range
    ** from 0 to STAB_SIZE-1.
    */
    i = shuffleSelect / STAB_DIV;
    /* Mix the element of the shuffle table with the current iterate of
    ** the second sub-generator, and replace the chosen element of the
    ** shuffle table with the current iterate of the first sub-generator.
    */
    shuffleSelect   = shuffleTable[i] - utilRand2;
    shuffleTable[i] = utilRand;
    shuffleSelect  += (shuffleSelect < 1) * (MODULUS1 - 1);
    /* Since shuffleSelect != 0, and we want to be able to return 0,
    ** here we subtract 1 before returning.
    */
    return(shuffleSelect - 1);

} /* end of util_random */


/*---------------------------------------------------------------------------*/
/* Definition of internal functions                                          */
/*---------------------------------------------------------------------------*/


/*---------------------------------------------------------------------------*/
/* Definition of static functions                                            */
/*---------------------------------------------------------------------------*/