TITLE Random Number Generator based on numerical rec. implementing ran1(),gasdev(), and bnldev()

COMMENT
Note several modifcations from the original Numerical Recipes source:
1. idum is a static variable of ran1() and is not sent as an argument.
   if called with argument 0, ran1(0) this means repick. If called with non-zero argument,
   the argument is treated as a new seed.
2. gasdev() and bnldev() does not pass idum to ran1().
3. Since we cannot control the order of which the mod files are 
   initialized the initialization of ran1() is checked thorugh a flag variable of this RNG.mod file,
   initran1, and if it was called before INITIAL block is called it runs the same code. 
4. Since bnldev have underflow errors when computing exp(oldg-gammln(em+1.0)
			      -gammln(en-em+1.0)+em*plog+(en-em)*pclog)  
    this set errno, and causes neuron to print an error message. Thus errno is checked
    before and after this command and if changed it being set back to its old value.
ENDCOMMENT

INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}

NEURON {
	POINT_PROCESS RNG
	NONSPECIFIC_CURRENT i
	RANGE seed
	GLOBAL dummy
}

ASSIGNED { 
     dummy
     i 
}

PARAMETER {
	seed  = 0
	initran1 = 0
}

INITIAL {
  VERBATIM
  if ((long) seed == 0) {
    fprintf(stderr,"Please provide a seed to random number generator!\n"); 
  }
  if ((long) seed > 0) {
	seed = -1*seed;
  }  
  ran1((long)seed);

  ENDVERBATIM    
}


BREAKPOINT {
 i = 0 
}

VERBATIM

#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 wise 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 
ENDVERBATIM

PROCEDURE test_ran1() {
    LOCAL j,x
      VERBATIM 
      FILE *out; 
      out = fopen("test_ran1.out","w");
      ENDVERBATIM 

    j = 0
    while (j<10000) {
      VERBATIM 
       fprintf(out,"%lg\n",ran1(0));
      ENDVERBATIM 
      j = j + 1
    }
      VERBATIM 
      fclose(out);
      ENDVERBATIM 
   
}       

PROCEDURE test_gasdev() {
    LOCAL j,x
      VERBATIM 
      FILE *out; 
      out = fopen("test_gasdev.out","w");
      ENDVERBATIM 

    j = 0
    while (j<10000) {
      VERBATIM 
       fprintf(out,"%lg\n",gasdev());
      ENDVERBATIM 
      j = j + 1
    }
      VERBATIM 
      fclose(out);
      ENDVERBATIM 
   
}