/*
  Copyright 2005 David Bateman

  I grant this code to the public domain
*/

#include <string.h>
#include <math.h>
#include "mex.h"
#include <stdint.h>

typedef mwSize octave_idx_type;

#include "randmtzig.c"

static int init_generator = 0;

static mxArray *
do_seed (int nargin, const mxArray *pargin [ ], const char *nm)
{

  char *s_arg;
  int buflen, i;
  double *out;
  mxArray *ret, *state;

  buflen = (mxGetM(pargin[0]) * mxGetN(pargin[0])) + 1;
  
  s_arg = (char *) mxCalloc (buflen, sizeof(char));

  if (mxGetString (pargin[0], s_arg, buflen))
    mexWarnMsgTxt ("Not enough space. String is truncated.");
    
  if (strcmp (s_arg, "seed") == 0)
    {
      uint32_t a = randi32();
      oct_init_by_int(a);
      init_generator = 1;

      ret = mxCreateDoubleMatrix (1, 1, mxREAL);
      out = mxGetPr (ret);
      out[0] = (double) a;
    }
  else if (strcmp (s_arg, "state") == 0)
    {
      if (!init_generator)
	{
	  init_generator = 1;
	  oct_init_by_entropy ();
	}

      uint32_t state[MT_N+1];
      oct_get_state(state);

      ret = mxCreateDoubleMatrix (MT_N+1, 1, mxREAL);
      out = mxGetPr (ret);
      for (i=0; i < MT_N+1; i++)
	out[i] = (double) (state[i]);

    }
  else
    {
      mexErrMsgIdAndTxt ("MATLAB:rand:UnrecognizedArg", "%s: unrecognized string argument", nm);
      return ret;
    }

  if (nargin == 1)
    return ret;

  /* Set the state from either a scalar or a previously returned state vector */
  state = (mxArray *)pargin[1];
  if (mxGetNumberOfElements (state) == 1)
    {
      uint32_t n = (uint32_t)((mxGetPr (state))[0]);
      oct_init_by_int(n);
    }
  else if (mxGetN (state) == 1 || mxGetM (state) == 1)
    {
      double *tmp;
      int input_len, n, i;
      uint32_t istate[MT_N+1];

      tmp = mxGetPr (state);
      input_len = mxGetNumberOfElements (state);
      n = input_len < MT_N+1 ? input_len : MT_N+1;

      for (i = 0; i < n; i++) 
	istate[i] = (uint32_t)tmp[i];

      if (input_len == MT_N+1 && istate[MT_N] <= MT_N && istate[MT_N] > 0)
	oct_set_state (istate);
      else
	oct_init_by_array (istate, n);
    }
  else
    error ("%s: not a state vector", nm);

  return ret;
}

void mexFunction
(    int nargout,		/* number of outputs */
    mxArray *pargout [ ],	/* output arguments */
    int nargin,			/* number of inputs */
    const mxArray *pargin [ ]	/* input arguments */
)
{
  int i, sz, *ndim, n;
  double *out;
  const char *nm;
  

  nm = mexFunctionName ();

  if (nargin >= 1 && mxIsChar (pargin[0]))
    pargout [0] = do_seed (nargin, pargin,nm);
  else
    {
      
      init_generator = 1;
      if (nargin == 0)
	{
	  ndim = (int *) mxMalloc (sizeof (int));
	  ndim[0] = 1;
	  sz = 1;
	  n = 1;
	}
      else 
	{
	  if (nargin == 1)
	    {
	      double *d;
	      ndim = (int *) mxMalloc (2 * sizeof (int));
	      d = mxGetPr (pargin [0]);
	      ndim[0] = ndim[1] = (int) d[0];
	      sz = ndim[1] * ndim[0];
	      n = 2;
	    }
	  else
	    {
	      n = nargin;
	      ndim = (int *) mxMalloc (nargin * sizeof (int));
	      sz = 1;
  
	      for (i = 0; i < nargin; i++)
		{
		  double *d;
		  d = mxGetPr (pargin [i]);
		  ndim[i] = (int) d[0];
		  sz *= ndim[i];
		}
	    }
	}

      pargout [0] = mxCreateNumericArray (n, ndim, mxDOUBLE_CLASS, mxREAL);
      out = mxGetPr (pargout [0]);
  
      if (strcmp (nm, "randn") == 0)
	oct_fill_randn (sz, out);
      else
	oct_fill_randu (sz, out);
    }
  
  return ;
}