/* Steven Andrews, 5/12/95.
Several random number generators.
See documentation called random_doc.doc.
Copyright 2003-2007 by Steven Andrews.  This work is distributed under the terms
of the Gnu Lesser General Public License (LGPL). */

#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 Rand_TriPtD(double *pt1,double *pt2,double *pt3,int dim,double *ans) {
	static double x,y;
	static int xd,yd;
	double x0,y0,x1,y1,x2,y2,m01,m02,m12,area,yx1,yy;
	int d,dsmall,zd;
	double range,smrange,swap,coefx,coefy,coefz,coefk;

	if(dim==2) {
		xd=0;
		yd=1;
		if(pt1[0]<=pt2[0]&&pt2[0]<=pt3[0]) {x0=pt1[0];y0=pt1[1];x1=pt2[0];y1=pt2[1];x2=pt3[0];y2=pt3[1];}
		else if(pt1[0]<=pt3[0]&&pt3[0]<=pt2[0]) {x0=pt1[0];y0=pt1[1];x1=pt3[0];y1=pt3[1];x2=pt2[0];y2=pt2[1];}
		else if(pt2[0]<=pt3[0]&&pt3[0]<=pt1[0]) {x0=pt2[0];y0=pt2[1];x1=pt3[0];y1=pt3[1];x2=pt1[0];y2=pt1[1];}
		else if(pt2[0]<=pt1[0]&&pt1[0]<=pt3[0]) {x0=pt2[0];y0=pt2[1];x1=pt1[0];y1=pt1[1];x2=pt3[0];y2=pt3[1];}
		else if(pt3[0]<=pt1[0]&&pt1[0]<=pt2[0]) {x0=pt3[0];y0=pt3[1];x1=pt1[0];y1=pt1[1];x2=pt2[0];y2=pt2[1];}
		else if(pt3[0]<=pt2[0]&&pt2[0]<=pt1[0]) {x0=pt3[0];y0=pt3[1];x1=pt2[0];y1=pt2[1];x2=pt1[0];y2=pt1[1];}
		else {x0=pt1[0];y0=pt1[1];x1=pt2[0];y1=pt2[1];x2=pt3[0];y2=pt3[1];}
		if(x0==x2) {
			ans[0]=x0;
			x0=(y0<y1)?(y0<y2?y0:y2):(y1<y2?y1:y2);
			x2=(y0>y1)?(y0>y2?y0:y2):(y1>y2?y1:y2);
			ans[1]=unirand30(x0,x2);
			return; }
		m01=(y1-y0)/(x1-x0);
		m02=(y2-y0)/(x2-x0);
		m12=(y2-y1)/(x2-x1);
		area=0.5*(y1-y0-m02*(x1-x0))*(x2-x0);
		if(x1-x0>x2-x1) yx1=0.5/area*(m01-m02)*(x1-x0)*(x1-x0);
		else yx1=1.0-0.5/area*(m02-m12)*(x1-x2)*(x1-x2);
		yy=unirand30(0,1);
		if(yy<yx1) {
			x=x0+sqrt(2.0*area*yy/(m01-m02));
			y=unirand30(y0+m02*(x-x0),y0+m01*(x-x0)); }
		else if(yy>yx1) {
			x=x2-sqrt(2.0*area*(1.0-yy)/(m02-m12));
			y=unirand30(y0+m02*(x-x0),y1+m12*(x-x1)); }
		else {
			x=x1;
			y=unirand30(y0+m02*(x-x0),y1); }
		ans[0]=x;
		ans[1]=y;
		return; }

	smrange=-1;
	dsmall=0;
	for(d=0;d<dim;d++) {									// find coordinate with smallest used range
		range=(pt1[d]>pt2[d])?(pt1[d]>pt3[d]?pt1[d]:pt3[d]):(pt2[d]>pt3[d]?pt2[d]:pt3[d]);
		range-=(pt1[d]<pt2[d])?(pt1[d]<pt3[d]?pt1[d]:pt3[d]):(pt2[d]<pt3[d]?pt2[d]:pt3[d]);
		if(range<smrange||smrange==-1) {
			smrange=range;
			dsmall=d; }}
	swap=pt1[dim-1];pt1[dim-1]=pt1[dsmall];pt1[dsmall]=swap;	// swap dsmall and last
	swap=pt2[dim-1];pt2[dim-1]=pt2[dsmall];pt2[dsmall]=swap;
	swap=pt3[dim-1];pt3[dim-1]=pt3[dsmall];pt3[dsmall]=swap;
	Rand_TriPtD(pt1,pt2,pt3,dim-1,ans);												// find reduced dimensional answer

	zd=dim-1;
	coefx=pt1[yd]*(pt2[zd]-pt3[zd])+pt2[yd]*(pt3[zd]-pt1[zd])+pt3[yd]*(pt1[zd]-pt2[zd]);
	coefy=pt1[zd]*(pt2[xd]-pt3[xd])+pt2[zd]*(pt3[xd]-pt1[xd])+pt3[zd]*(pt1[xd]-pt2[xd]);
	coefz=pt1[xd]*(pt2[yd]-pt3[yd])+pt2[xd]*(pt3[yd]-pt1[yd])+pt3[xd]*(pt1[yd]-pt2[yd]);
	coefk=-(pt1[xd]*(pt2[yd]*pt3[zd]-pt3[yd]*pt2[zd])+pt2[xd]*(pt3[yd]*pt1[zd]-pt1[yd]*pt3[zd])+pt3[xd]*(pt1[yd]*pt2[zd]-pt2[yd]*pt1[zd]));
	if(coefz!=0) ans[zd]=-(coefk+coefx*x+coefy*y)/coefz;
	else ans[zd]=pt1[zd];

	swap=pt1[dim-1];pt1[dim-1]=pt1[dsmall];pt1[dsmall]=swap;	// swap dsmall and last
	swap=pt2[dim-1];pt2[dim-1]=pt2[dsmall];pt2[dsmall]=swap;
	swap=pt3[dim-1];pt3[dim-1]=pt3[dsmall];pt3[dsmall]=swap;
	swap=ans[dim-1];ans[dim-1]=ans[dsmall];ans[dsmall]=swap;

	if(xd==dsmall) xd=dim-1;
	if(yd==dsmall) yd=dim-1;

	return; }



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

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


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

	if(eq==1) {
		dy=2.0/n;
	  for(i=0;i<n/2;i++)
	    a[i]=SQRT2*inversefnD(erfnD,(i+0.5)*dy-1.0,-20,20,30);
	  for(i=n/2;i<n;i++)
	    a[i]=-a[n-i-1]; }
	else if(eq==2) {
		dy=1.0/SQRTPI/n;
		for(i=0;i<n;i++)
			a[i]=SQRT2*inversefnD(erfcintegralD,(i+0.5)*dy,0,20,30); }
	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; }