/* Steven Andrews, 5/12/95; modified 8/6/96, converted to #define 11/98 */
/* See documentation called random doc */
/* Copyright 2003 by Steven Andrews.  Permission is granted
   for non-commercial use of and modifications to the code. */

#ifndef __random_h
#define __random_h

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

#define RAND_MAX_30 1073741823
#if RAND_MAX==32767
	#define rand30() ((long)rand()<<15|rand())
#elif RAND_MAX<RAND_MAX_30
	#define rand30() ((rand()&32767L)<<15|rand()&32767L)
#else
	#define rand30() (rand()&RAND_MAX_30)
#endif

#define exprand(a) (-log(((float)rand()+1.0)/(RAND_MAX+1.0))*(a))
#define exprand30(a) (-log(((double)rand30()+1.0)/(RAND_MAX_30+1.0))*(a))
#define unirand(lo,hi) ((float)rand()/RAND_MAX*((hi)-(lo))+(lo))
#define unirand30(lo,hi) ((double)rand30()/RAND_MAX_30*((hi)-(lo))+(lo))
#define signrand() (rand()&1?1:-1)
#define coinrand(p) (rand()<(RAND_MAX+1.0)*(p))
#define coinrand30(p) (rand30()<(RAND_MAX_30+1.0)*(p))
#define intrand(n) (rand()%(n))
#define intrand30(n) (rand30()%(n))
#define thetarand() (acos(1.0-2.0*rand()/RAND_MAX))
#define radrand1(r) ((r)*sqrt((float)rand()/RAND_MAX))
#define radrand2(r) ((r)*pow((float)rand()/RAND_MAX,0.333333333333))
#define powrand(xmin,power) (xmin*pow(((float)rand()+1.0)/(RAND_MAX+1.0),1.0/(1.0+power)))
#define powrand30(xmin,power) (xmin*pow(((double)rand30()+1.0)/(RAND_MAX_30+1.0),1.0/(1.0+power)))

unsigned int randomize();
float binomrand(int n,float m,float s);
int intrandp(int n,float *p);
int poisrand(float xm);
int poisrandD(double xm);
float binomialrand(float p,int n);
float gaussrand();
void sphererand(float *x,float rad1,float rad2);
void sphererandd(double *x,double rad1,double rad2);
void randtable(float *a,int n,int eq);
void randtableD(double *a,int n,int eq);
void randshuffletable(float *a,int n);
void randshuffletableD(double *a,int n);
void showdist(int n,float low,float high,int bin);

#endif