//This library of functions needs to be linked
//with arrays.cpp library functions.


//#define UNIX
#ifdef UNIX
#include <unistd.h>
#endif
#include <iostream.h>
#include <math.h>
#include "arrays.h"
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
double gammln(double xx);
double lfact(double n) ;
double Binomial(int n, int k) ;
double Multinomial(Vector n) ;
double random_num(int&);
void set_seed(void);
int seed;
/*
	double gammln(double x)
	returns the log of gamma(xx)
	most common use is because
	log(n!) = gammln(n+1)
*/
double gammln(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);
}
/*
	double lfact(double n)
	does table lookup to get log(n!) for
	small values of n, and calls
	gammln() for large values of n.
	It does table lookup for small factorials.
	Depending on the nature of your problem, you
	might want to change LFACT_TABLE_SIZE
	so you can always look up values.
	Even inline it if you use it a lot.
*/
#define LFACT_TABLE_SIZE 30
double lfact(double n)
{
  static double table[LFACT_TABLE_SIZE];
  static int ntop=0;
  static int test=LFACT_TABLE_SIZE-1;

  int j;

  if(n < 0.0) cerr << "ERROR: argument n= " << n << " to lfact() is negative\n";
  if(n > test) return(gammln(n+1.0));
  while(ntop < n)
	 {
		j=ntop++;
		table[ntop]=table[j]+log((double)ntop);
	 }
  return table[(int)n];
}

double Binomial(int n, int k)
{
  return (exp( lfact((double)n)-lfact((double)k)-lfact((double)n-k)));
}
double Multinomial(Vector n)
{
	double N=0.0,d=0.0;
	int i;
	for(i=0;i<n.length();i++) N+= n[i];
	for( i=0;i<n.length();i++) if(n[i] >0.0) d += lfact(n[i]);
	return exp(lfact(N) -d);
}
#ifdef UNIX

void set_seed(void)
{
	seed= -getpid();
	return;
}

double random_num(int& idum)
{
	static double y,v[100];
	static int RAND_MAX = (2 << 30)-1;
	static double maxran=((double)RAND_MAX +1.0);
	static int iff=0;
	int j;
	if (idum < 0 || iff == 0) {
		iff=1;
		srand(idum);
		idum=1;
		for (j=0;j<100;j++) v[j]=rand();
		y=rand();
	}
	j= (int)floor(100.0*y/maxran);
	if (j >=100 || j < 0) cerr << "ERROR random_num(int&): This cannot happen." << endl;
	y=v[j];
	v[j]=rand();
	return y/maxran;
}
#else
void set_seed(void)
{
	char timebuf[20];
	unsigned ta,tb,tc;
	_strtime(timebuf);
	sscanf(timebuf,"%d:%d:%d",&ta,&tb,&tc);
	seed=   -1*( (ta+ 20*(tb+30*tc))%1024   );
	if (seed > 0) seed=(-1*seed);
	return;
}
double random_num(int& idum)
{
	static double y,v[100];
	static double maxran=((double)RAND_MAX +1.0);
	static int iff=0;
	int j;
	if (idum < 0 || iff == 0) {
		iff=1;
		srand(idum);
		idum=1;
		for (j=0;j<100;j++) v[j]=rand();
		y=rand();
	}
	j= (int)floor(100.0*y/maxran);
	if (j >=100 || j < 0) cerr << "ERROR random_num(int&): This cannot happen." << endl;
	y=v[j];
	v[j]=rand();
	return y/maxran;
}
#endif