/***********************************************************
 Basic statistical, input-output and matrix manipulation

 Author. Peter Mueller, Stephen Morris, David Rossell
 Last modified. 04 2006
***********************************************************/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "cstat.h"

static const char interface_c_sccs_id[] = "%W%";
static const char mess_c_sccs_id[] = "%W%";
static const char nrutil_c_sccs_id[] = "%W%"; 
static const char vector_c_sccs_id[] = "%W%";
static const char css_c_sccs_id[] = "@(#)$Workfile: rand.c$ $Revision: 5$";
static long is1=123456789, is2=981963; 
static int set=0;

FILE	*ifile, *fopen(), *ofile;
int nv = 0;

long Xm1,Xm2,Xa1,Xa2,Xcg1[32],Xcg2[32],Xa1w,Xa2w,Xig1[32],Xig2[32],Xlg1[32],
    Xlg2[32],Xa1vw,Xa2vw;
long Xqanti[32];


/************************************************************************
                         STATISTICAL FUNCTIONS
************************************************************************/

/* Sample mean of elements 0 through lim of vector x */
double meani(int *x, int lim) 
{
    int i;
    double value;
    for(i=0,value=0; i<=lim; i++) { value += x[i]; }
    value *= 1.0/(lim+1.0);
    return value;
}

double wmeani(int *x, int lim, double *w)
{
    int i;
    double value, wtot;
    for(i=0,value=0,wtot=0; i<=lim; i++) { value += w[i]*x[i]; wtot += w[i]; }
    value *= 1.0/wtot;
    return value;
}

/* Sample mean of elements 0 through lim of vector x */
double meanx(double *x, int lim) 
{
    int i;
    double value;
    for(i=0,value=0; i<=lim; i++) value += x[i];
    value *= 1.0/(lim+1.0);
    return value;
}

double wmeanx(double *x, int lim, double *w)
{
    int i;
    double value, wtot;
    for(i=0,value=0,wtot=0; i<=lim; i++) {
        value += w[i]*x[i]; wtot += w[i];
    }
    value *= 1.0/wtot;
    return value;
}

/* Sample variance of elements 0 through lim of vector x, unb=1 returns unbiased version */
double vari(int *x, int lim, int unb) 
{
    int i;
    double value;
    for(i=0,value=0; i<=lim; i++) { value += pow(x[i],2) / (1.0+lim); }
    value -= pow(meani(x,lim),2);
    if (unb==1 && lim>0) value *= (1.0+lim)/(0.0+lim);
    return value;
}

double wvari(int *x, int lim, double *w) 
{
    int i;
    double value, wtot;
    for(i=0,value=0,wtot=0; i<=lim; i++) { value += w[i]*pow(x[i],2); wtot += w[i]; }
    value = value/wtot - pow(wmeani(x,lim,w),2);
    return value;
}


/* Sample variance of elements 0 through lim of vector x, unb=1 returns unbiased version */
double varx(double *x, int lim, int unb) 
{
    int i;
    double value;
    for(i=0,value=0; i<=lim; i++) { value += pow(x[i],2) / (1.0+lim); }
    value -= pow(meanx(x,lim),2);
    if (unb==1 && lim>0) value *= (1.0+lim)/(lim+0.0);
    return value;
}

double wvarx(double *x, int lim, double *w) 
{
    int i;
    double value, wtot;
    for(i=0,value=0,wtot=0; i<=lim; i++) { value += w[i]*pow(x[i],2); wtot += w[i]; }
    value = value/wtot - pow(wmeanx(x,lim,w),2);
    return value;
}


/************************************************************************
                         INPUT/OUTPUT FUNCTIONS
************************************************************************/

/* open file for input */
FILE *openIn(char *name)
{
  if ((ifile=fopen(name,"r"))==NULL){
    fserror("openIn","open file",name);
  }
  return(ifile);
}

/* open file for output */

FILE *openOut(char *name)
{
  if ((ofile=fopen(name,"w"))==NULL){
    fserror("openOut","open file",name);
  }
  return ofile;
}


/* --------------------   read in   ------------------------- */

void scanFloat(char *txt, float *f)
{
  fscanf(ifile,txt);
  if (fscanf(ifile," %f ",f) != 1) {
    fserror ("scanFloat","read float",txt);
  }
}

void scanDouble(char *txt,double *f)
{
  fscanf(ifile,txt);
  if (fscanf(ifile," %lf ",f) != 1) {
    fserror ("scanDouble","read double",txt);
  }
}

void fscanDouble(FILE *ifile, char *txt, double *f)
{
  fscanf(ifile,txt);
  if (fscanf(ifile," %lf ",f) != 1) {
    fserror ("fscanDouble","read double",txt);
  }
}

void scanInt(char *txt, int *n)
{
  int 	s;

  fscanf(ifile,txt);
  if ((s = fscanf(ifile," %d ",n)) != 1) {
    fserror ("scanInt","read int",txt);
  }
}

void fscanInt(FILE *ifile, char *txt, int *n)
{
	int 	s;

  fscanf(ifile,txt);
  if ((s = fscanf(ifile," %d ",n)) != 1) {
    fserror ("fscanInt","read int",txt);
  }
}


void scanLong(char *txt, long *n)
{
        int     s;

  fscanf(ifile,txt);
  if ((s = fscanf(ifile," %ld ",n)) != 1) {
    fserror ("scanLong","read long",txt);
  }
}

/* --------------------   read arrays   ------------------------- */

void scanFloatArray(char *txt, float *x, int n)
{
	scanArray(txt,x,n);
}

void scanArray(char *txt, float *x, int n)
{
	int	i; 

  fscanf(ifile,txt);
  for(i=0;i<n;i++){
    if (fscanf(ifile," %f ",&x[i]) != 1) {
      fserror ("scanArray","read float array",txt);
    }
  }
}

void scanDoubleArray(char *txt, double *x, int n)
{
  int	i;

  fscanf(ifile,txt);
  for(i=0;i<n;i++){
    if (fscanf(ifile," %lg ",&x[i]) != 1) {
      fserror ("scanDoubleArray",
	       "read double array",txt);
    }
  }
}

void fscanDoubleArray(FILE *in, double *x, int n)
{
  int	i;

  for(i=0;i<n;i++){
    if (fscanf(in," %lg ",&x[i]) != 1) {
      /* printf("i=%d\n",i); */
      fserror("fscanDoubleArray","read double array","");
    }
  }
}

void scanString(char *txt, char *s, int n)
{
  fgets(s,n,ifile);
}
  

void fscanString(FILE *ifile, char *txt, char *s, int n)
{
  fgets(s,n,ifile);
}
  
void scanDoubleMatrix(char *txt,double **x,int r,int c)
{
  int	i,j;
  
  fscanf(ifile,txt);
  for(i=0;i<r;i++)
    for(j=0;j<c;j++){
      if (fscanf(ifile," %lg ",&x[i][j]) != 1) {
	fserror ("scanDoubleMatrix","read double matrix",txt);
      }
    }
}

void fscanDoubleMatrix(FILE *ifile, double **x,int r,int c)
{
  int	i,j;
  
  for(i=0;i<r;i++)
    for(j=0;j<c;j++){
      if (fscanf(ifile," %lg ",&x[i][j]) != 1) {
	exit(1);
      }
    }
}

void scanIntArray(char *txt, int *x, int n)
{
  int	i;

  fscanf(ifile,txt);
  for(i=0;i<n;i++){
    if (fscanf(ifile," %d ",&x[i]) != 1) {
      fserror ("scanIntArray","read int array",txt);
    }
  }
}

void fscanIntArray(FILE *ifile, int *x, int n)
{
  int	i;

  for(i=0;i<n;i++){
    if (fscanf(ifile," %d ",&x[i]) != 1) {
      fserror("fscanIntArray","read int array","");
    }
  }
}


/* ------------------------  write scalars  ------------------------ */
void writeInt(int i)
{
  int s;
  s=fprintf(ofile,"%d\n",i);
  if(s<0)
    fserror("writeInt","write int","");
}

void writeLong(long i)
{
  int s;
  s=fprintf(ofile,"%ld\n",i);
  if (s<0)
    fserror("writeLong","write long","");
    
}


void writeFloat(float x)
{
  int s;
  s=fprintf(ofile,"%f\n",x);
  if (s<0)
    fserror("writeFloat","write float","");
  
}

void writeDouble(double x)
{
  int s;
  s=fprintf(ofile,"%5.3e\n",x);
  if (s<0)
    fserror("writeDouble","write double","");
}


/* -----------------------  write arrays   --------------------- */

void fwriteDoubleArray(FILE *f, double *x, int rows, int cols)
{
  int	i,j,s1,s2;
  
  s1 = 0;
  for(i=0;i<rows;i++){
    for(j=0;j<cols;j++){
      if(j%10 == 9)
	fprintf(f,"\n\t");
      s1=fprintf(f,"%5.3e ",x[i*cols+j]);
      if (s1<0) break;
    }
    s2=fprintf(f,"\n");
    if ((s2<0)|(s1<0))
      fserror("fwriteDoubleArray","write double array","");
  }
}

void fwriteIntArray(FILE *f, int *x, int rows, int cols)
{
  int	i,j,s1,s2;
  
  s1 = 0;
  for(i=0;i<rows;i++){
    for(j=0;j<cols;j++){
      if(j%10 == 9)
	fprintf(f,"\n\t");
      s1=fprintf(f,"%d\t",x[i*cols+j]);
      if (s1<0) break;
    }
    s2=fprintf(f,"\n");
    if ((s2<0)|(s1<0))
      fserror("fwriteIntArray","write int array","");
  }
}


void writeIntArray(int *x, int rows, int cols)
{
  fwriteIntArray(ofile,x,rows,cols);
}

void fwriteIntMatrix(FILE *f, int **x, int rows, int cols)
{
  int	i,j,s1;
  
  for(i=0;i<rows;i++){
    for(j=0;j<cols;j++){
      if(j%10 == 9)
	fprintf(f,"\n\t");
      s1=fprintf(f,"%d\t",x[i][j]);
      if (s1<0) 
	fserror("fwriteIntMatrix","write int matrix","");
    }
    fprintf(f,"\n");
  }
}

void writeIntMatrix(int **x, int rows, int cols)
{
  fwriteIntMatrix(ofile,x,rows,cols);
}

void writeDoubleArray(double *x,int rows,int cols)
{
  fwriteDoubleArray(ofile,x,rows,cols);
}

void fwriteDoubleMatrix2(FILE *f, double **x, int rows, int cols)
{
  int	i, j, s;
  
  for(i=0;i<rows;i++){
    for(j=0;j<cols;j++){
      if(j%10 == 9){
	fprintf(f,"\n\t");
      }
      s=fprintf(f,"%5.3e ",x[i][j]);
      if (s<0)
	fserror("fwriteDoubleMatrix2","write double matrix","");
    }
    fprintf(f,"\n");
  }
}

void writeDoubleMatrix2(double **x, int rows, int cols)
{
  fwriteDoubleMatrix2(ofile,x,rows,cols);
}

void writeDoubleMatrix(double **x, int rows, int cols)
{
  int	i,j, c,s;
  
  for(i=0;i<rows;i++){
    for(j=0, c=0;j<cols;j++){
      if(++c > 10){
	fprintf(ofile,"\n\t");
	c = 0;
      }
      s=fprintf(ofile,"%5.3e ",x[i][j]);
      if (s<0)
	fserror("fwriteDoubleMatrix","write double matrix","");
    }
    fprintf(ofile,"\n");
  }
}

void writeFloatArray(float *x,int rows,int cols)
{
 writeArray(x,rows,cols);
}

void writeArray(float *x,int rows,int cols)
{
  int	
    i,j, c,s;
  
  for(i=0;i<rows;i++){
    for(j=0,c=0;j<cols;j++){
      if (c++>9){
	fprintf(ofile,"\n\t");
	c = 0;
      }
      s=fprintf(ofile,"%5.3e ",x[i*cols+j]);
      if (s<0)
	fserror("fwriteDoubleMatrix","write double matrix","");
    }
    fprintf(ofile,"\n");
  }
}

/***************************   error   *************************** */ 

void fserror(char *proc, char *act, char *what){
  /* writes error message and aborts */
  
  fprintf(stderr, "\n ** Error ");
  if (proc[0]!='\0') /* not empty */
    fprintf(stderr, " in function '%s', ", proc);
  if (act[0]!='\0') /* not empty */
    fprintf(stderr, " trying to %s", act);
  if (what[0]!='\0') /* not empty */
    fprintf(stderr, " '%s'", what);
  fprintf(stderr, "\n ** .. exiting program");
  fprintf(stderr, " (from a function in 'interface.c').\n");
  
  exit(1);
}

/******************************************************************************************
                                      DEBUG MESSAGES ETC.
******************************************************************************************/

void error(char *module, char *mess, int nr)              
{
   	printf("\n *** ERROR # %d in %s***\n %s",nr,module, mess);
	printf(  " exiting program \n");
	exit(1);
}

void err_msg(char *fct, char *txt, int n1, int n2, int n3)
{ /* print error message */
  printf("\n\n *** Error in %s \n", fct);
  printf(txt,n1,n2,n3); /* n1,n2 allows to include numbers in txt */
  printf("\n");
  exit(1);
}



/******************************************************************************************
                          MEMORY ALLOCATION, FROM NUMERICAL RECIPIES
******************************************************************************************/

double  *dvector(int nl,int nh) 
{ 
        double  *v; 
 
	nv += (nh-nl+1); 
        v=(double  *)calloc((unsigned) (nh-nl+1),sizeof(double)); 
        if (!v)  
	  nrerror("dvector","allocate a double vector",""); 
        return v-nl; 
} 

double  **dmatrix(int nrl,int nrh,int ncl,int nch) 
{ 
        int i; 
        double  **m; 
	 
	nv += (nrh-nrl+1)*(nch-ncl+1); 
        m=(double  **)calloc((unsigned) (nrh-nrl+1), 
				   sizeof(double  *)); 
        if (!m)  
	  nrerror("dmatrix","allocate a double matrix (1st dim)",""); 
        m -= nrl; 
 
        for(i=nrl;i<=nrh;i++) { 
                m[i]=(double  *)calloc((unsigned) (nch-ncl+1), 
					     sizeof(double)); 
                if (!m[i])  
	  	  nrerror("dmatrix","allocate a double matrix (2nd dim)",""); 
                m[i] -= ncl; 
        } 
        return m; 
} 

double ***darray_3(int lo, int hi) 
{ 
  double ***m; 
 
  nv += (hi-lo+1); 
  m=(double ***)calloc((unsigned) (hi-lo+1),sizeof(double **)); 
  if (!m)  
    nrerror("darray_3","allocate a 3dim double array ",""); 
  m -= lo; 
  return m; 
} 

double ***darray3(int n, int p, int q) 
{ 
  double ***a; 
  int i; 
  a = darray_3(0,n); 
  for(i=0;i<n;i++) 
    a[i] = dmatrix(0,p,0,q); 
  return(a); 
} 

int  *ivector(int nl,int nh) 
{ 
        int  *v; 
 
	nv += (nh-nl+1); 
        v=(int  *)calloc((unsigned) (nh-nl+1),sizeof(int)); 
        if (!v) nrerror("ivector","allocate an int vector",""); 
        return v-nl; 
} 

int  **imatrix(int nrl,int nrh,int ncl,int nch) 
{ 
        int i, **m; 
 
	nv += (nrh-nrl+1)*(nch-ncl+1); 
        m=(int  **)calloc((unsigned) (nrh-nrl+1),sizeof(int  *)); 
        if (!m)  
	  nrerror("imatrix","allocate a int matrix (1st dim).",""); 
        m -= nrl; 
 
        for(i=nrl;i<=nrh;i++) { 
                m[i]=(int  *)calloc((unsigned) (nch-ncl+1),sizeof(int)); 
                if (!m[i])  
	  	  nrerror("imatrix","allocate a int matrix (2nd dim).",""); 
                m[i] -= ncl; 
        } 
        return m; 
} 

int ***iarray_3(int lo, int hi) 
{ 
  int ***m; 
 
  nv += (hi-lo+1); 
  m=(int ***)calloc((unsigned) (hi-lo+1),sizeof(int **)); 
  if (!m)  
    nrerror("iarray_3","allocate a 3dim int array ",""); 
  m -= lo; 
  return m; 
}

int ***iarray3(int p1, int p2, int p3) 
{ 
  int ***m, i; 
 
  m = iarray_3(0,p1); 
  for (i=0;i<p1;i++) 
    m[i] = imatrix(0,p2,0,p3); 
  return m; 
} 

void free_dvector(double  *v,int nl,int nh) 
{ 
        if( (v+nl) != NULL ) free((char  *) (v+nl)); 
	nv -= (nh-nl+1); 
} 

void free_ivector(int  *v,int nl,int nh) 
{ 
        if( (v+nl) != NULL ) free((char  *) (v+nl)); 
	nv -= (nh-nl+1); 
} 

void free_dmatrix(double  **m,int nrl,int nrh,int ncl,int nch) 
{ 
        int i; 
 
        for(i=nrh;i>=nrl;i--) {if( (m[i]+ncl) != NULL )  
				 free((char  *) (m[i]+ncl));} 
        if( (m+nrl) != NULL ) free((char  *) (m+nrl)); 
        nv -= (nch-ncl+1)*(nrh-nrl+1); 
} 
 
void free_imatrix(int  **m,int nrl,int nrh,int ncl,int nch) 
{ 
        int i; 
 
        for(i=nrh;i>=nrl;i--) {if( (m[i]+ncl) != NULL )  
				 free((char  *) (m[i]+ncl));} 
        if( (m+nrl) != NULL ) free((char  *) (m+nrl)); 
        nv -= (nch-ncl+1)*(nrh-nrl+1); 
} 


void nrerror(char *proc, char *act, char *what) 
{ 
  void exit(); 
 
  fprintf(stderr, "\n ** Error "); 
  if (proc[0]!='\0') /* not empty */ 
    fprintf(stderr, " in function '%s', ", proc); 
  if (act[0]!='\0') /* not empty */ 
    fprintf(stderr, " trying to %s", act); 
  if (what[0]!='\0') /* not empty */ 
    fprintf(stderr, " '%s',", what); 
  else 
    fprintf(stderr, ", "); 
  fprintf(stderr, "\n ** .. exiting program.\n"); 
  fprintf(stderr, "\n ** (a function in 'nrutil.c').\n"); 
  exit(1); 
} 


/************************************************************************
                          MATHEMATICAL FUNCTIONS
************************************************************************/


// digamma(x) returns the derivative of the natural log of the gamma function, e.g. gamma'(x)/gamma(x)
//  x must be positive and finite

double digamma(double x) {

  double stirling[] = {
    -8.333333333333333e-02, 8.333333333333333e-03,
    -3.968253968253968e-03, 4.166666666666667e-03,
    -7.575757575757576e-03, 2.109279609279609e-02,
    -8.333333333333334e-02, 4.432598039215686e-01,
    -3.053954330270120e+00, 2.645621212121212e+01,
    -2.814601449275362e+02, 3.607510546398047e+03 
  };

  long i;
  double lower= 1.0e-8, upper= 19.5, euler_one= .422784335098467139393488, ans, x_inv, x_pow;


  if (x<=0) error("digamma",  "digamma must be given a positive argument", 1);;

  if (x<lower) {
    ans = -1.0 / x - 1.0/(1.0+x) + euler_one;
    return(ans);
  }

  ans= 0.0;
  while (x<upper) {
    ans= ans - 1.0/x;
    x= x + 1.0;
  }

  x_inv= 1.0/x;
  ans= ans + log(x) - 0.5*x_inv;

  x_inv= x_inv*x_inv;
  x_pow= x_inv;

  for (i=0; i<12; i++) {
    ans= ans + stirling[i] * x_pow;
    x_pow= x_pow * x_inv;
  }
  return(ans);

}

/************************************************************************
                            VECTOR ALGEBRA
************************************************************************/

void grid (double x0, double xn, int n, double *x) 
{ 
  int i; 
  double dx,xi; 
  dx = (xn-x0)/(n-1.0); 
  for(i=0,xi=x0;i<n;i++,xi+=dx) 
    x[i]=xi; 
} 

void  rA(double r,double **A, double **B, int p, int q)
{
  int _i, _j;
  for(_i=0; _i<(p); _i++){			 
    for(_j=0; _j<(q); _j++)			 
      B[_i][_j] = r* A[_i][_j];	 
  } 
} 

void R_zero(double **A, int p,int q) 
{
  int _i, _j;
  for(_i=0;_i<p;_i++) 
    for(_j=0;_j<q;_j++) 
      A[_i][_j] = 0.0; 
} 

void R_izero(int **A, int p,int q) 
{
  int _i, _j;
  for(_i=0;_i<p;_i++) 
    for(_j=0;_j<q;_j++) 
      A[_i][_j] = 0; 
} 

double max_xy(double x, double y) 
{ 
  return (x>y) ? x : y; 
} 

int dcompare (const void *a, const void *b) {
const double *da = (const double *) a;
const double *db = (const double *) b;
     
return (*da > *db) - (*da < *db);
}

void dvecsort(double *v, int size) {
qsort (v, size, sizeof (double), dcompare);
}


void dindexsort(double *x, int *index, int ilo, int ihi, int incr) {
/* Sorts vector of doubles x by rearranging values in index with quicksort algorithm e.g. x[index[ilo]], x[index[ilo+1]]... x[index[ihi]] is ordered */
/* Input
   - x: vector of doubles that we want to order from position ilo to ihi
   - index: vector of integers indexing the values of x
   - ilo: first element of x we want to order
   - ihi: last element of x we want to order
   - incr: for incr==1 x is returned in increasing order; incr==-1 in decreasing order
   Output: vector index rearranged so that x[index[lo]], x[index[lo+1]]... x[index[ihi]] is ordered
*/

int pivot;              // pivot value for partitioning array
int ulo, uhi;           // indices at ends of unpartitioned region
int tempEntry;          // temporary entry used for swapping
int sortlo, sortup;      // indicate if sub-vectors are sorted so no further subdivision is needed

if (ilo >= ihi) { return; }

sortlo= sortup= 1;
pivot = (ilo + ihi)/2;                                                // Select a pivot value
ulo = ilo; uhi = ihi;                                                 // Initialize ends of unpartitioned region
// While the unpartitioned region is not empty, try to reduce its size.
while (ulo < uhi) {
  if ((x[index[uhi]]*incr) > (x[index[pivot]]*incr)) {                // Here, we can reduce the size of the unpartitioned region and try again.
    if ((uhi<ihi) && ((x[index[uhi]]*incr)>(x[index[uhi+1]]*incr))) sortup= 0;      // Check if upper subvector is ordered
    uhi--;
    if ((uhi==pivot) && (ulo<pivot))  { tempEntry= index[pivot]; index[pivot]= index[pivot-1]; index[pivot-1]= tempEntry; pivot--; } 
  } else {                                                            // Here, x[index[uhi]] <= x[index[pivot]], so swap entries at indices ulo and uhi.
    tempEntry = index[ulo]; index[ulo] = index[uhi]; index[uhi] = tempEntry;
    if (pivot==ulo) pivot= uhi;
    if ((ulo>ilo) && ((x[index[ulo]]*incr)<(x[index[ulo-1]]*incr))) sortlo= 0;
    ulo++;                                                            // Reduce the size of the unpartitioned region
    if ((ulo==pivot) && (uhi>(pivot+1)))  { tempEntry= index[pivot]; index[pivot]= index[pivot+1]; index[pivot+1]= tempEntry; pivot++; }
  }
}

// Entries from ilo to pivot - 1 are < or > pivot and from pivot+1 to ihi are > or < pivot. The two regions can be sorted recursively.
if ((sortlo==0) && (ilo<(pivot-1))) dindexsort(x, index, ilo, pivot - 1, incr);
if ((sortup==0) && (ihi>(pivot+1))) dindexsort(x, index, pivot + 1, ihi, incr);

}


/************************************************************************
                       RANDOM VARIATE GENERATION
************************************************************************/

/* call setall(is1,is2) */ 
void setseed(long is1, long is2) 
{ 
  set=1;   
  setall(is1,is2); 
} 

double runif() 
{ 
  double  
    x; 
 
  if (set==0){ 
    setall(is1,is2); 
    set=1; 
  } 
   
  /* assign to double x for conversion */ 
  x = genunf(0.0,1.0); 
  return(x); 
} 

double betadev(double alpha, double beta) 
{ 
      double gamdev(); 
      double x, y; 
       
      x = gamdev(alpha);               /* X ~ gamma(alpha) */ 
      y = gamdev(beta);                /* Y ~ gamma(beta)  */ 
 
      return (x/(x+y));     /* X/(X+y) ~ ebta(alpha, beta) */ 
} 

double gamdev(double alpha) 
{ 
  double  
    value; 
  float 
    a; 
  float  
    gengam(float,float); 
   
  a = alpha; /* type conversion */ 
  value = gengam(1.0,a); 
  return value; 
} 

/* *************************************************' 
   normal cdf and inv cdf 
 ************************************************* */ 
double	cdfnorm(double y, double m, double s) 
/* returns cdf of normal N(m,s^2) at x */ 
{	 
  double  /* used to all be float... */
    p, q, mean,sd,bound,x,z; 
  double  
    cdf; 
  int  
    status,which; 
 
  /* primitive type conversion */ 
  x = y;  
  mean = m; 
  sd = s; 
  which = 1; 
  z = (x-mean)/sd; 
 
  if (z < -5.0) 
    p = 2.86e-7F; 
  else if (z > 5.0) 
    p = 0.9999997F; 
  else 
    cdfnor(&which,&p, &q,     &x,     &mean,  &sd,    &status,    &bound); 
   
  cdf = p; /* another primitive type conversion */ 
  return cdf; 
}

double	invcdfn (double cdf, double m, double s) 
/* returns inv cdf of normal N(m,s^2) at p */ 
{	 
  double /* used to all be floats.. */
    p,q,mean,sd,bound,x; 
  double  
    y; 
  int  
    status,which; 
 
  if( (cdf < 0.0) | (cdf > 1.0) ){ 
    error("invcdfn",  "Tried inverse cdf with p<0 or p>1", 1); 
  }  
  /* par check */ 
  if (cdf <= 2.86e-07) 
    y = -5.0*s+m; 
  else if (cdf >= 0.9999997) 
    y = 5.0*s+m; 
  else{ 
    /* primitive type conversion */ 
    p = cdf; 
    q = 1.0-p;
    mean = m; 
    sd = s; 
    which = 2; 
     
    cdfnor(&which,&p,&q, &x,&mean,&sd,&status,&bound); 
   
    y = x; /* another primitive type conversion */ 
  } 
  return y; 
} 


/* returns draw from binomial(n,p) */  
int	binomial(int n, double p) 
{ 	 
  int i,x; 
  for(i=0,x=0;i<n;i++) 
    x += (runif() < p) ? 1 : 0; 
  return x; 
} 


/* returns draw from multinomial with cell prob pr  
  ------------------------------------------   
  Value: 
    x:   vector of indices indicating draws 
         x in [0,..n_cell-1] 
  ------------------------------------------   
  Args: 
    n_draw: number of draws 
    n_cell: number of cells 
    pr:     n_cell vector of cell probs 
            (not necessarily standardized) 
    x:      n_draw vector of indices 
            (OUTPUT) 
*/ 
void multinomial(int n_draw, int n_cell, double *pr, int *x) 
{ 
  double *cum_p, uj; 
  int i,j; 
 
  cum_p = dvector(0,n_cell); 
 
  for(i=1,cum_p[0]=pr[0];i<n_cell;i++)     
    cum_p[i] = cum_p[i-1]+pr[i]; 
  for(j=0;j<n_draw;j++){ 
    uj = runif()*cum_p[n_cell-1]; 
    for(i=0; ((uj > cum_p[i]) & (i<n_cell)); i++); 
    x[j] = i; 
  } 
 
  free_dvector(cum_p,0,n_cell); 
} 



/************************************************************************
                       MORE RANDOM VARIATE STUFF
************************************************************************/

/************************************************************************
FIFDINT:
Truncates a double precision number to an integer and returns the
value in a double.
************************************************************************/
double fifdint(double a)
/* a     -     number to be truncated */
{
  long temp;
  temp = (long)(a);
  return (double)(temp);
}


void cdfnor(int *which,double *p,double *q,double *x,double *mean,
	    double *sd,int *status,double *bound)
/**********************************************************************

      void cdfnor(int *which,double *p,double *q,double *x,double *mean,
            double *sd,int *status,double *bound)

               Cumulative Distribution Function
               NORmal distribution


                              Function


     Calculates any one parameter of the normal
     distribution given values for the others.


                              Arguments


     WHICH  --> Integer indicating  which of the  next  parameter
     values is to be calculated using values  of the others.
     Legal range: 1..4
               iwhich = 1 : Calculate P and Q from X,MEAN and SD
               iwhich = 2 : Calculate X from P,Q,MEAN and SD
               iwhich = 3 : Calculate MEAN from P,Q,X and SD
               iwhich = 4 : Calculate SD from P,Q,X and MEAN

     P <--> The integral from -infinity to X of the normal density.
            Input range: (0,1].

     Q <--> 1-P.
            Input range: (0, 1].
            P + Q = 1.0.

     X < --> Upper limit of integration of the normal-density.
             Input range: ( -infinity, +infinity)

     MEAN <--> The mean of the normal density.
               Input range: (-infinity, +infinity)

     SD <--> Standard Deviation of the normal density.
             Input range: (0, +infinity).

     STATUS <-- 0 if calculation completed correctly
               -I if input parameter number I is out of range
                1 if answer appears to be lower than lowest
                  search bound
                2 if answer appears to be higher than greatest
                  search bound
                3 if P + Q .ne. 1

     BOUND <-- Undefined if STATUS is 0

               Bound exceeded by parameter number I if STATUS
               is negative.

               Lower search bound if STATUS is 1.

               Upper search bound if STATUS is 2.


                              Method




     A slightly modified version of ANORM from

     Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
     Package of Special Function Routines and Test Drivers"
     acm Transactions on Mathematical Software. 19, 22-32.

     is used to calulate the  cumulative standard normal distribution.

     The rational functions from pages  90-95  of Kennedy and Gentle,
     Statistical  Computing,  Marcel  Dekker, NY,  1980 are  used  as
     starting values to Newton's Iterations which compute the inverse
     standard normal.  Therefore no  searches  are necessary for  any
     parameter.

     For X < -15, the asymptotic expansion for the normal is used  as
     the starting value in finding the inverse standard normal.
     This is formula 26.2.12 of Abramowitz and Stegun.


                              Note


      The normal density is proportional to
      exp( - 0.5 * (( X - MEAN)/SD)**2)

**********************************************************************/
{
static int K1 = 1;
static double z,pq;
/*
     ..
     .. Executable Statements ..
*/
/*
     Check arguments
*/
    *status = 0;
    if(!(*which < 1 || *which > 4)) goto S30;
    if(!(*which < 1)) goto S10;
    *bound = 1.0e0;
    goto S20;
S10:
    *bound = 4.0e0;
S20:
    *status = -1;
    return;
S30:
    if(*which == 1) goto S70;
/*
     P
*/
    if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
    if(!(*p <= 0.0e0)) goto S40;
    *bound = 0.0e0;
    goto S50;
S40:
    *bound = 1.0e0;
S50:
    *status = -2;
    return;
S70:
S60:
    if(*which == 1) goto S110;
/*
     Q
*/
    if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
    if(!(*q <= 0.0e0)) goto S80;
    *bound = 0.0e0;
    goto S90;
S80:
    *bound = 1.0e0;
S90:
    *status = -3;
    return;
S110:
S100:
    if(*which == 1) goto S150;
/*
     P + Q
*/
    pq = *p+*q;
    if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S140;
    if(!(pq < 0.0e0)) goto S120;
    *bound = 0.0e0;
    goto S130;
S120:
    *bound = 1.0e0;
S130:
    *status = 3;
    return;
S150:
S140:
    if(*which == 4) goto S170;
/*
     SD
*/
    if(!(*sd <= 0.0e0)) goto S160;
    *bound = 0.0e0;
    *status = -6;
    return;
S170:
S160:
/*
     Calculate ANSWERS
*/
    if(1 == *which) {
/*
     Computing P
*/
        z = (*x-*mean)/ *sd;
        cumnor(&z,p,q);
    }
    else if(2 == *which) {
/*
     Computing X
*/
        z = dinvnr(p,q);
        *x = *sd*z+*mean;
    }
    else if(3 == *which) {
/*
     Computing the MEAN
*/
        z = dinvnr(p,q);
        *mean = *x-*sd*z;
    }
    else if(4 == *which) {
/*
     Computing SD
*/
        z = dinvnr(p,q);
        *sd = (*x-*mean)/z;
    }
    return;
}

double spmpar(int *i)
/*
-----------------------------------------------------------------------
 
     SPMPAR PROVIDES THE SINGLE PRECISION MACHINE CONSTANTS FOR
     THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT
     I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE
     SINGLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND
     ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN
 
        SPMPAR(1) = B**(1 - M), THE MACHINE PRECISION,
 
        SPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE,
 
        SPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE.
 
-----------------------------------------------------------------------
     WRITTEN BY
        ALFRED H. MORRIS, JR.
        NAVAL SURFACE WARFARE CENTER
        DAHLGREN VIRGINIA
-----------------------------------------------------------------------
-----------------------------------------------------------------------
     MODIFIED BY BARRY W. BROWN TO RETURN DOUBLE PRECISION MACHINE
     CONSTANTS FOR THE COMPUTER BEING USED.  THIS MODIFICATION WAS
     MADE AS PART OF CONVERTING BRATIO TO DOUBLE PRECISION
-----------------------------------------------------------------------
*/
{
static int K1 = 4;
static int K2 = 8;
static int K3 = 9;
static int K4 = 10;
static double spmpar,b,binv,bm1,one,w,z;
static int emax,emin,ibeta,m;
/*
     ..
     .. Executable Statements ..
*/
    if(*i > 1) goto S10;
    b = ipmpar(&K1);
    m = ipmpar(&K2);
    spmpar = pow(b,(double)(1-m));
    return spmpar;
S10:
    if(*i > 2) goto S20;
    b = ipmpar(&K1);
    emin = ipmpar(&K3);
    one = 1.0;
    binv = one/b;
    w = pow(b,(double)(emin+2));
    spmpar = w*binv*binv*binv;
    return spmpar;
S20:
    ibeta = ipmpar(&K1);
    m = ipmpar(&K2);
    emax = ipmpar(&K4);
    b = ibeta;
    bm1 = ibeta-1;
    one = 1.0;
    z = pow(b,(double)(m-1));
    w = ((z-one)*b+bm1)/(b*z);
    z = pow(b,(double)(emax-2));
    spmpar = w*z*b*b;
    return spmpar;
}


void cumnor(double *arg,double *result,double *ccum)
/*
**********************************************************************
 
     void cumnor(double *arg,double *result,double *ccum)
 
 
                              Function
 
 
     Computes the cumulative  of    the  normal   distribution,   i.e.,
     the integral from -infinity to x of
          (1/sqrt(2*pi)) exp(-u*u/2) du
 
     X --> Upper limit of integration.
                                        X is DOUBLE PRECISION
 
     RESULT <-- Cumulative normal distribution.
                                        RESULT is DOUBLE PRECISION
 
     CCUM <-- Compliment of Cumulative normal distribution.
                                        CCUM is DOUBLE PRECISION
 
     Renaming of function ANORM from:

     Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
     Package of Special Function Routines and Test Drivers"
     acm Transactions on Mathematical Software. 19, 22-32.

     with slight modifications to return ccum and to deal with
     machine constants.
 
**********************************************************************
  Original Comments:
------------------------------------------------------------------
 
 This function evaluates the normal distribution function:
 
                              / x
                     1       |       -t*t/2
          P(x) = ----------- |      e       dt
                 sqrt(2 pi)  |
                             /-oo
 
   The main computation evaluates near-minimax approximations
   derived from those in "Rational Chebyshev approximations for
   the error function" by W. J. Cody, Math. Comp., 1969, 631-637.
   This transportable program uses rational functions that
   theoretically approximate the normal distribution function to
   at least 18 significant decimal digits.  The accuracy achieved
   depends on the arithmetic system, the compiler, the intrinsic
   functions, and proper selection of the machine-dependent
   constants.
 
*******************************************************************
*******************************************************************
 
 Explanation of machine-dependent constants.
 
   MIN   = smallest machine representable number.
 
   EPS   = argument below which anorm(x) may be represented by
           0.5  and above which  x*x  will not underflow.
           A conservative value is the largest machine number X
           such that   1.0 + X = 1.0   to machine precision.
*******************************************************************
*******************************************************************
 
 Error returns
 
  The program returns  ANORM = 0     for  ARG .LE. XLOW.
 
 
 Intrinsic functions required are:
 
     ABS, AINT, EXP
 
 
  Author: W. J. Cody
          Mathematics and Computer Science Division
          Argonne National Laboratory
          Argonne, IL 60439
 
  Latest modification: March 15, 1992
 
------------------------------------------------------------------
*/
{
static double a[5] = {
    2.2352520354606839287e00,1.6102823106855587881e02,1.0676894854603709582e03,
    1.8154981253343561249e04,6.5682337918207449113e-2
};
static double b[4] = {
    4.7202581904688241870e01,9.7609855173777669322e02,1.0260932208618978205e04,
    4.5507789335026729956e04
};
static double c[9] = {
    3.9894151208813466764e-1,8.8831497943883759412e00,9.3506656132177855979e01,
    5.9727027639480026226e02,2.4945375852903726711e03,6.8481904505362823326e03,
    1.1602651437647350124e04,9.8427148383839780218e03,1.0765576773720192317e-8
};
static double d[8] = {
    2.2266688044328115691e01,2.3538790178262499861e02,1.5193775994075548050e03,
    6.4855582982667607550e03,1.8615571640885098091e04,3.4900952721145977266e04,
    3.8912003286093271411e04,1.9685429676859990727e04
};
static double half = 0.5e0;
static double p[6] = {
    2.1589853405795699e-1,1.274011611602473639e-1,2.2235277870649807e-2,
    1.421619193227893466e-3,2.9112874951168792e-5,2.307344176494017303e-2
};
static double one = 1.0e0;
static double q[5] = {
    1.28426009614491121e00,4.68238212480865118e-1,6.59881378689285515e-2,
    3.78239633202758244e-3,7.29751555083966205e-5
};
static double sixten = 1.60e0;
static double sqrpi = 3.9894228040143267794e-1;
static double thrsh = 0.66291e0;
static double root32 = 5.656854248e0;
static double zero = 0.0e0;
static int K1 = 1;
static int K2 = 2;
static int i;
static double del,eps,temp,x,xden,xnum,y,xsq,min;
/*
------------------------------------------------------------------
  Machine dependent constants
------------------------------------------------------------------
*/
    eps = spmpar(&K1)*0.5e0;
    min = spmpar(&K2);
    x = *arg;
    y = fabs(x);
    if(y <= thrsh) {
/*
------------------------------------------------------------------
  Evaluate  anorm  for  |X| <= 0.66291
------------------------------------------------------------------
*/
        xsq = zero;
        if(y > eps) xsq = x*x;
        xnum = a[4]*xsq;
        xden = xsq;
        for(i=0; i<3; i++) {
            xnum = (xnum+a[i])*xsq;
            xden = (xden+b[i])*xsq;
        }
        *result = x*(xnum+a[3])/(xden+b[3]);
        temp = *result;
        *result = half+temp;
        *ccum = half-temp;
    }
/*
------------------------------------------------------------------
  Evaluate  anorm  for 0.66291 <= |X| <= sqrt(32)
------------------------------------------------------------------
*/
    else if(y <= root32) {
        xnum = c[8]*y;
        xden = y;
        for(i=0; i<7; i++) {
            xnum = (xnum+c[i])*y;
            xden = (xden+d[i])*y;
        }
        *result = (xnum+c[7])/(xden+d[7]);
        xsq = fifdint(y*sixten)/sixten;
        del = (y-xsq)*(y+xsq);
        *result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
        *ccum = one-*result;
        if(x > zero) {
            temp = *result;
            *result = *ccum;
            *ccum = temp;
        }
    }
/*
------------------------------------------------------------------
  Evaluate  anorm  for |X| > sqrt(32)
------------------------------------------------------------------
*/
    else  {
        *result = zero;
        xsq = one/(x*x);
        xnum = p[5]*xsq;
        xden = xsq;
        for(i=0; i<4; i++) {
            xnum = (xnum+p[i])*xsq;
            xden = (xden+q[i])*xsq;
        }
        *result = xsq*(xnum+p[4])/(xden+q[4]);
        *result = (sqrpi-*result)/y;
        xsq = fifdint(x*sixten)/sixten;
        del = (x-xsq)*(x+xsq);
        *result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
        *ccum = one-*result;
        if(x > zero) {
            temp = *result;
            *result = *ccum;
            *ccum = temp;
        }
    }
    if(*result < min) *result = 0.0e0;
/*
------------------------------------------------------------------
  Fix up for negative argument, erf, etc.
------------------------------------------------------------------
----------Last card of ANORM ----------
*/
    if(*ccum < min) *ccum = 0.0e0;
}

double dinvnr(double *p,double *q)
/*
**********************************************************************
 
     double dinvnr(double *p,double *q)
     Double precision NoRmal distribution INVerse
 
 
                              Function
 
 
     Returns X  such that CUMNOR(X)  =   P,  i.e., the  integral from -
     infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P
 
 
                              Arguments
 
 
     P --> The probability whose normal deviate is sought.
                    P is DOUBLE PRECISION
 
     Q --> 1-P
                    P is DOUBLE PRECISION
 
 
                              Method
 
 
     The  rational   function   on  page 95    of Kennedy  and  Gentle,
     Statistical Computing, Marcel Dekker, NY , 1980 is used as a start
     value for the Newton method of finding roots.
 
 
                              Note
 
 
     If P or Q .lt. machine EPS returns +/- DINVNR(EPS)
 
**********************************************************************
*/
{
#define maxit 100
#define eps 1.0e-13
#define r2pi 0.3989422804014326e0
#define nhalf -0.5e0
#define dennor(x) (r2pi*exp(nhalf*(x)*(x)))
static double dinvnr,strtx,xcur,cum,ccum,pp,dx;
static int i;
static unsigned long qporq;
/*
     ..
     .. Executable Statements ..
*/
/*
     FIND MINIMUM OF P AND Q
*/
    qporq = *p <= *q;
    if(!qporq) goto S10;
    pp = *p;
    goto S20;
S10:
    pp = *q;
S20:
/*
     INITIALIZATION STEP
*/
    strtx = stvaln(&pp);
    xcur = strtx;
/*
     NEWTON INTERATIONS
*/
    for(i=1; i<=maxit; i++) {
        cumnor(&xcur,&cum,&ccum);
        dx = (cum-pp)/dennor(xcur);
        xcur -= dx;
        if(fabs(dx/xcur) < eps) goto S40;
    }
    dinvnr = strtx;
/*
     IF WE GET HERE, NEWTON HAS FAILED
*/
    if(!qporq) dinvnr = -dinvnr;
    return dinvnr;
S40:
/*
     IF WE GET HERE, NEWTON HAS SUCCEDED
*/
    dinvnr = xcur;
    if(!qporq) dinvnr = -dinvnr;
    return dinvnr;
#undef maxit
#undef eps
#undef r2pi
#undef nhalf
#undef dennor
}


double stvaln(double *p)
/*
**********************************************************************
 
     double stvaln(double *p)
                    STarting VALue for Neton-Raphon
                calculation of Normal distribution Inverse
 
 
                              Function
 
 
     Returns X  such that CUMNOR(X)  =   P,  i.e., the  integral from -
     infinity to X of (1/SQRT(2*PI)) EXP(-U*U/2) dU is P
 
 
                              Arguments
 
 
     P --> The probability whose normal deviate is sought.
                    P is DOUBLE PRECISION
 
 
                              Method
 
 
     The  rational   function   on  page 95    of Kennedy  and  Gentle,
     Statistical Computing, Marcel Dekker, NY , 1980.
 
**********************************************************************
*/
{
static double xden[5] = {
    0.993484626060e-1,0.588581570495e0,0.531103462366e0,0.103537752850e0,
    0.38560700634e-2
};
static double xnum[5] = {
    -0.322232431088e0,-1.000000000000e0,-0.342242088547e0,-0.204231210245e-1,
    -0.453642210148e-4
};
static int K1 = 5;
static double stvaln,sign,y,z;
/*
     ..
     .. Executable Statements ..
*/
    if(!(*p <= 0.5e0)) goto S10;
    sign = -1.0e0;
    z = *p;
    goto S20;
S10:
    sign = 1.0e0;
    z = 1.0e0-*p;
S20:
    y = sqrt(-(2.0e0*log(z)));
    stvaln = y+devlpl(xnum,&K1,&y)/devlpl(xden,&K1,&y);
    stvaln = sign*stvaln;
    return stvaln;
}

double devlpl(double a[],int *n,double *x)
/*
**********************************************************************
 
     double devlpl(double a[],int *n,double *x)
              Double precision EVALuate a PoLynomial at X
 
 
                              Function
 
 
     returns
          A(1) + A(2)*X + ... + A(N)*X**(N-1)
 
 
                              Arguments
 
 
     A --> Array of coefficients of the polynomial.
                                        A is DOUBLE PRECISION(N)
 
     N --> Length of A, also degree of polynomial - 1.
                                        N is INTEGER
 
     X --> Point at which the polynomial is to be evaluated.
                                        X is DOUBLE PRECISION
 
**********************************************************************
*/
{
static double devlpl,term;
static int i;
/*
     ..
     .. Executable Statements ..
*/
    term = a[*n-1];
    for(i= *n-1-1; i>=0; i--) term = a[i]+term**x;
    devlpl = term;
    return devlpl;
}


double gamln(double *a)
/*
-----------------------------------------------------------------------
            EVALUATION OF LN(GAMMA(A)) FOR POSITIVE A
-----------------------------------------------------------------------
     WRITTEN BY ALFRED H. MORRIS
          NAVAL SURFACE WARFARE CENTER
          DAHLGREN, VIRGINIA
--------------------------
     D = 0.5*(LN(2*PI) - 1)
--------------------------
*/
{
static double c0 = .833333333333333e-01;
static double c1 = -.277777777760991e-02;
static double c2 = .793650666825390e-03;
static double c3 = -.595202931351870e-03;
static double c4 = .837308034031215e-03;
static double c5 = -.165322962780713e-02;
static double d = .418938533204673e0;
static double gamln,t,w;
static int i,n;
static double T1;
/*
     ..
     .. Executable Statements ..
*/
    if(*a > 0.8e0) goto S10;
    gamln = gamln1(a)-log(*a);
    return gamln;
S10:
    if(*a > 2.25e0) goto S20;
    t = *a-0.5e0-0.5e0;
    gamln = gamln1(&t);
    return gamln;
S20:
    if(*a >= 10.0e0) goto S40;
    n = (long)(*a - 1.25e0);
    t = *a;
    w = 1.0e0;
    for(i=1; i<=n; i++) {
        t -= 1.0e0;
        w = t*w;
    }
    T1 = t-1.0e0;
    gamln = gamln1(&T1)+log(w);
    return gamln;
S40:
    t = pow(1.0e0/ *a,2.0);
    w = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/ *a;
    gamln = d+w+(*a-0.5e0)*(log(*a)-1.0e0);
    return gamln;
}

double gamln1(double *a)
/*
-----------------------------------------------------------------------
     EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 .LE. A .LE. 1.25
-----------------------------------------------------------------------
*/
{
static double p0 = .577215664901533e+00;
static double p1 = .844203922187225e+00;
static double p2 = -.168860593646662e+00;
static double p3 = -.780427615533591e+00;
static double p4 = -.402055799310489e+00;
static double p5 = -.673562214325671e-01;
static double p6 = -.271935708322958e-02;
static double q1 = .288743195473681e+01;
static double q2 = .312755088914843e+01;
static double q3 = .156875193295039e+01;
static double q4 = .361951990101499e+00;
static double q5 = .325038868253937e-01;
static double q6 = .667465618796164e-03;
static double r0 = .422784335098467e+00;
static double r1 = .848044614534529e+00;
static double r2 = .565221050691933e+00;
static double r3 = .156513060486551e+00;
static double r4 = .170502484022650e-01;
static double r5 = .497958207639485e-03;
static double s1 = .124313399877507e+01;
static double s2 = .548042109832463e+00;
static double s3 = .101552187439830e+00;
static double s4 = .713309612391000e-02;
static double s5 = .116165475989616e-03;
static double gamln1,w,x;
/*
     ..
     .. Executable Statements ..
*/
    if(*a >= 0.6e0) goto S10;
    w = ((((((p6**a+p5)**a+p4)**a+p3)**a+p2)**a+p1)**a+p0)/((((((q6**a+q5)**a+
      q4)**a+q3)**a+q2)**a+q1)**a+1.0e0);
    gamln1 = -(*a*w);
    return gamln1;
S10:
    x = *a-0.5e0-0.5e0;
    w = (((((r5*x+r4)*x+r3)*x+r2)*x+r1)*x+r0)/(((((s5*x+s4)*x+s3)*x+s2)*x+s1)*x
      +1.0e0);
    gamln1 = x*w;
    return gamln1;
}



int ipmpar(int*);
/*
-----------------------------------------------------------------------
 
     IPMPAR PROVIDES THE INTEGER MACHINE CONSTANTS FOR THE COMPUTER
     THAT IS USED. IT IS ASSUMED THAT THE ARGUMENT I IS AN INTEGER
     HAVING ONE OF THE VALUES 1-10. IPMPAR(I) HAS THE VALUE ...
 
  INTEGERS.
 
     ASSUME INTEGERS ARE REPRESENTED IN THE N-DIGIT, BASE-A FORM
 
               SIGN ( X(N-1)*A**(N-1) + ... + X(1)*A + X(0) )
 
               WHERE 0 .LE. X(I) .LT. A FOR I=0,...,N-1.
 
     IPMPAR(1) = A, THE BASE.
 
     IPMPAR(2) = N, THE NUMBER OF BASE-A DIGITS.
 
     IPMPAR(3) = A**N - 1, THE LARGEST MAGNITUDE.
 
  FLOATING-POINT NUMBERS.
 
     IT IS ASSUMED THAT THE SINGLE AND DOUBLE PRECISION FLOATING
     POINT ARITHMETICS HAVE THE SAME BASE, SAY B, AND THAT THE
     NONZERO NUMBERS ARE REPRESENTED IN THE FORM
 
               SIGN (B**E) * (X(1)/B + ... + X(M)/B**M)
 
               WHERE X(I) = 0,1,...,B-1 FOR I=1,...,M,
               X(1) .GE. 1, AND EMIN .LE. E .LE. EMAX.
 
     IPMPAR(4) = B, THE BASE.
 
  SINGLE-PRECISION
 
     IPMPAR(5) = M, THE NUMBER OF BASE-B DIGITS.
 
     IPMPAR(6) = EMIN, THE SMALLEST EXPONENT E.
 
     IPMPAR(7) = EMAX, THE LARGEST EXPONENT E.
 
  DOUBLE-PRECISION
 
     IPMPAR(8) = M, THE NUMBER OF BASE-B DIGITS.
 
     IPMPAR(9) = EMIN, THE SMALLEST EXPONENT E.
 
     IPMPAR(10) = EMAX, THE LARGEST EXPONENT E.
 
-----------------------------------------------------------------------
 
     TO DEFINE THIS FUNCTION FOR THE COMPUTER BEING USED REMOVE
     THE COMMENT DELIMITORS FROM THE DEFINITIONS DIRECTLY BELOW THE NAME
     OF THE MACHINE
 
-----------------------------------------------------------------------
 
     IPMPAR IS AN ADAPTATION OF THE FUNCTION I1MACH, WRITTEN BY
     P.A. FOX, A.D. HALL, AND N.L. SCHRYER (BELL LABORATORIES).
     IPMPAR WAS FORMED BY A.H. MORRIS (NSWC). THE CONSTANTS ARE
     FROM BELL LABORATORIES, NSWC, AND OTHER SOURCES.
 
-----------------------------------------------------------------------
     .. Scalar Arguments ..
*/
int ipmpar(int *i)
{
static int imach[11];
static int ipmpar;
/*     MACHINE CONSTANTS FOR AMDAHL MACHINES. */
/*
   imach[1] = 2;
   imach[2] = 31;
   imach[3] = 2147483647;
   imach[4] = 16;
   imach[5] = 6;
   imach[6] = -64;
   imach[7] = 63;
   imach[8] = 14;
   imach[9] = -64;
   imach[10] = 63;
*/
/*     MACHINE CONSTANTS FOR THE AT&T 3B SERIES, AT&T
       PC 7300, AND AT&T 6300. */
/*
   imach[1] = 2;
   imach[2] = 31;
   imach[3] = 2147483647;
   imach[4] = 2;
   imach[5] = 24;
   imach[6] = -125;
   imach[7] = 128;
   imach[8] = 53;
   imach[9] = -1021;
   imach[10] = 1024;
*/
/*     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. */
/*
   imach[1] = 2;
   imach[2] = 33;
   imach[3] = 8589934591;
   imach[4] = 2;
   imach[5] = 24;
   imach[6] = -256;
   imach[7] = 255;
   imach[8] = 60;
   imach[9] = -256;
   imach[10] = 255;
*/
/*     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. */
/*
   imach[1] = 2;
   imach[2] = 39;
   imach[3] = 549755813887;
   imach[4] = 8;
   imach[5] = 13;
   imach[6] = -50;
   imach[7] = 76;
   imach[8] = 26;
   imach[9] = -50;
   imach[10] = 76;
*/
/*     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. */
/*
   imach[1] = 2;
   imach[2] = 39;
   imach[3] = 549755813887;
   imach[4] = 8;
   imach[5] = 13;
   imach[6] = -50;
   imach[7] = 76;
   imach[8] = 26;
   imach[9] = -32754;
   imach[10] = 32780;
*/
/*     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES
       60 BIT ARITHMETIC, AND THE CDC CYBER 995 64 BIT
       ARITHMETIC (NOS OPERATING SYSTEM). */
/*
   imach[1] = 2;
   imach[2] = 48;
   imach[3] = 281474976710655;
   imach[4] = 2;
   imach[5] = 48;
   imach[6] = -974;
   imach[7] = 1070;
   imach[8] = 95;
   imach[9] = -926;
   imach[10] = 1070;
*/
/*     MACHINE CONSTANTS FOR THE CDC CYBER 995 64 BIT
       ARITHMETIC (NOS/VE OPERATING SYSTEM). */
/*
   imach[1] = 2;
   imach[2] = 63;
   imach[3] = 9223372036854775807;
   imach[4] = 2;
   imach[5] = 48;
   imach[6] = -4096;
   imach[7] = 4095;
   imach[8] = 96;
   imach[9] = -4096;
   imach[10] = 4095;
*/
/*     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. */
/*
   imach[1] = 2;
   imach[2] = 63;
   imach[3] = 9223372036854775807;
   imach[4] = 2;
   imach[5] = 47;
   imach[6] = -8189;
   imach[7] = 8190;
   imach[8] = 94;
   imach[9] = -8099;
   imach[10] = 8190;
*/
/*     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200. */
/*
   imach[1] = 2;
   imach[2] = 15;
   imach[3] = 32767;
   imach[4] = 16;
   imach[5] = 6;
   imach[6] = -64;
   imach[7] = 63;
   imach[8] = 14;
   imach[9] = -64;
   imach[10] = 63;
*/
/*     MACHINE CONSTANTS FOR THE HARRIS 220. */
/*
   imach[1] = 2;
   imach[2] = 23;
   imach[3] = 8388607;
   imach[4] = 2;
   imach[5] = 23;
   imach[6] = -127;
   imach[7] = 127;
   imach[8] = 38;
   imach[9] = -127;
   imach[10] = 127;
*/
/*     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000
       AND DPS 8/70 SERIES. */
/*
   imach[1] = 2;
   imach[2] = 35;
   imach[3] = 34359738367;
   imach[4] = 2;
   imach[5] = 27;
   imach[6] = -127;
   imach[7] = 127;
   imach[8] = 63;
   imach[9] = -127;
   imach[10] = 127;
*/
/*     MACHINE CONSTANTS FOR THE HP 2100
       3 WORD DOUBLE PRECISION OPTION WITH FTN4 */
/*
   imach[1] = 2;
   imach[2] = 15;
   imach[3] = 32767;
   imach[4] = 2;
   imach[5] = 23;
   imach[6] = -128;
   imach[7] = 127;
   imach[8] = 39;
   imach[9] = -128;
   imach[10] = 127;
*/
/*     MACHINE CONSTANTS FOR THE HP 2100
       4 WORD DOUBLE PRECISION OPTION WITH FTN4 */
/*
   imach[1] = 2;
   imach[2] = 15;
   imach[3] = 32767;
   imach[4] = 2;
   imach[5] = 23;
   imach[6] = -128;
   imach[7] = 127;
   imach[8] = 55;
   imach[9] = -128;
   imach[10] = 127;
*/
/*     MACHINE CONSTANTS FOR THE HP 9000. */
/*
   imach[1] = 2;
   imach[2] = 31;
   imach[3] = 2147483647;
   imach[4] = 2;
   imach[5] = 24;
   imach[6] = -126;
   imach[7] = 128;
   imach[8] = 53;
   imach[9] = -1021;
   imach[10] = 1024;
*/
/*     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
       THE ICL 2900, THE ITEL AS/6, THE XEROX SIGMA
       5/7/9 AND THE SEL SYSTEMS 85/86. */
/*
   imach[1] = 2;
   imach[2] = 31;
   imach[3] = 2147483647;
   imach[4] = 16;
   imach[5] = 6;
   imach[6] = -64;
   imach[7] = 63;
   imach[8] = 14;
   imach[9] = -64;
   imach[10] = 63;
*/
/*     MACHINE CONSTANTS FOR THE IBM PC. */
/*
   imach[1] = 2;
   imach[2] = 31;
   imach[3] = 2147483647;
   imach[4] = 2;
   imach[5] = 24;
   imach[6] = -125;
   imach[7] = 128;
   imach[8] = 53;
   imach[9] = -1021;
   imach[10] = 1024;
*/
/*     MACHINE CONSTANTS FOR THE MACINTOSH II - ABSOFT
       MACFORTRAN II. */
/*
   imach[1] = 2;
   imach[2] = 31;
   imach[3] = 2147483647;
   imach[4] = 2;
   imach[5] = 24;
   imach[6] = -125;
   imach[7] = 128;
   imach[8] = 53;
   imach[9] = -1021;
   imach[10] = 1024;
*/
/*     MACHINE CONSTANTS FOR THE MICROVAX - VMS FORTRAN. */
/*
   imach[1] = 2;
   imach[2] = 31;
   imach[3] = 2147483647;
   imach[4] = 2;
   imach[5] = 24;
   imach[6] = -127;
   imach[7] = 127;
   imach[8] = 56;
   imach[9] = -127;
   imach[10] = 127;
*/
/*     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). */
/*
   imach[1] = 2;
   imach[2] = 35;
   imach[3] = 34359738367;
   imach[4] = 2;
   imach[5] = 27;
   imach[6] = -128;
   imach[7] = 127;
   imach[8] = 54;
   imach[9] = -101;
   imach[10] = 127;
*/
/*     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). */
/*
   imach[1] = 2;
   imach[2] = 35;
   imach[3] = 34359738367;
   imach[4] = 2;
   imach[5] = 27;
   imach[6] = -128;
   imach[7] = 127;
   imach[8] = 62;
   imach[9] = -128;
   imach[10] = 127;
*/
/*     MACHINE CONSTANTS FOR THE PDP-11 FORTRAN SUPPORTING
       32-BIT INTEGER ARITHMETIC. */
/*
   imach[1] = 2;
   imach[2] = 31;
   imach[3] = 2147483647;
   imach[4] = 2;
   imach[5] = 24;
   imach[6] = -127;
   imach[7] = 127;
   imach[8] = 56;
   imach[9] = -127;
   imach[10] = 127;
*/
/*     MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000. */
/*
   imach[1] = 2;
   imach[2] = 31;
   imach[3] = 2147483647;
   imach[4] = 2;
   imach[5] = 24;
   imach[6] = -125;
   imach[7] = 128;
   imach[8] = 53;
   imach[9] = -1021;
   imach[10] = 1024;
*/
/*     MACHINE CONSTANTS FOR THE SILICON GRAPHICS IRIS-4D
       SERIES (MIPS R3000 PROCESSOR). */
/*
   imach[1] = 2;
   imach[2] = 31;
   imach[3] = 2147483647;
   imach[4] = 2;
   imach[5] = 24;
   imach[6] = -125;
   imach[7] = 128;
   imach[8] = 53;
   imach[9] = -1021;
   imach[10] = 1024;
*/
/*     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T
       3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T
       PC 7300), AND 8087 BASED MICROS (E.G. IBM PC AND AT&T 6300). */

   imach[1] = 2;
   imach[2] = 31;
   imach[3] = 2147483647;
   imach[4] = 2;
   imach[5] = 24;
   imach[6] = -125;
   imach[7] = 128;
   imach[8] = 53;
   imach[9] = -1021;
   imach[10] = 1024;

/*     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. */
/*
   imach[1] = 2;
   imach[2] = 35;
   imach[3] = 34359738367;
   imach[4] = 2;
   imach[5] = 27;
   imach[6] = -128;
   imach[7] = 127;
   imach[8] = 60;
   imach[9] = -1024;
   imach[10] = 1023;
*/
/*     MACHINE CONSTANTS FOR THE VAX 11/780. */
/*
   imach[1] = 2;
   imach[2] = 31;
   imach[3] = 2147483647;
   imach[4] = 2;
   imach[5] = 24;
   imach[6] = -127;
   imach[7] = 127;
   imach[8] = 56;
   imach[9] = -127;
   imach[10] = 127;
*/
    ipmpar = imach[*i];
    return ipmpar;
}


/************************************************************************
                        EVEN MORE STUFF
************************************************************************/

float genunf(float low,float high)
/*
**********************************************************************
     float genunf(float low,float high)
               GeNerate Uniform Real between LOW and HIGH
                              Function
     Generates a real uniformly distributed between LOW and HIGH.
                              Arguments
     low --> Low bound (exclusive) on real value to be generated
     high --> High bound (exclusive) on real value to be generated
**********************************************************************
*/
{
static float genunf;

    if(!(low > high)) goto S10;
    printf("LOW > HIGH in GENUNF: LOW %16.6E HIGH: %16.6E\n",low,high);
    puts("Abort");
    exit(1);
S10:
    genunf = low+(high-low)*ranf();
    return genunf;
}

float gengam(float a,float r)
/*
**********************************************************************
     float gengam(float a,float r)
           GENerates random deviates from GAMma distribution
                              Function
     Generates random deviates from the gamma distribution whose
     density is
          (A**R)/Gamma(R) * X**(R-1) * Exp(-A*X)
                              Arguments
     a --> Location parameter of Gamma distribution
     r --> Shape parameter of Gamma distribution
     CAREFUL: order of parameters is reversed wrt usual notation. mean is r/a; var is r/a^2
                              Method
     Renames SGAMMA from TOMS as slightly modified by BWB to use RANF
     instead of SUNIF.
     For details see:
               (Case R >= 1.0)
               Ahrens, J.H. and Dieter, U.
               Generating Gamma Variates by a
               Modified Rejection Technique.
               Comm. ACM, 25,1 (Jan. 1982), 47 - 54.
     Algorithm GD
               (Case 0.0 <= R <= 1.0)
               Ahrens, J.H. and Dieter, U.
               Computer Methods for Sampling from Gamma,
               Beta, Poisson and Binomial Distributions.
               Computing, 12 (1974), 223-246/
     Adapted algorithm GS.
**********************************************************************
*/
{
static float gengam;

    gengam = sgamma(r);
    gengam /= a;
    return gengam;
}


float sgamma(float a)
/*
**********************************************************************
                                                                      
                                                                      
     (STANDARD-)  G A M M A  DISTRIBUTION                             
                                                                      
                                                                      
**********************************************************************
**********************************************************************
                                                                      
               PARAMETER  A >= 1.0  !                                 
                                                                      
**********************************************************************
                                                                      
     FOR DETAILS SEE:                                                 
                                                                      
               AHRENS, J.H. AND DIETER, U.                            
               GENERATING GAMMA VARIATES BY A                         
               MODIFIED REJECTION TECHNIQUE.                          
               COMM. ACM, 25,1 (JAN. 1982), 47 - 54.                  
                                                                      
     STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER     
                                 (STRAIGHTFORWARD IMPLEMENTATION)     
                                                                      
     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   
     SUNIF.  The argument IR thus goes away.                          
                                                                      
**********************************************************************
                                                                      
               PARAMETER  0.0 < A < 1.0  !                            
                                                                      
**********************************************************************
                                                                      
     FOR DETAILS SEE:                                                 
                                                                      
               AHRENS, J.H. AND DIETER, U.                            
               COMPUTER METHODS FOR SAMPLING FROM GAMMA,              
               BETA, POISSON AND BINOMIAL DISTRIBUTIONS.              
               COMPUTING, 12 (1974), 223 - 246.                       
                                                                      
     (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)    
                                                                      
**********************************************************************
     INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
     OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
     COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
     COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
     COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
     PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
     SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
*/
{
extern float fsign( float num, float sign );
static float q1 = 4.166669E-2;
static float q2 = 2.083148E-2;
static float q3 = 8.01191E-3;
static float q4 = 1.44121E-3;
static float q5 = -7.388E-5;
static float q6 = 2.4511E-4;
static float q7 = 2.424E-4;
static float a1 = 0.3333333;
static float a2 = -0.250003;
static float a3 = 0.2000062;
static float a4 = -0.1662921;
static float a5 = 0.1423657;
static float a6 = -0.1367177;
static float a7 = 0.1233795;
static float e1 = 1.0;
static float e2 = 0.4999897;
static float e3 = 0.166829;
static float e4 = 4.07753E-2;
static float e5 = 1.0293E-2;
static float aa = 0.0;
static float aaa = 0.0;
static float sqrt32 = 5.656854;
static float sgamma,s2,s,d,t,x,u,r,q0,b,si,c,v,q,e,w,p;
    if(a == aa) goto S10;
    if(a < 1.0) goto S120;
/*
     STEP  1:  RECALCULATIONS OF S2,S,D IF A HAS CHANGED
*/
    aa = a;
    s2 = a-0.5;
    s = sqrt(s2);
    d = sqrt32-12.0*s;
S10:
/*
     STEP  2:  T=STANDARD NORMAL DEVIATE,
               X=(S,1/2)-NORMAL DEVIATE.
               IMMEDIATE ACCEPTANCE (I)
*/
    t = snorm();
    x = s+0.5*t;
    sgamma = x*x;
    if(t >= 0.0) return sgamma;
/*
     STEP  3:  U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
*/
    u = ranf();
    if(d*u <= t*t*t) return sgamma;
/*
     STEP  4:  RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
*/
    if(a == aaa) goto S40;
    aaa = a;
    r = 1.0/ a;
    q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r;
/*
               APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
               THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
               C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
*/
    if(a <= 3.686) goto S30;
    if(a <= 13.022) goto S20;
/*
               CASE 3:  A .GT. 13.022
*/
    b = 1.77;
    si = 0.75;
    c = 0.1515/s;
    goto S40;
S20:
/*
               CASE 2:  3.686 .LT. A .LE. 13.022
*/
    b = 1.654+7.6E-3*s2;
    si = 1.68/s+0.275;
    c = 6.2E-2/s+2.4E-2;
    goto S40;
S30:
/*
               CASE 1:  A .LE. 3.686
*/
    b = 0.463+s-0.178*s2;
    si = 1.235;
    c = 0.195/s-7.9E-2+1.6E-2*s;
S40:
/*
     STEP  5:  NO QUOTIENT TEST IF X NOT POSITIVE
*/
    if(x <= 0.0) goto S70;
/*
     STEP  6:  CALCULATION OF V AND QUOTIENT Q
*/
    v = t/(s+s);
    if(fabs(v) <= 0.25) goto S50;
    q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
    goto S60;
S50:
    q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
S60:
/*
     STEP  7:  QUOTIENT ACCEPTANCE (Q)
*/
    if(log(1.0-u) <= q) return sgamma;
S70:
/*
     STEP  8:  E=STANDARD EXPONENTIAL DEVIATE
               U= 0,1 -UNIFORM DEVIATE
               T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
*/
    e = sexpo();
    u = ranf();
    u += (u-1.0);
    t = b+fsign(si*e,u);
/*
     STEP  9:  REJECTION IF T .LT. TAU(1) = -.71874483771719
*/
    if(t < -0.7187449) goto S70;
/*
     STEP 10:  CALCULATION OF V AND QUOTIENT Q
*/
    v = t/(s+s);
    if(fabs(v) <= 0.25) goto S80;
    q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
    goto S90;
S80:
    q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
S90:
/*
     STEP 11:  HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
*/
    if(q <= 0.0) goto S70;
    if(q <= 0.5) goto S100;
    w = exp(q)-1.0;
    goto S110;
S100:
    w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q;
S110:
/*
               IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
*/
    if(c*fabs(u) > w*exp(e-0.5*t*t)) goto S70;
    x = s+0.5*t;
    sgamma = x*x;
    return sgamma;
S120:
/*
     ALTERNATE METHOD FOR PARAMETERS A BELOW 1  (.3678794=EXP(-1.))
*/
    aa = 0.0;
    b = 1.0+0.3678794*a;
S130:
    p = b*ranf();
    if(p >= 1.0) goto S140;
    sgamma = exp(log(p)/ a);
    if(sexpo() < sgamma) goto S130;
    return sgamma;
S140:
    sgamma = -log((b-p)/ a);
    if(sexpo() < (1.0-a)*log(sgamma)) goto S130;
    return sgamma;
}


float snorm(void)
/*
**********************************************************************
                                                                      
                                                                      
     (STANDARD-)  N O R M A L  DISTRIBUTION                           
                                                                      
                                                                      
**********************************************************************
**********************************************************************
                                                                      
     FOR DETAILS SEE:                                                 
                                                                      
               AHRENS, J.H. AND DIETER, U.                            
               EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM             
               SAMPLING FROM THE NORMAL DISTRIBUTION.                 
               MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.          
                                                                      
     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'  
     (M=5) IN THE ABOVE PAPER     (SLIGHTLY MODIFIED IMPLEMENTATION)  
                                                                      
     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   
     SUNIF.  The argument IR thus goes away.                          
                                                                      
**********************************************************************
     THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
     H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
*/
{
static float a[32] = {
    0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,
    0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
    0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
    1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
    1.862732,2.153875
};
static float d[31] = {
    0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
    0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
    0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
    0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
};
static float t[31] = {
    7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
    1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
    2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
    4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
    9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
};
static float h[31] = {
    3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
    4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
    4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
    5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
    8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
};
static long i;
static float snorm,u,s,ustar,aa,w,y,tt;
    u = ranf();
    s = 0.0;
    if(u > 0.5) s = 1.0;
    u += (u-s);
    u = 32.0*u;
    i = (long) (u);
    if(i == 32) i = 31;
    if(i == 0) goto S100;
/*
                                START CENTER
*/
    ustar = u-(float)i;
    aa = *(a+i-1);
S40:
    if(ustar <= *(t+i-1)) goto S60;
    w = (ustar-*(t+i-1))**(h+i-1);
S50:
/*
                                EXIT   (BOTH CASES)
*/
    y = aa+w;
    snorm = y;
    if(s == 1.0) snorm = -y;
    return snorm;
S60:
/*
                                CENTER CONTINUED
*/
    u = ranf();
    w = u*(*(a+i)-aa);
    tt = (0.5*w+aa)*w;
    goto S80;
S70:
    tt = u;
    ustar = ranf();
S80:
    if(ustar > tt) goto S50;
    u = ranf();
    if(ustar >= u) goto S70;
    ustar = ranf();
    goto S40;
S100:
/*
                                START TAIL
*/
    i = 6;
    aa = *(a+31);
    goto S120;
S110:
    aa += *(d+i-1);
    i += 1;
S120:
    u += u;
    if(u < 1.0) goto S110;
    u -= 1.0;
S140:
    w = u**(d+i-1);
    tt = (0.5*w+aa)*w;
    goto S160;
S150:
    tt = u;
S160:
    ustar = ranf();
    if(ustar > tt) goto S50;
    u = ranf();
    if(ustar >= u) goto S150;
    u = ranf();
    goto S140;
}

float fsign( float num, float sign )
/* Transfers sign of argument sign to argument num */
{
if ( ( sign > 0.0 && num < 0.0) || ( sign < 0.0 && num > 0.0) )
    return -num;
else return num;
}



float sexpo(void)
/*
**********************************************************************
                                                                      
                                                                      
     (STANDARD-)  E X P O N E N T I A L   DISTRIBUTION                
                                                                      
                                                                      
**********************************************************************
**********************************************************************
                                                                      
     FOR DETAILS SEE:                                                 
                                                                      
               AHRENS, J.H. AND DIETER, U.                            
               COMPUTER METHODS FOR SAMPLING FROM THE                 
               EXPONENTIAL AND NORMAL DISTRIBUTIONS.                  
               COMM. ACM, 15,10 (OCT. 1972), 873 - 882.               
                                                                      
     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM       
     'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)       
                                                                      
     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   
     SUNIF.  The argument IR thus goes away.                          
                                                                      
**********************************************************************
     Q(N) = SUM(ALOG(2.0)**K/K!)    K=1,..,N ,      THE HIGHEST N
     (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
*/
{
static float q[8] = {
    0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,1.0
};
static long i;
static float sexpo,a,u,ustar,umin;
static float *q1 = q;
    a = 0.0;
    u = ranf();
    goto S30;
S20:
    a += *q1;
S30:
    u += u;
    if(u <= 1.0) goto S20;
    u -= 1.0;
    if(u > *q1) goto S60;
    sexpo = a+u;
    return sexpo;
S60:
    i = 1;
    ustar = ranf();
    umin = ustar;
S70:
    ustar = ranf();
    if(ustar < umin) umin = ustar;
    i += 1;
    if(u > *(q+i-1)) goto S70;
    sexpo = a+umin**q1;
    return sexpo;
}

long mltmod(long a,long s,long m)
/*
**********************************************************************
     long mltmod(long a,long s,long m)
                    Returns (A*S) MOD M
     This is a transcription from Pascal to Fortran of routine
     MULtMod_Decompos from the paper
     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
     with Splitting Facilities." ACM Transactions on Mathematical
     Software, 17:98-111 (1991)
                              Arguments
     a, s, m  -->
**********************************************************************
*/
{
#define h 32768L
static long mltmod,a0,a1,k,p,q,qh,rh;
/*
     H = 2**((b-2)/2) where b = 32 because we are using a 32 bit
      machine. On a different machine recompute H
*/
    if(!(a <= 0 || a >= m || s <= 0 || s >= m)) goto S10;
    puts(" a, m, s out of order in mltmod - ABORT!");
    printf(" a = %12ld s = %12ld m = %12ld\n",a,s,m);
    puts(" mltmod requires: 0 < a < m; 0 < s < m");
    exit(1);
S10:
    if(!(a < h)) goto S20;
    a0 = a;
    p = 0;
    goto S120;
S20:
    a1 = a/h;
    a0 = a-h*a1;
    qh = m/h;
    rh = m-h*qh;
    if(!(a1 >= h)) goto S50;
    a1 -= h;
    k = s/qh;
    p = h*(s-k*qh)-k*rh;
S30:
    if(!(p < 0)) goto S40;
    p += m;
    goto S30;
S40:
    goto S60;
S50:
    p = 0;
S60:
/*
     P = (A2*S*H)MOD M
*/
    if(!(a1 != 0)) goto S90;
    q = m/a1;
    k = s/q;
    p -= (k*(m-a1*q));
    if(p > 0) p -= m;
    p += (a1*(s-k*q));
S70:
    if(!(p < 0)) goto S80;
    p += m;
    goto S70;
S90:
S80:
    k = p/qh;
/*
     P = ((A2*H + A1)*S)MOD M
*/
    p = h*(p-k*qh)-k*rh;
S100:
    if(!(p < 0)) goto S110;
    p += m;
    goto S100;
S120:
S110:
    if(!(a0 != 0)) goto S150;
/*
     P = ((A2*H + A1)*H*S)MOD M
*/
    q = m/a0;
    k = s/q;
    p -= (k*(m-a0*q));
    if(p > 0) p -= m;
    p += (a0*(s-k*q));
S130:
    if(!(p < 0)) goto S140;
    p += m;
    goto S130;
S150:
S140:
    mltmod = p;
    return mltmod;
#undef h
}

float ranf(void)
/*
**********************************************************************
     float ranf(void)
                RANDom number generator as a Function
     Returns a random floating point number from a uniform distribution
     over 0 - 1 (endpoints of this interval are not returned) using the
     current generator
     This is a transcription from Pascal to Fortran of routine
     Uniform_01 from the paper
     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
     with Splitting Facilities." ACM Transactions on Mathematical
     Software, 17:98-111 (1991)
**********************************************************************
*/
{
static float ranf;
/*
     4.656613057E-10 is 1/M1  M1 is set in a data statement in IGNLGI
      and is currently 2147483563. If M1 changes, change this also.
*/
    ranf = ignlgi()*4.656613057E-10;
    return ranf;
}

void gscgn(long getset,long *g)
/*
**********************************************************************
     void gscgn(long getset,long *g)
                         Get/Set GeNerator
     Gets or returns in G the number of the current generator
                              Arguments
     getset --> 0 Get
                1 Set
     g <-- Number of the current random number generator (1..32)
**********************************************************************
*/
{
#define numg 32L
static long curntg = 1;
    if(getset == 0) *g = curntg;
    else  {
        if(*g < 0 || *g > numg) {
            puts(" Generator number out of range in GSCGN");
            exit(0);
        }
        curntg = *g;
    }
#undef numg
}

void gsrgs(long getset,long *qvalue)
/*
**********************************************************************
     void gsrgs(long getset,long *qvalue)
               Get/Set Random Generators Set
     Gets or sets whether random generators set (initialized).
     Initially (data statement) state is not set
     If getset is 1 state is set to qvalue
     If getset is 0 state returned in qvalue
**********************************************************************
*/
{
static long qinit = 0;

    if(getset == 0) *qvalue = qinit;
    else qinit = *qvalue;
}

void gssst(long getset,long *qset)
/*
**********************************************************************
     void gssst(long getset,long *qset)
          Get or Set whether Seed is Set
     Initialize to Seed not Set
     If getset is 1 sets state to Seed Set
     If getset is 0 returns T in qset if Seed Set
     Else returns F in qset
**********************************************************************
*/
{
static long qstate = 0;
    if(getset != 0) qstate = 1;
    else  *qset = qstate;
}

void setall(long iseed1,long iseed2)
/*
**********************************************************************
     void setall(long iseed1,long iseed2)
               SET ALL random number generators
     Sets the initial seed of generator 1 to ISEED1 and ISEED2. The
     initial seeds of the other generators are set accordingly, and
     all generators states are set to these seeds.
     This is a transcription from Pascal to Fortran of routine
     Set_Initial_Seed from the paper
     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
     with Splitting Facilities." ACM Transactions on Mathematical
     Software, 17:98-111 (1991)
                              Arguments
     iseed1 -> First of two integer seeds
     iseed2 -> Second of two integer seeds
**********************************************************************
*/
{
#define numg 32L
extern void gsrgs(long getset,long *qvalue);
extern void gssst(long getset,long *qset);
extern void gscgn(long getset,long *g);
extern long Xm1,Xm2,Xa1vw,Xa2vw,Xig1[],Xig2[];
static long T1;
static long g,ocgn;
static long qrgnin;
    T1 = 1;
/*
     TELL IGNLGI, THE ACTUAL NUMBER GENERATOR, THAT THIS ROUTINE
      HAS BEEN CALLED.
*/
    gssst(1,&T1);
    gscgn(0L,&ocgn);
/*
     Initialize Common Block if Necessary
*/
    gsrgs(0L,&qrgnin);
    if(!qrgnin) inrgcm();
    *Xig1 = iseed1; 
    *Xig2 = iseed2;
    initgn(-1L);
    for(g=2; g<=numg; g++) {
        *(Xig1+g-1) = mltmod(Xa1vw,*(Xig1+g-2),Xm1);
        *(Xig2+g-1) = mltmod(Xa2vw,*(Xig2+g-2),Xm2);
        gscgn(1L,&g);
        initgn(-1L);
    }
    gscgn(1L,&ocgn);
#undef numg
}

void initgn(long isdtyp)
/*
**********************************************************************
     void initgn(long isdtyp)
          INIT-ialize current G-e-N-erator
     Reinitializes the state of the current generator
     This is a transcription from Pascal to Fortran of routine
     Init_Generator from the paper
     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
     with Splitting Facilities." ACM Transactions on Mathematical
     Software, 17:98-111 (1991)
                              Arguments
     isdtyp -> The state to which the generator is to be set
          isdtyp = -1  => sets the seeds to their initial value
          isdtyp =  0  => sets the seeds to the first value of
                          the current block
          isdtyp =  1  => sets the seeds to the first value of
                          the next block
**********************************************************************
*/
{
#define numg 32L
extern void gsrgs(long getset,long *qvalue);
extern void gscgn(long getset,long *g);
extern long Xm1,Xm2,Xa1w,Xa2w,Xig1[],Xig2[],Xlg1[],Xlg2[],Xcg1[],Xcg2[];
static long g;
static long qrgnin;
/*
     Abort unless random number generator initialized
*/
    gsrgs(0L,&qrgnin);
    if(qrgnin) goto S10;
    printf(
      " INITGN called before random number generator  initialized -- abort!\n");
    exit(1);
S10:
    gscgn(0L,&g);
    if(-1 != isdtyp) goto S20;
    *(Xlg1+g-1) = *(Xig1+g-1);
    *(Xlg2+g-1) = *(Xig2+g-1);
    goto S50;
S20:
    if(0 != isdtyp) goto S30;
    goto S50;
S30:
/*
     do nothing
*/
    if(1 != isdtyp) goto S40;
    *(Xlg1+g-1) = mltmod(Xa1w,*(Xlg1+g-1),Xm1);
    *(Xlg2+g-1) = mltmod(Xa2w,*(Xlg2+g-1),Xm2);
    goto S50;
S40:
    printf("isdtyp not in range in INITGN");
    exit(1);
S50:
    *(Xcg1+g-1) = *(Xlg1+g-1);
    *(Xcg2+g-1) = *(Xlg2+g-1);
#undef numg
}


long ignlgi(void)
/*
**********************************************************************
     long ignlgi(void)
               GeNerate LarGe Integer
     Returns a random integer following a uniform distribution over
     (1, 2147483562) using the current generator.
     This is a transcription from Pascal to Fortran of routine
     Random from the paper
     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
     with Splitting Facilities." ACM Transactions on Mathematical
     Software, 17:98-111 (1991)
**********************************************************************
*/
{
#define numg 32L
extern void gsrgs(long getset,long *qvalue);
extern void gssst(long getset,long *qset);
extern void gscgn(long getset,long *g);
extern long Xm1,Xm2,Xa1,Xa2,Xcg1[],Xcg2[];
extern long Xqanti[];
static long ignlgi,curntg,k,s1,s2,z;
static long qqssd,qrgnin;
/*
     IF THE RANDOM NUMBER PACKAGE HAS NOT BEEN INITIALIZED YET, DO SO.
     IT CAN BE INITIALIZED IN ONE OF TWO WAYS : 1) THE FIRST CALL TO
     THIS ROUTINE  2) A CALL TO SETALL.
*/
    gsrgs(0L,&qrgnin);
    if(!qrgnin) inrgcm();
    gssst(0,&qqssd);
    if(!qqssd) setall(1234567890L,123456789L);
/*
     Get Current Generator
*/
    gscgn(0L,&curntg);
    s1 = *(Xcg1+curntg-1);
    s2 = *(Xcg2+curntg-1);
    k = s1/53668L;
    s1 = Xa1*(s1-k*53668L)-k*12211;
    if(s1 < 0) s1 += Xm1;
    k = s2/52774L;
    s2 = Xa2*(s2-k*52774L)-k*3791;
    if(s2 < 0) s2 += Xm2;
    *(Xcg1+curntg-1) = s1;
    *(Xcg2+curntg-1) = s2;
    z = s1-s2;
    if(z < 1) z += (Xm1-1);
    if(*(Xqanti+curntg-1)) z = Xm1-z;
    ignlgi = z;
    return ignlgi;
#undef numg
}

void inrgcm(void)
/*
**********************************************************************
     void inrgcm(void)
          INitialize Random number Generator CoMmon
                              Function
     Initializes common area  for random number  generator.  This saves
     the  nuisance  of  a  BLOCK DATA  routine  and the  difficulty  of
     assuring that the routine is loaded with the other routines.
**********************************************************************
*/
{
#define numg 32L
extern void gsrgs(long getset,long *qvalue);
extern long Xm1,Xm2,Xa1,Xa2,Xa1w,Xa2w,Xa1vw,Xa2vw;
extern long Xqanti[];
static long T1;
static long i;
/*
     V=20;                            W=30;
     A1W = MOD(A1**(2**W),M1)         A2W = MOD(A2**(2**W),M2)
     A1VW = MOD(A1**(2**(V+W)),M1)    A2VW = MOD(A2**(2**(V+W)),M2)
   If V or W is changed A1W, A2W, A1VW, and A2VW need to be recomputed.
    An efficient way to precompute a**(2*j) MOD m is to start with
    a and square it j times modulo m using the function MLTMOD.
*/
    Xm1 = 2147483563L;
    Xm2 = 2147483399L;
    Xa1 = 40014L;
    Xa2 = 40692L;
    Xa1w = 1033780774L;
    Xa2w = 1494757890L;
    Xa1vw = 2082007225L;
    Xa2vw = 784306273L;
    for(i=0; i<numg; i++) *(Xqanti+i) = 0;
    T1 = 1;
/*
     Tell the world that common has been initialized
*/
    gsrgs(1L,&T1);
#undef numg
}
