#include "neededroutines.h"

float gammln(float xx)
{
// This routine, from Numerical Recipes in C, calculates log of the gamma 
// function of x (remember that gamma(x) = (x-1)!, if x is an integer.
  double x, z, tmp, ser;
  double cof[6] = {76.18009173, -86.50532033, 24.01409822, -1.231739516, 
		   0.120858003e-2, -0.536382e-5};
  int i;

  if (xx < 1)
    {z = 1 - xx;
    return log(PI2) + log(z) - gammln(1+z) - log(sin(PI2 * z)); }
  x = xx - 1.0;
  tmp = x+5.5;
  tmp  -= (x+0.5) * log(tmp);
  ser = 1.0;
  for (i=0; i<=5; i++)
    {x += 1.0;
    ser += cof[i] / x; }
  return -tmp + log(2.50662827465*ser);
}

double randunif()
{
  JC = (JC * 69069) & 037777777777;  // congruential part 
  JT ^= JT >> 15;	  // tausworthe part
  JT ^= (JT << 17) & 037777777777;
  
  return(((JT ^ JC) >> 1) * Norm);
}

void sdrni(unsigned long* i)
{
  unsigned long k=*i;
//  if(k==0) k=time(0);
// I have had to comment out the line above and replace it with the line
// below because Visual C++ does not seem to recognise the time() command.
  k = 1;
  JT = k/65536; JC = k-65536*JT;
  JT = 65536*JT+1; JC = 32768*JC+1;
}


#define SORT2SWAPf(a,b) tempf=(a); (a)=(b); (b)=tempf;
#define SORT2SWAPi(a,b) tempi=(a); (a)=(b); (b)=tempi;
const int SORT2NSTACK=50;
const int SORT2M=7;

void sort2(unsigned long n, float arr[], int brr[])
{
  // Sorts array arr[1:n] into ascending order, while making the corresponding
  // rearrangement of the array brr[1:n]
  // (from p334 of Numerical Recipes in C)

  unsigned long i, ir=n, j, k, l=1;
  int *istack, jstack = 0;
  float a, tempf;
  int b, tempi;

  istack = new int [SORT2NSTACK+1];

  for (;;)
    {if (ir-l < SORT2M)
      {for (j=l+1; j<=ir; j++)
	{a = arr[j];
	b = brr[j];
	for (i=j-1; i>=1; i--)
	  {if (arr[i] <= a)
	    break;
	  arr[i+1] = arr[i];
	  brr[i+1] = brr[i];
	  }
	  arr[i+1] = a;
	  brr[i+1] = b;
	}
      if (!jstack)
	{delete [] istack;
	return;
	}
      ir = istack[jstack];
      l = istack[jstack-1];
      jstack -= 2;
      }
    else
      {k=(l+ir) >> 1; // This means (l+ir)/2 (>> is right shift operator)
      SORT2SWAPf(arr[k], arr[l+1]);
      SORT2SWAPi(brr[k], brr[l+1]);
      if (arr[l+1] > arr[ir])
	{SORT2SWAPf(arr[l+1], arr[ir]);
	SORT2SWAPi(brr[l+1], brr[ir]);
	}
      if (arr[l] > arr[ir])
	{SORT2SWAPf(arr[l], arr[ir]);
	SORT2SWAPi(brr[l], brr[ir]);
	}
      if (arr[l+1] > arr[l])
	{SORT2SWAPf(arr[l+1], arr[l]);
	SORT2SWAPi(brr[l+1], brr[l]);
	}
      i=l+1;
      j = ir;
      a = arr[l];
      b = brr[l];
      for (;;)
	{do {i++;} while (arr[i]<a);
	do {j--;} while (arr[j]>a);
	if (j<i) break;
	SORT2SWAPf(arr[i], arr[j]);
	SORT2SWAPi(brr[i], brr[j]);
	}
      arr[l] = arr[j];
      arr[j] = a;
      brr[l] = brr[j];
      brr[j] = b;
      jstack += 2;
      if (jstack > SORT2NSTACK)
	{cout << "ERROR in sort2: SORT2NSTACK too small";
	exit(0);
	}
      if (ir-i+1 >= j-l)
	{istack[jstack] = ir;
	istack[jstack-1] = i;
	ir = j-1;
	}
      else
	{istack[jstack] = j-1;
	istack[jstack-1] = l;
	l = i;
	}
      }
    }
}


void rank(unsigned long len, float arr[], int ranks[])
{
  // Ranks array arr[1:n] (based on sort2(), 
  // taken from p334 of Numerical Recipes in C)

  int i;
  float *arrcopy;
  int *index;

  arrcopy = new float [len+1];
  index = new int [len+1];
  for (i=1; i<=len; i++)
    {arrcopy[i] = (float) arr[i];
    index[i] = i;
    }

  sort2(len, arrcopy, index);
  for (i=1; i<=len; i++)
    ranks[index[i]] = i;
}


float rankties(unsigned long len, float arr[], float ranks[])
{
  // This is very similar to rank().  The difference is that this assigns
  // non-integer ranks when there are tied data.
  // It also returns a correction factor, CKW, for tied data.  This is used by
  // the Kruskal-Wallis test (see Biostatistics, Fisher and van Belle, p478).

  int i, begin, end, TL;
  float mid, CKW;
  int *intranks;

  intranks = new int [len+1];
  rank(len, arr, intranks);

  for (i=1; i<=len; i++)
    ranks[i] = (float) intranks[i];

  begin = 1;
  CKW = 0;

  while(begin<len)
    {end = begin;
    while (intranks[end+1]==intranks[begin] && end<len)
      end++;
    if (end>begin)
      {mid = (begin + end) * 1.0 / 2;
      for (i=begin; i<=end; i++)
	ranks[i] = mid;
      TL = end - begin + 1;
      CKW += pow(TL*1.0, 3) - TL;
      }
    begin = end+1;
    }

  CKW /= pow(len*1.0, 3) - len;

  return CKW;
}


float pchisq(int df, float x)
{
  // This returns the probability that X>x where X~chisq(df)

  return gammq(df*0.5, x/2);
}


float gammp(float a, float x)
{
  // Returns the incomplete gamma function P(a,x)
  // Function taken from Numerical Recipes in C, p218-9
  
  float gamser, gammcf, gln;
  
  if (x<0.0 || a<=0.0)
    {cout << "Invalid arguments in routine gammp\n";
    exit(0);
    }
  if (x < (a+1.0) )
    {// Use the series representation
      gser(&gamser, a, x, &gln);
      return gamser;
    }
  else
    {// Use the continued fraction representation
      gcf(&gammcf, a, x, &gln);
      return 1.0 - gammcf;
      // and take its complement
    }
}


float gammq(float a, float x)
{
  // Returns the incomplete gamma function Q(a,x) = 1 - P(a,x)
  // Function taken from Numerical Recipes in C, p218-9
  
  float gamser, gammcf, gln;
  
  if (x<0.0 || a<=0.0)
    {cout << "Invalid arguments in routine gammp\n";
    exit(0);
    }
  if (x < (a+1.0) )
    {// Use the series representation
      gser(&gamser, a, x, &gln);
      return 1.0 - gamser;
      // and take its complement
    }
  else
    {// Use the continued fraction representation
      gcf(&gammcf, a, x, &gln);
      return gammcf;
    }
}


void gser(float *gamser, float a, float x, float *gln)
{
  // Returns the incomplete gamma function P(a,x) evaluated by its series
  // representation as gamser.  Also returns ln Gamma(a) as gln.
  // Function taken from Numerical Recipes in C, p218-9
  
  int n;
  float sum, del, ap;
  
  *gln = gammln(a);
  if (x<=0.0)
    {if (x<0.0)
      {cout << "x less than 0 in routine gser\n";
      exit(0);
      }
    *gamser = 0.0;
    return;
    }
  else
    {ap = a;
    del = sum = 1.0/a;
    for (n=1; n<=ITMAX; n++)
      {++ap;
      del *= x/ap;
      sum += del;
      if (fabs(del) < fabs(sum)*EPS)
	{*gamser = sum*exp(-x+a*log(x)-(*gln));
	return;
	}
      }
    cout << "a too large, or ITMAX (set to 100) too small in routine gser";
    exit(1);
    return;
    }
}
      

void gcf(float *gammcf, float a, float x, float *gln)
{
  // Returns the incomplete gamma function Q(a,x) evaluated by its continued
  // fraction representation as gammcf.  Also returns ln Gamma(a) as gln.
  // Function taken from Numerical Recipes in C, p218-9
 
  int i;
  float an, b, c, d, del, h;

  *gln = gammln(a);
  b = x+1.0-a;
  c = 1.09/FPMIN;
  d=1.0/b;
  h=d;

  for (i=1; i<=ITMAX; i++)
    {an = -i*(i-a);
    b += 2.0;
    d = an*d+b;
    if (fabs(d) < FPMIN)
      d = FPMIN;
    c = b+an/c;
    if (fabs(c) < FPMIN)
      c = FPMIN;
    d = 1.0/d;
    del = d*c;
    h *= del;
    if (fabs(del-1.0) < EPS)
      break;
    }
 
  if (i>ITMAX)
    {cout << "a too large, ITMAX too small in gcf";
    exit(0);
    }
  *gammcf = exp(-x+a*log(x)-(*gln))*h;
}