/* The internal random generators need to be converted to double precision... */#include "math.h"#include "matrix.h"#include "mex.h"#include "string.h"#include <stdlib.h>#include <stdio.h>#include "randlib.h"#define USAGE "\n P = RT_rand(pdftype,seed1,seed2,param1,param2...,{mrows,ncols})\n A mex interface for the RANDLIB libraries by Glen Davidson.\n P is a random variate from distribution PDFTYPE with parameters param* \n THE SEEDs SHOULD BE PASSED AS rand FROM Matlab \n see randlib.c.fdoc for usage. \n PDFTYPE = [1..13] as beta(1), binomial(2), chisquare(3), \n non-central chisquare(4), F(5), non-central F(6), gamma(7), negative binomial(8), normal(9), \n poisson(10), not_implemented(11), not_implemented(12), exponential(13). \n If the parameters are not all scalar, the common size will be returned, else use MROWS and NCOLS.\n use RT_rand(pdftype) to display required parameters.\n\n NOTE THE VARIATES ARE SINGLE PRECISION (CONVERTED FROM A C float)"#define BETAdist 1#define BINOMIALdist 2#define CHISQUAREdist 3/* non-central chisquare */#define NCCHISQUAREdist 4/* non-central F distribution */#define Fdist 5#define NCFdist 6#define GAMMAdist 7/* negative binomial */#define NEGBINOMIALdist 8#define NORMALdist 9#define POISSONdist 10/* Student's T distribution NOT IMPLEMENTED*/#define STUDENTTdist 11/* Non-central T distribution NOT IMPLEMENTED*/#define NCTdist 12#define EXPdist 13void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ){  double p;			/* outputs from core routines */  double *p_outputPtr = NULL; /* pointer to outputs */  int nel, mrows, ncols, pdftype = 0;  int do_mrows_ncols = -1;  int mrows_input = -1, ncols_input = -1; /* matrix output from scalar parameters */  int nel_check_mrows, nel_check_ncols; /* matrix output from scalar parameters */  int nel_check, mrows_check, ncols_check;  int rloop, tloop;  long seed1_long, seed2_long;	/* RNG needs LONG seed initialiser */  double seed1_double, seed2_double; /* rand passed from Matlab */  int NPDF=13;			/* total of 10 pdfs in the package */  /* following is number of params for each distribution, but indexed from 1 so NparamsPtr[0] invalid */  int NparamsPtr[] = {-1, 2, 2, 1, 2, 2, 3, 2, 2, 2, 1, 0, 0, 1};  /* following is number of elements within each parameter to allow matrix inputs, 10 allocated for future use */  int NelementsPtr[10];  /* pointer to the actual elements */  double *elementsPtrPtr[10];  /* FOLLOWING USES FLOAT FOR THE SINGLE-PRECISION RANDOM NUMBER LIBRARY */  /* array space to pass to a single routine */  float param_passPtr[10];  long binparam;		/* binomial & negative binomial require a long input */   if (nrhs < 1) {		/* minimum number of 1 inputs */    mexPrintf("\n Incorrect number of input arguments");    mexErrMsgTxt(USAGE);  }  if (nlhs > 1) {    mexPrintf("\n Incorrect number of output arguments");    mexErrMsgTxt(USAGE);  }  nel_check = mxGetNumberOfElements(prhs[0]);  if (nel_check == 1) {    pdftype = (int)mxGetScalar(prhs[0]);  } else {    mexPrintf("\n PDFTYPE should be a scalar for the particular cdf");    mexErrMsgTxt(USAGE);  }  if ( (pdftype < 1) || (pdftype > NPDF)) {    mexPrintf("\n PDFTYPE should range from 1 to %d",NPDF);    mexErrMsgTxt(USAGE);  }  /* check for correct number of inputs (subtracting 1 for the pdftype input and 2 for the seeds)*/  if ( (nrhs-3) == NparamsPtr[pdftype] ) {/* no mrows or ncols */    do_mrows_ncols = 0;  } else if ( (nrhs-5) == NparamsPtr[pdftype] ) { /* mrows and ncols input */    do_mrows_ncols = 1;    nel_check_mrows = mxGetNumberOfElements(prhs[nrhs-2]);    nel_check_ncols = mxGetNumberOfElements(prhs[nrhs-1]);    if ( (nel_check_mrows != 1) || (nel_check_ncols != 1) ) {      mexPrintf("\n Need scalar inputs for MROWS and NCOLS");      mexErrMsgTxt(USAGE);    }    mrows_input = (int)mxGetScalar(prhs[nrhs-2]);    ncols_input = (int)mxGetScalar(prhs[nrhs-1]);    if ( (mrows_input < 1) || (ncols_input < 1) ) {      mexPrintf("\n Need scalar inputs > 0 for MROWS and NCOLS");      mexErrMsgTxt(USAGE);    }  } else {    switch(pdftype) {    case BETAdist: mexPrintf("\n Need A and B for BETA distribution");break;    case BINOMIALdist: mexPrintf("\n Need XN, and PR for BINOMIAL distribution");break;    case CHISQUAREdist: mexPrintf("\n Need DF for CHISQUARE distribution");break;    case NCCHISQUAREdist: mexPrintf("\n Need DF and PNONC for NON CENTRAL CHISQUARE distribution");break;    case Fdist: mexPrintf("\n Need DFN and DFD for F distribution");break;    case NCFdist: mexPrintf("\n Need DFN, DFD and PNONC for NON CENTRAL F distribution");break;    case GAMMAdist: mexPrintf("\n Need SHAPE and SCALE for GAMMA distribution");break;    case NEGBINOMIALdist: mexPrintf("\n Need XN and PR for NEGATIVE BINOMIAL distribution");break;    case NORMALdist: mexPrintf("\n Need MEAN, and SD for NORMAL distribution");break;    case POISSONdist: mexPrintf("\n Need XLAM for POISSON distribution");break;    case EXPdist: mexPrintf("\n Need MEAN for EXPONENTIAL distribution");break;    default:      mexPrintf("\n Case not implemented");      mexErrMsgTxt(USAGE);    }    mexPrintf("\n with optional MROW and NCOL input");    mexErrMsgTxt(USAGE);  }  nel_check = mxGetNumberOfElements(prhs[1]); /* seed1 */  if (nel_check == 1) {    seed1_double = mxGetScalar(prhs[1]);    if ( (seed1_double <= 0.0) | (seed1_double >= 1.0) ) {      mexPrintf("\n SEED1 should be a rand output from 0.0 < SEED1 < 1.0");      mexErrMsgTxt(USAGE);    }  } else {    mexPrintf("\n SEED1 should be a scalar 0.0 < SEED1 < 1.0 ");    mexErrMsgTxt(USAGE);    seed1_double = -1;  }  nel_check = mxGetNumberOfElements(prhs[2]); /* seed2 */  if (nel_check == 1) {    seed2_double = mxGetScalar(prhs[2]);    if ( (seed2_double <= 0.0) | (seed2_double >= 1.0) ) {      mexPrintf("\n SEED2 should be a rand output from 0.0 < SEED2 < 1.0");      mexErrMsgTxt(USAGE);    }  } else {    mexPrintf("\n SEED2 should be a scalar 0.0 < SEED2 < 1.0 ");    mexErrMsgTxt(USAGE);    seed2_double = -1;  }    /* very laborious test to determine if input sizes are valid,     this means they should be scalars or the size of mrows_out, ncols_out */  nel = 1;		/* assumed number of output elements */  mrows = 1;		/* assumed rows of output elements */  ncols = 1;		/* assumed columns of output elements */    for (tloop = 1; tloop <= NparamsPtr[pdftype]; tloop++) { /* [pdftype,seed1,seed2,1,2,3..N,{mrows,ncols}] */    nel_check = mxGetNumberOfElements(prhs[tloop+2]); /* tloop+2 skips pdftype, seed1 and seed2 */    if (nel_check == 0) {      mexPrintf("\n Can't process parameters passed as []");      mexErrMsgTxt(USAGE);    }    NelementsPtr[tloop] = nel_check; /* store number of elements */    elementsPtrPtr[tloop] = mxGetPr(prhs[tloop+2]); /* and actual element pointer reference */    if (nel_check > 1) {		/* is this input a matrix */      if (nel > 1) {		/* is there already a matrix input */	mrows_check = mxGetM(prhs[tloop+2]); /* tloop+2 skips pdftype, seed1 and seed2 */	ncols_check = mxGetN(prhs[tloop+2]); /* tloop+2 skips pdftype, seed1 and seed2 */	if ( (mrows_check != mrows) || (ncols_check != ncols) ) {	  mexPrintf("\n Input sizes are invalid");	  mexErrMsgTxt(USAGE);	}      } else {			/* this is the first matrix input */	nel = mxGetNumberOfElements(prhs[tloop+2]); /* tloop+2 skips pdftype, seed1 and seed2 */	mrows = mxGetM(prhs[tloop+2]); /* tloop+2 skips pdftype, seed1 and seed2 */	ncols = mxGetN(prhs[tloop+2]); /* tloop+2 skips pdftype, seed1 and seed2 */      }    }  }  /* Only if all the parameters are scalar and an mrows_input and ncols_input is present, can     a matrix of similarly distributed variates be output */  if (do_mrows_ncols == 1) {	/* mrows_input and ncols_input */    if ( (mrows == 1) && (ncols == 1) ) { /* scalar parameters */      mrows = mrows_input;      ncols = ncols_input;      nel = mrows_input * ncols_input;    } else {      mexPrintf("\n Parameters must be scalar for a matric MROWS x NCOLS output");      mexErrMsgTxt(USAGE);    }  }      				/* create output (p) */  plhs[0] = mxCreateDoubleMatrix(mrows, ncols, mxREAL);  p_outputPtr = mxGetPr(plhs[0]);  /* Note problem with all this is that it's single precision... param_passPtr modified for this */  /* Initialise random number generator: From Basegen.c.doc in randlib*/  /* The state of a generator  is determined by two  integers called seeds.     The seeds can be  initialized  by the  user; the initial values of the     first must lie between 1 and 2,147,483,562, that of the second between     1 and 2,147,483,398.  Each time a number is generated,  the  values of     the seeds changes.  */  seed1_long = (long)(seed1_double * 1073741824 + 1);	/* a seed of 1..2^30 is large enough */  seed2_long = (long)(seed2_double * 1073741824 + 1);	/* a seed of 1..2^30 is large enough */  setall( seed1_long, seed2_long );  for (rloop=0; rloop < nel; rloop++) {    for (tloop=1; tloop <= NparamsPtr[pdftype]; tloop++) {      if (NelementsPtr[tloop] == 1) { /* if only one element input, always pass this */	param_passPtr[tloop] = (float)elementsPtrPtr[tloop][0];      } else {			/* else pass the array over the loop */	param_passPtr[tloop] = (float)elementsPtrPtr[tloop][rloop];      }    } /* finished forming inputs */    switch (pdftype) {    case BETAdist:      p = (double)genbet(param_passPtr[1], param_passPtr[2]); break;    case BINOMIALdist:      binparam = (long)param_passPtr[1]; /* need to pass a long integer */      p = (double)ignbin(binparam, param_passPtr[2]); break;    case CHISQUAREdist:      p = (double)genchi(param_passPtr[1]); break;    case NCCHISQUAREdist:      p = (double)gennch(param_passPtr[1], param_passPtr[2]); break;    case Fdist:      p = (double)genf(param_passPtr[1], param_passPtr[2]); break;    case NCFdist:      p = (double)gennf(param_passPtr[1], param_passPtr[2], param_passPtr[3]); break;    case GAMMAdist:       p = (double)gengam(param_passPtr[1], param_passPtr[2]); break;    case NEGBINOMIALdist:       binparam = (long)param_passPtr[1]; /* need to pass a long integer */      p = (double)ignnbn(binparam, param_passPtr[2]); break;    case NORMALdist:       p = (double)gennor(param_passPtr[1], param_passPtr[2]); break;    case POISSONdist:       p = (double)ignpoi(param_passPtr[1]); break;    case EXPdist:       p = (double)genexp(param_passPtr[1]); break;    default:      mexErrMsgTxt("\n Case not implemented");      p = 0;    }    p_outputPtr[rloop] = p;  }}