/* The ziggurat method for RNOR and REXP
Combine the code below with the main program in which you want
normal or exponential variates.   Then use of RNOR in any expression
will provide a standard normal variate with mean zero, variance 1,
while use of REXP in any expression will provide an exponential variate
with density exp(-x),x>0.
Before using RNOR or REXP in your main, insert a command such as
zigset(86947731 );
with your own choice of seed value>0, rather than 86947731.
(If you do not invoke zigset(...) you will get all zeros for RNOR and REXP.)
For details of the method, see Marsaglia and Tsang, "The ziggurat method
for generating random variables", Journ. Statistical Software.
*/

#include <ctime>
#include <climits>
#include <iostream>
#include <octave/oct.h>
#include <octave/lo-mappers.h>

static bool initialized = false;
static unsigned long jz,jsr=123456789;

#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
#define UNI (.5 + (signed) SHR3*.2328306e-9)
#define IUNI SHR3

static long hz;
static unsigned long iz, kn[128], ke[256];
static double wn[128],fn[128], we[256],fe[256];

#define RNOR (hz=SHR3, iz=hz&127, (abs(hz)<kn[iz])? hz*wn[iz] : nfix())
#define REXP (jz=SHR3, iz=jz&255, (    jz <ke[iz])? jz*we[iz] : efix())

/* nfix() generates variates from the residue when rejection in RNOR occurs. */

double nfix(void)
{
const double r = 3.442620;     /* The start of the right tail */
static double x, y;
 for(;;)
  {  x=hz*wn[iz];      /* iz==0, handles the base strip */
     if(iz==0)
       { do{ x=-log(UNI)*0.2904764; y=-log(UNI);}	/* .2904764 is 1/r */
        while(y+y<x*x);
        return (hz>0)? r+x : -r-x;
       }
                         /* iz>0, handle the wedges of other strips */
      if( fn[iz]+UNI*(fn[iz-1]-fn[iz]) < exp(-.5*x*x) ) return x;

     /* initiate, try to exit for(;;) for loop*/
      hz=SHR3;
      iz=hz&127;
      if(fabs(hz)<kn[iz]) return (hz*wn[iz]);
  }

}

/* efix() generates variates from the residue when rejection in REXP occurs. */
double efix(void)
{ 
  double x;
  for(;;)
    {  
      if(iz==0) 
	return (7.69711-log(UNI));          /* iz==0 */
      x=jz*we[iz]; 
      if( fe[iz]+UNI*(fe[iz-1]-fe[iz]) < exp(-x) ) 
	return (x);

      /* initiate, try to exit for(;;) loop */
      jz=SHR3;
      iz=(jz&255);
      if(jz<ke[iz]) return (jz*we[iz]);
    }
}
/*--------This procedure sets the seed and creates the tables------*/

void zigset(unsigned long jsrseed)
{  
  const double m1 = 2147483648.0, m2 = 4294967296.;
  double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
  double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
  int i;

  // Want seed to be easily set
  // jsr^=jsrseed;
  jsr = jsrseed;

  /* Set up tables for RNOR */
  q=vn/exp(-.5*dn*dn);
  kn[0]=(unsigned long int)((dn/q)*m1);
  kn[1]=0;

  wn[0]=q/m1;
  wn[127]=dn/m1;

  fn[0]=1.;
  fn[127]=exp(-.5*dn*dn);

  for(i=126;i>=1;i--)
    {
      dn=sqrt(-2.*log(vn/dn+exp(-.5*dn*dn)));
      kn[i+1]=(unsigned long int)((dn/tn)*m1);
      tn=dn;
      fn[i]=exp(-.5*dn*dn);
      wn[i]=dn/m1;
    }

  /* Set up tables for REXP */
  q = ve/exp(-de);
  ke[0]=(unsigned long int)((de/q)*m2);
  ke[1]=0;

  we[0]=q/m2;
  we[255]=de/m2;

  fe[0]=1.;
  fe[255]=exp(-de);

  for(i=254;i>=1;i--)
    {
      de=-log(ve/de+exp(-de));
      ke[i+1]= (unsigned long int)((de/te)*m2);
      te=de;
      fe[i]=exp(-de);
      we[i]=de/m2;
    }
}

// Octave interface starts here
static inline unsigned long hash( time_t t, clock_t c )
{
  // Get a unsigned long from t and c
  // Better than (unsigned long)(x) in case x is floating point in [0,1]
  // Based on code by Lawrence Kirby (address@hidden)

  static unsigned long differ = 0;  // guarantee time-based seeds will change

  unsigned long h1 = 0;
  unsigned char *p = (unsigned char *) &t;
  for( size_t i = 0; i < sizeof(t); ++i )
    {
      h1 *= UCHAR_MAX + 2U;
      h1 += p[i];
    }
  unsigned long h2 = 0;
  p = (unsigned char *) &c;
  for( size_t j = 0; j < sizeof(c); ++j )
    {
      h2 *= UCHAR_MAX + 2U;
      h2 += p[j];
    }
  return ( h1 + differ++ ) ^ h2;
}

static unsigned long
get_seed (void)
{
  FILE* urandom = fopen( "/dev/urandom", "rb" );
  bool success = false;
  unsigned long s;
  if( urandom )
    {
      success = fread (&s, sizeof(unsigned long), 1, urandom); 
      fclose(urandom);
    }
  if ( !success )
    // Was not successful, so use time() and clock() instead
    s = (unsigned long int) hash( time(NULL), clock() );

  // There a chance that urandom gave zero -> make it 1
  if (s == 0)
    s = 1;

  return s;
}

static octave_value 
do_seed (octave_value_list args)
{
  octave_value retval;

  // Check if they said the magic words
  std::string s_arg = args(0).string_value ();
  if (s_arg == "seed")
    retval = (double)jsr;
  else if (s_arg == "state")
    // Don't have state vector, return seed here as well
    retval = (double)jsr;
  else
    {
      error ("rand: unrecognized string argument");
      return retval;
    }

  // Check if just getting seed
  if (args.length() == 1)
    return retval;

  // Set the seed from a scalar
  octave_value tmp = args(1);
  if (tmp.is_scalar_type ()) 
    {
      unsigned long int seed = (unsigned long)tmp.double_value();
      if (!initialized)
	{
	  zigset ((seed == 0 ? 1 : seed));
	  initialized = true;
	}
      else
	jsr = seed;
    }
  else
    error ("randn: not a valid seed");
  
  return retval;
}

static void
do_size(octave_value_list args, int& nr, int& nc)
{
  int nargin = args.length();

  if (nargin == 0)
    {
      nr = nc = 1;
    }
  else if (nargin == 1)
    {
      octave_value tmp = args(0);

      if (tmp.is_scalar_type ())
	{
	  double dval = tmp.double_value ();
	  
	  if (xisnan (dval))
	    {
	      error ("randn: NaN is invalid a matrix dimension");
	    }
	  else
	    {
	      nr = nc = NINT (tmp.double_value ());
	    }
	}
      else if (tmp.is_range ())
	{
	  Range rng = tmp.range_value ();
	  nr = 1;
	  nc = rng.nelem ();
	}
      else if (tmp.is_matrix_type ())
	{
	  // XXX FIXME XXX -- this should probably use the function
	  // from data.cc.

	  Matrix a = args(0).matrix_value ();

	  if (error_state)
	    return;
	  
	  nr = a.rows ();
	  nc = a.columns ();
	  
	  if (nr == 1 && nc == 2)
	    {
	      nr = NINT (a (0, 0));
	      nc = NINT (a (0, 1));
	    }
	  else if (nr == 2 && nc == 1)
	    {
	      nr = NINT (a (0, 0));
	      nc = NINT (a (1, 0));
	    }
	  else
	    warning ("randn (A): use randn (size (A)) instead");
	}
      else
	{
	  gripe_wrong_type_arg ("randn", tmp);
	}
    }
  else if (nargin == 2)
    {
      double rval = args(0).double_value ();
      double cval = args(1).double_value ();
      if (! error_state)
	{
	  if (xisnan (rval) || xisnan (cval))
	    {
	      error ("randn: NaN is invalid as a matrix dimension");
	    }
	  else
	    {
	      nr = NINT (rval);
	      nc = NINT (cval);
	    }
	}
    }
}

DEFUN_DLD (zigurat, args, nargout, 
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {} randn (@var{x})\n\
@deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\
@deftypefnx {Loadable Function} address@hidden =} randn (\"state\", @var{x})\n\
@deftypefnx {Loadable Function} address@hidden =} randn (\"seed\", @var{x})\n\
Return a matrix with normally distributed random elements.  The\n\
arguments are handled the same as the arguments for @code{rand}.\n\
\n\
@code{randn} uses Ahrens and Dieter (1973) to transform from U to N(0,1).\n\n\
@end deftypefn\n\
@seealso{rand}\n")
{
  octave_value_list retval;	// list of return values

  int nargin = args.length ();	// number of arguments supplied
  if (nargin > 2) 
    print_usage("randn");

  else if (nargin > 0 && args(0).is_string())
    retval(0) = do_seed (args);

  else
    {
      if (!initialized)
	{
	  zigset(get_seed());
	  initialized = true;
	}

      int nr=0, nc=0;
      do_size (args, nr, nc);

      if (! error_state)
	{
	  Matrix X(nr, nc);

	  for (int c=0; c < nc; c++)
	    for (int r=0; r < nr; r++)
	      X(r,c) = RNOR ;
	  retval(0) = X;
	}
    }

  return retval;
}

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/