/* Steven Andrews, 5/12/95, modified 11/98 */
/* Several random number generators */
/* See documentation called random doc */
/* Copyright 2003 by Steven Andrews.  Permission is granted
   for non-commercial use of and modifications to the code. */

#include <stdio.h>
#include "random.h"
#include "math2.h"


unsigned int randomize() {
	unsigned int seed;

	seed=(unsigned int) time(NULL);
	srand(seed);
	return seed; }


float binomrand(int n,float m,float s) {
	float x=0;
	int i;
	
	for(i=1;i<=n;i++)
		x+=(float)rand()/RAND_MAX;
	x=(x-0.5*n)/sqrt(n/12.0);
	return(x*s+m); }


int intrandp(int n,float *p) {
	float r;
	int jl,ju,jm;
	
	r=(float)rand()/(RAND_MAX+1.0);
	jl=-1;
	ju=n-1;
	while(ju-jl>1) {
		jm=(ju+jl)>>1;
		if(p[jm]<=r) jl=jm;
		else ju=jm; }
	return ju; }


int poisrand(float xm) {
	static float sq,alxm,g,oldm=-1.0;
	float em,t,y;

	if(xm<=0) return 0;
	else if(xm<12) {
		if(xm!=oldm) {
			oldm=xm;
			g=exp(-xm); }
		em=0;
		for(t=(float)rand()/RAND_MAX;t>g;t*=(float)rand()/RAND_MAX) {
			em+=1.0; }}
	else {
		if(xm!=oldm) {
			oldm=xm;
			sq=sqrt(2.0*xm);
			alxm=log(xm);
			g=xm*alxm-gammaln(xm+1.0); }
		do {
			do {
				y=tan(PI*rand()/RAND_MAX);
				em=sq*y+xm; } while(em<0);
			em=floor(em);
			t=0.9*(1.0+y*y)*exp(em*alxm-gammaln(em+1.0)-g); } while((float)rand()/RAND_MAX>t); }
	return (int) em; }


int poisrandD(double xm) {
	static double sq,alxm,g,oldm=-1.0;
	float em,t,y;

	if(xm<=0) return 0;
	else if(xm<12) {
		if(xm!=oldm) {
			oldm=xm;
			g=exp(-xm); }
		em=0;
		for(t=(double)rand()/RAND_MAX;t>g;t*=(double)rand()/RAND_MAX) {
			em+=1.0; }}
	else {
		if(xm!=oldm) {
			oldm=xm;
			sq=sqrt(2.0*xm);
			alxm=log(xm);
			g=xm*alxm-gammaln(xm+1.0); }
		do {
			do {
				y=tan(PI*rand()/RAND_MAX);
				em=sq*y+xm; } while(em<0);
			em=floor(em);
			t=0.9*(1.0+y*y)*exp(em*alxm-gammaln(em+1.0)-g); } while((double)rand()/RAND_MAX>t); }
	return (int) em; }


float binomialrand(float p,int n) {
	int swap,j;
	float am,bnl,pscale,g,t,sq,angle,y,em;
	static float nold=-1,pold=-1;
	static float en,oldg,pc,plog,pclog;

	if(n<1) return 0;
	if(p>1) return n;
	if(p<0) return 0;
	swap=0;
	if(p>0.5) {
		swap=1;
		p=1.0-p; }
	am=n*p;
	if(n<25) {
		pscale=p*RAND_MAX;
		bnl=0;
		for(j=0;j<n;j++)
			if(rand()<pscale) bnl+=1.0; }
	else if(am<1.0)	{
		g=exp(-am);
		t=1.0;
		for(j=0;j<=n;j++) {
			t*=(float)rand()/RAND_MAX;
			if(t<g) break; }
		bnl=j<=n?j:n; }
	else {
		if(n!=nold) {
			en=n;
			oldg=gammaln(en+1.0);
			nold=n; }
		if(p!=pold) {
			pc=1.0-p;
			plog=log(p);
			pclog=log(pc);
			pold=p; }
		sq=sqrt(2.0*am*pc);
		do {
			do {
				angle=PI*(float)rand()/RAND_MAX;
				y=tan(angle);
				em=sq*y+am; }
				while(em<0||em>=(en+1.0));
			em=floor(em);
			t=1.2*sq*(1.0+y*y)*exp(oldg-gammaln(em+1.0)-gammaln(en-em+1.0)+em*plog+(en-em)*pclog); }
			while((float)rand()/RAND_MAX>t);
		bnl=em; }
	if(swap) bnl=n-bnl;
	return bnl; }


float gaussrand() {
	static int iset=0;
	static float gset;
	float fac,r,v1,v2;

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


void sphererand(float *x,float rad1,float rad2) {
	float th,phi;

	th=thetarand();
	phi=unirand(0,2.0*PI);
	if(rad1==rad2)
		;
	else if(rad1==0) rad1=radrand2(rad2);
	else rad1=pow((float)rand()/RAND_MAX*(rad2*rad2*rad2-rad1*rad1*rad1)+rad1*rad1*rad1,0.333333333333);
	x[0]=rad1*sin(th)*cos(phi);
	x[1]=rad1*sin(th)*sin(phi);
	x[2]=rad1*cos(th);
	return; }


void sphererandd(double *x,double rad1,double rad2) {
	double th,phi;

	th=thetarand();
	phi=unirand(0,2*PI);
	if(rad1==rad2)
		;
	else if(rad1==0) rad1=radrand2(rad2);
	else rad1=pow((double)rand()/RAND_MAX*(rad2*rad2*rad2-rad1*rad1*rad1)+rad1*rad1*rad1,0.333333333333);
	x[0]=rad1*sin(th)*cos(phi);
	x[1]=rad1*sin(th)*sin(phi);
	x[2]=rad1*cos(th);
	return; }


void randtable(float *a,int n,int eq) {
	int i;

	if(eq==1) {
	  for(i=0;i<n/2;i++)
	    a[i]=sqrt(2)*inversefn(erfn,2.0*(i+0.5)/n-1.0,-20,20,30);
	  for(i=n/2;i<n;i++)
	    a[i]=-a[n-i-1]; }
	return; }


void randtableD(double *a,int n,int eq) {
	int i;

	if(eq==1) {
	  for(i=0;i<n/2;i++)
	    a[i]=sqrt(2)*inversefn(erfn,2.0*(i+0.5)/n-1.0,-20,20,30);
	  for(i=n/2;i<n;i++)
	    a[i]=-a[n-i-1]; }
	return; }


void randshuffletable(float *a,int n) {
	int i,j;
	float x;

	for(i=0;i<n;i++) {
		j=intrand(n);
		x=a[i];
		a[i]=a[j];
		a[j]=x; }
	return; }


void randshuffletableD(double *a,int n) {
	int i,j;
	double x;

	for(i=0;i<n;i++) {
		j=intrand(n);
		x=a[i];
		a[i]=a[j];
		a[j]=x; }
	return; }


void showdist(int n,float low,float high,int bin) {
	int i,a[100],uflow=0,oflow=0,b;
	float x,sum=0,sum2=0;
//	float p[5];
//	float gauss[4096];

//	p[0]=0.1;p[1]=0.2;p[2]=0.7;p[3]=0.7;p[4]=1;
//	randtable(gauss,4096,1);
	if(bin>=100) bin=99;
	for(b=0;b<bin;b++)	a[b]=0;
	for(i=1;i<=n;i++)	{

//		x=binomrand(10,0,1);
//		x=thetarand();
//		x=intrandp(5,p);
//		x=poisrand(20);
		x=10*gaussrand();
//		x=10*gauss[rand()&4095];

		b=floor((x-low)*(bin-1)/(high-low)+0.5);
		if(b<0) uflow++;
		else if(b>=bin) oflow++;
		else a[b]++;
		sum+=x;
		sum2+=x*x; 	}
	printf("<%0.2f\t:",low-(high-low)/(bin-1)/2);
	for(i=1;i<=uflow;i++)	printf("x");
	for(b=0;b<bin;b++)	{
		printf("\n %0.2f\t:",b*(high-low)/(bin-1)+low);
		for(i=1;i<=a[b];i++)	printf("x");	}	
	printf("\n>%0.2f\t:",high+(high-low)/(bin-1)/2);
	for(i=1;i<=oflow;i++)	printf("x");
	printf("\n");
	printf("mean: %f\tstandard deviation: %f\n",sum/n,sqrt(sum2/n-sum/n*sum/n));
	return; }