/* $Id: xsh_model_r250.c,v 1.6 2011-12-02 14:15:28 amodigli Exp $
 *
 *Not sure about the copyright stuff here!
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

/*
 * $Author: amodigli $
 * $Date: 2011-12-02 14:15:28 $
 * $Revision: 1.6 $
 * $Name: not supported by cvs2svn $
 */
/* r250.c	the r250 uniform random number algorithm

		Kirkpatrick, S., and E. Stoll, 1981; "A Very Fast
		Shift-Register Sequence Random Number Generator",
		Journal of Computational Physics, V.40

		also:

		see W.L. Maier, DDJ May 1991



*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
/*---------------------------------------------------------------------------*/
/**
 * @addtogroup xsh_model    xsh_model_r250
 *
 * Physical model random generator
 * See man-page for details.
 */
/*---------------------------------------------------------------------------*/
/**@{*/
/*-----------------------------------------------------------------------------
   								Includes
 -----------------------------------------------------------------------------*/

#include <cpl.h>

#include "xsh_model_kernel.h"

/*----------------------------------------------------------------------------*/

//static char rcsid[] = "@(#)r250.c	1.2 15:50:31 11/21/94   EFC";


#include <limits.h>

#include "xsh_model_r250.h"

/* set the following if you trust rand(), otherwise the minimal standard
   generator is used
*/
/* #define TRUST_RAND */


#ifndef TRUST_RAND
#include "xsh_model_randlcg.h"
#endif

/* defines to allow for 16 or 32 bit integers */
#define BITS 31
#define WORD_BIT 32

#if WORD_BIT == 32
#ifndef BITS
#define BITS	32
#endif
#else
#ifndef BITS
#define BITS    16
#endif
#endif

#if BITS == 31
#define MSB          0x40000000L
#define ALL_BITS     0x7fffffffL
#define HALF_RANGE   0x20000000L
#define STEP         7
#endif

#if BITS == 32
#define MSB          0x80000000L
#define ALL_BITS     0xffffffffL
#define HALF_RANGE   0x40000000L
#define STEP         7
#endif

#if BITS == 16
#define MSB         0x8000
#define ALL_BITS    0xffff
#define HALF_RANGE  0x4000
#define STEP        11
#endif

static unsigned int r250_buffer[ 250 ];
static int r250_index;

#ifdef NO_PROTO
void xsh_r250_init(sd)
int seed;
#else
void xsh_r250_init(int sd)
#endif
{
	int j, k;
	unsigned int mask, msb;

#ifdef TRUST_RAND        

#if BITS == 32 || BITS == 31       
	srand48( sd );
#else
	srand( sd );
#endif	


#else
	xsh_set_seed( sd );
#endif
	
	r250_index = 0;
	for (j = 0; j < 250; j++)      /* fill r250 buffer with BITS-1 bit values */
#ifdef TRUST_RAND
#if BITS == 32 || BITS == 31
		r250_buffer[j] = (unsigned int)lrand48();
#else
		r250_buffer[j] = rand();
#endif
#else
		r250_buffer[j] = xsh_randlcg();
#endif


	for (j = 0; j < 250; j++)	/* set some MSBs to 1 */
#ifdef TRUST_RAND
		if ( rand() > HALF_RANGE )
			r250_buffer[j] |= MSB;
#else
		if ( xsh_randlcg() > HALF_RANGE )
			r250_buffer[j] |= MSB;
#endif


	msb = MSB;	        /* turn on diagonal bit */
	mask = ALL_BITS;	/* turn off the leftmost bits */

	for (j = 0; j < BITS; j++)
	{
		k = STEP * j + 3;	/* select a word to operate on */
		r250_buffer[k] &= mask; /* turn off bits left of the diagonal */
		r250_buffer[k] |= msb;	/* turn on the diagonal bit */
		mask >>= 1;
		msb  >>= 1;
	}

}

unsigned int r250(void)		/* returns a random unsigned integer */
{
	register int	j;
	register unsigned int new_rand;

	if ( r250_index >= 147 )
		j = r250_index - 147;	/* wrap pointer around */
	else
		j = r250_index + 103;

	new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
	r250_buffer[ r250_index ] = new_rand;

	if ( r250_index >= 249 )	/* increment pointer for next time */
		r250_index = 0;
	else
		r250_index++;

	return new_rand;

}


double xsh_dr250(void)		/* returns a random double in range 0..1 */
{
	register int	j;
	register unsigned int new_rand;

	if ( r250_index >= 147 )
		j = r250_index - 147;	/* wrap pointer around */
	else
		j = r250_index + 103;

	new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
	r250_buffer[ r250_index ] = new_rand;

	if ( r250_index >= 249 )	/* increment pointer for next time */
		r250_index = 0;
	else
		r250_index++;

	return (double)new_rand / ALL_BITS;

}

#ifdef MAIN

/* test driver	prints out either NMR_RAND values or a histogram	*/

#include <stdio.h>

#define NMR_RAND	5000
#define MAX_BINS	500

#ifdef NO_PROTO
void main(argc, argv)
int argc;
char **argv;
#else
void main(int argc, char **argv)
#endif
{
	int j,k,nmr_bins,seed;
	int bins[MAX_BINS];
	double randm, bin_inc;
	double bin_limit[MAX_BINS];

	if ( argc != 3 )
	{
		printf("Usage -- %s nmr_bins seed\n", argv[0]);
		abort();
	}

	nmr_bins = atoi( argv[1] );
	if ( nmr_bins > MAX_BINS )
	{
		printf("ERROR -- maximum number of bins is %d\n", MAX_BINS);
		abort();
	}

	seed = atoi( argv[2] );

	xsh_r250_init( seed );

	if ( nmr_bins < 1 )	/* just print out the numbers */
	{
		for (j = 0; j < NMR_RAND; j++)
			printf("%f\n", xsh_dr250() );
		abort();
	}
	
	bin_inc = 1.0 / nmr_bins;
	for (j = 0; j < nmr_bins; j++)	/* initialize bins to zero */
	{
		bins[j] = 0;
		bin_limit[j] = (j + 1) * bin_inc;
	}

	bin_limit[nmr_bins-1] = 1.0e7;	/* make sure all others are in last bin */

	for (j = 0; j < NMR_RAND; j++)
	{
		randm = r250() / (double)ALL_BITS;
		for (k = 0; k < nmr_bins; k++)
			if ( randm < bin_limit[k] )
			{
				bins[k]++;
				break;
			}
	}


	for (j = 0; j < nmr_bins; j++)
		printf("%d\n", bins[j]);
	
}

#endif
/**@}*/