/* tdens.c      Last Modified 12/01/99  */
/* Evaluate t density.                  */
#include <math.h>

double tdens(double *x, double *mu, double *Cinv, int d, int df, double det)
{
  double c;
  double y;
  int inc=1;
  double *w;
  double pi=3.1415926535897931;
  void daxpy_(int *dim, double *alpha, double *xv, int *incx,
	      double *yv, int *incy);
  void dtrmv_(char *uplow, char *transf, char *udiag, int *dim, 
	      double *a , int *lda , double *xv,int *incr);
  double ddot_(int *nn, double *xv, int *incx, double *yv, 
	       int *incy);
  double dcopy_(int *nn, double *xv, int *incx, double *yv, 
	       int *incy);
  double a=-1.0;
  double gamma(double xa);
  char uplo[]="U";
  char trans[]="N";
  char diag[]="N";
  double *vecalloc(int nv);

  w=vecalloc(d);
  c=gamma(((double)(df+d))/2.0)/gamma(((double)(df))/2.0);
  c=c/pow(pi*((double)df),((double)d)/2.0);
  c=c*det;
  /* copy x to work vector w: w<-x */
  dcopy_(&d,x,&inc,w,&inc);
  /* The quadratic form: */
  /*  1) center w: w<-w-mu */ 
  daxpy_(&d,&a,mu,&inc,w,&inc);
  /*  2) standardize w: w<-C%*%w */
  dtrmv_(uplo,trans,diag,&d,&Cinv[0],&d,w,&inc);
  /*  3) form inner product */
  y=ddot_(&d,w,&inc,w,&inc); 
  y=c*pow(1 + (1/((double)df))*y,-1.0*(((double)(df+d))/2.0));
  return(y);
}  /* end of lfunc */