/*
	random.h
	created: June 1, 2010
	
	author: gachaz
	
	function: store all functions to generate random numbers
*/


#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>


#define EPS FLT_MIN 
/*
in the original numer it was
#define EPS 1.2e-7 
*/

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define RNMX (1.0-EPS)


/***
	Random number generation
****/

static unsigned long idum_ran1;   /* the key value of ran1() */



/*
	Generate a random number between 0.00 and 1.00
	Should be called with a negative long the first time,
	then should not be reseeded.
	From Numerical recipies - ran1.c
*/
float ran1(long *idum_ran1){

	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	float temp;

	if (*idum_ran1 <= 0 || !iy) {
		if (-(*idum_ran1) < 1) *idum_ran1=1;
		else *idum_ran1 = -(*idum_ran1);
		for (j=NTAB+7;j>=0;j--) {
			k=(*idum_ran1)/IQ;
			*idum_ran1=IA*(*idum_ran1-k*IQ)-IR*k;
			if (*idum_ran1 < 0) *idum_ran1 += IM;
			if (j < NTAB) iv[j] = *idum_ran1;
		}
		iy=iv[0];
	}
	k=(*idum_ran1)/IQ;
	*idum_ran1=IA*(*idum_ran1-k*IQ)-IR*k;
	if (*idum_ran1 < 0) *idum_ran1 += IM;
	j=iy/NDIV;
	iy=iv[j];
	iv[j] = *idum_ran1;
	if ((temp=AM*iy) > RNMX) return RNMX;
	else return temp;
}
/*
	Seed the function ran1() with a negative long
	Should be used only ONCE before using ran1()
*/
void seed_ran1( long seed ){

	extern unsigned long idum_ran1;            /* argument of ran1() */
	long tmp;                                  /* a tmp value fom seeding ran1();  */
	
	
	if(seed == -1)
		idum_ran1 = (unsigned long)time(NULL);     /* time_t is a long in macosx */
	else
		idum_ran1 =(unsigned long)seed;
		
	tmp = -idum_ran1;
	ran1(&tmp);                                /* seed the function */
}


/* 
	uniform_dev return a number [0,1[
*/
double uniform_dev( void ){

	extern unsigned long idum_ran1;            /* argument of ran1() */

	double dum= ran1((long *)&idum_ran1);
	
	if(dum == 0.0)
		fprintf(stderr, "BUG: ran1() can also return 0\n"),exit(1);

	if(dum == 1.0)dum=0;
	
	return dum;
	
}
/*
	uniform_int is [min , max]
*/
long uniform_int(long min, long max) 
{ 
	double dum;
 	long diff=max-min+1;
	
	dum = uniform_dev();
	 
	return (min + (int)floor(dum*(double)diff)); 
}
/*
	Get two different long in a range
*/
void uniform_2DiffInt(long min, long max, long *int1, long *int2){

	*int1 = uniform_int( min,  max);
	do{
		*int2 = uniform_int( min,  max);
	}while(*int1 == *int2);

}


/*
	Get an array of different long in a range with an exception (i.e. an excluded value);
*/
long *uniform_ArrayDiffInt(long min, long max, long n, long exception ){

	long *array,       /* returned array */
	     *array_t;;    /* a tmp array that is used to avoid redrawing the same int */

	long x,            /* a counter */
	     r;            /* random int */
	
	char e=0;
	
	
	if( exception >= min && exception <= max  )
		e=1;
	

	if(n>max-min+1-e){
		fprintf(stderr, "uniform_ArrayDiffInt: array size should be less (or equal) than the number of different values\n");
		return NULL;
	}
	
	array = (long *)malloc(  (size_t) n*sizeof(long) );
	array_t = (long *)malloc(  (size_t) (max-min+1-e)*sizeof(long) );
	if(! array || !array_t )fprintf(stderr, "uniform_ArrayDiffInt: cannot allocate array, bye\n"), exit(3);

	/*
		Fill array_t with all possible values
	*/
	e=0;
	for(x=0; x < max-min+1; x++){
	
		if( x == exception )
			e=1;
		else
			array_t[x-e] = min+x;
	}

	/*
		draw from array_t
	*/
	x=0;
	do{
		r = uniform_int( 0,  max-min-e-x );
		array[ x ] = array_t[ r ];
		array_t[ r ] = array_t[ max-min-e-x  ];
		x++;
	}while( x < n );


	free(array_t);

	return array;
}


float *uniform_array( long n  ){


	float * array, *p;
	
	array = (float *) malloc( (size_t) n*sizeof(float) );
	if(!array)fprintf(stderr,"uniform_array: cannot allocate array, bye\n"), exit(3);
	
	
	for( p=array; p-array<n; p++)
		*p=uniform_dev(  );

	return array;

}




/*
	inspired from numrec
	normal_CR_dev is mean 0 and var 1
*/
float normal_CR_dev( )
{
	static int iset=0;
	static float gset;
	float fac,rsq,v1,v2;


	if  (iset == 0) {
	
		do {
			v1=2.0*uniform_dev()-1.0;
			v2=2.0*uniform_dev()-1.0;
			rsq=v1*v1+v2*v2;
		} while (rsq >= 1.0 || rsq == 0.0);
		
		fac=sqrt(-2.0*log(rsq)/rsq);
		gset=v1*fac;
		iset=1;
		
		return v2*fac;
	}
	else
	{
		iset=0;
		return gset;
	}
}

float normal_dev( double mean, double stdev ){
	return normal_CR_dev( )*stdev + mean;
}


float *normal_array( int n, double mean, double stdev ){


	float * array, *p;
	
	array = (float *) malloc( (size_t) n*sizeof(float) );
	if(!array)fprintf(stderr,"normal_array: cannot allocate array, bye\n"), exit(3);
	
	
	for( p=array; p-array<n; p++)
		*p=normal_dev( mean, stdev );

	return array;

}