/* rant.c                              Last Modified 12/01/99  */
/* Generate pseudo-random deviate from multivariate t density. */
#include <math.h>

void rant(double *z, double *mu, double *Cmat, int d, int xdf)
{
  double y;
  int j, inc=1;
  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);
  float genchi(float dfs);
  float gennor(float mus, float sds);
  double a=1.0;
  char uplo[]="U";
  char trans[]="T";
  char diag[]="N";

  /* Chi-squared deviate, df degrees of freedom */
  y=((double)genchi(((float)xdf)));
  /* independent normal deviates */
  for (j=0;j<d;j++) z[j]=((double)gennor(0.0,1.0))*sqrt(((double)xdf)/y);
  /*  Scale z: z<-C%*%z */
  dtrmv_(uplo,trans,diag,&d,&Cmat[0],&d,z,&inc);
  /*  Center z: z<-z+mu */ 
  daxpy_(&d,&a,mu,&inc,z,&inc);
}  /* end of rant */