#include <math.h>
#include <limits.h>
#include <errno.h>

#define        PI 3.141592654
#define        r_ia     16807
#define        r_im     2147483647
#define        r_am     (1.0/r_im)
#define        r_iq     127773
#define        r_ir     2836
#define        r_ntab   32 
#define        r_ndiv   (1+(r_im-1)/r_ntab)
#define        r_eps    1.2e-7
#define        r_rnmx   (1.0-r_eps)

static double gammln(xx) /* ln of gamma function */
double xx;
{
	double x,tmp,ser;
	static double cof[6]={76.18009173,-86.50532033,24.01409822,
		-1.231739516,0.120858003e-2,-0.536382e-5};
	int j;
	x=xx-1.0;
	tmp=x+5.5;
	tmp -= (x+0.5)*log(tmp);
	ser=1.0;
	for (j=0;j<=5;j++) {
		x += 1.0;
		ser += cof[j]/x;
	}
	return -tmp+log(2.50662827465*ser);
}


/* This is a modification of ran1() to be stand alone with its own seed */
/* NOTE: calling with 0, means continue to run, other wiTH the argument is */
/* considered as seed */
double ran1(seed_or_cont)
long seed_or_cont;
{     
       int  j;
       long k;
       static long idum;
       static long iy = 0; 
       static long iv[r_ntab];
       double temp; 

        if (seed_or_cont) {
  	 idum = fabs(seed_or_cont);
	 fprintf(stderr,"RNG: seed is set to %ld\n",idum); 

	 for ( j = r_ntab + 7; j>=0;j--){
	   k = (idum/r_iq);
           idum = r_ia*(idum-k*r_iq)-r_ir*k;
	   if (idum<0) idum += r_im;
	   if (j<r_ntab) iv[j] = idum;
	 }
	 iy = iv[0];
       }
       k = (idum/r_iq);
       idum = r_ia*(idum-k*r_iq)-r_ir*k;
       if (idum<0) idum += r_im;
       j = iy/r_ndiv;
       iy = iv[j];
       iv[j] = idum;
       if ((temp = r_am*iy) >r_rnmx) return r_rnmx;
       else return temp;
}

double gasdev()
{
        static int iset=0;
        static double gset;
        double fac,r,v1,v2;

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

double bnldev(ppr,nnr)
double ppr;
int nnr;
{
	int j,errnosv;
	static int nold=(-1);
	double am,em,g,angle,p,bnl,sq,bt,y;
	static double pold=(-1.0),pc,plog,pclog,en,oldg;
	p=(ppr <= 0.5 ? ppr : 1.0-ppr);
	am=nnr*p;
        if (nnr < 25) {
		bnl=0.0;
		for (j=1;j<=nnr;j++)
			if (ran1(0) < p) bnl += 1.0;
	} else if (am < 1.0) {
		g=exp(-am);
		bt=1.0;
		for (j=0;j<=nnr;j++) {
			bt *= ran1(0);
			if (bt < g) break;
		}
		bnl=(j <= nnr ? j : nnr);
	} else {
		if (nnr != nold) {
			en=nnr;
			oldg=gammln(en+1.0);
			nold=nnr;
		} 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*ran1(0);
					angle=PI*ran1(0);
 				y=tan(angle);

			em=sq*y+am;
			} while (em < 0.0 || em >= (en+1.0));
			em=floor(em);
                        errnosv = errno;
                   		bt=1.2*sq*(1.0+y*y)*exp(oldg-gammln(em+1.0)
				      -gammln(en-em+1.0)+em*plog+(en-em)*pclog); 
                        if (errno != errnosv) { 
/*  fprintf(stderr,"RNG:bnldev: errno was set during call to exp:errno%d\terrnosv%d\n",errno,errnosv);  */
                               errno = errnosv;
			}
		} while (ran1(0) > bt);
		bnl=em;
	}
	if (p != ppr) bnl=nnr-bnl;
	return bnl;
}

#undef PI 
#undef r_ia 
#undef r_im 
#undef r_am 
#undef r_iq  
#undef r_ir 
#undef r_ntab 
#undef r_ndiv 
#undef r_eps 
#undef r_rnmx 


int test_ran1() {
   double j,x;

      FILE *out; 
      out = fopen("test_ran1.out","w");

    ran1(2323723);
    j = 0;
    while (j<10000) {

       fprintf(out,"%lg\n",ran1(0));

      j = j + 1;
    }
      fclose(out);
return 0;
   
}       

int test_gasdev(long seed) {
    double j,x;

      FILE *out; 
      out = fopen("test_gasdev.out","w");

    ran1(seed);
    j = 0;
    while (j<10000) {

       fprintf(out,"%lg\n",gasdev());

      j = j + 1;
    }

      fclose(out);

  return 0; 
}