/*
 * Copyright (c) 2000-2003 Stephen Williams (steve@icarus.com)
 *
 *    This source code is free software; you can redistribute it
 *    and/or modify it in source code form under the terms of the GNU
 *    General Public License as published by the Free Software
 *    Foundation; either version 2 of the License, or (at your option)
 *    any later version.
 *
 *    This program is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with this program; if not, write to the Free Software
 *    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
 */
#ifdef HAVE_CVS_IDENT
#ident "$Id: sys_random.c,v 1.11 2004/06/09 22:14:10 steve Exp $"
#endif

# include "sys_priv.h"

# include  <vpi_user.h>
# include  <assert.h>
# include  <stdlib.h>
# include  <math.h>
# include  <limits.h>

#if ULONG_MAX > 4294967295UL
# define UNIFORM_MAX INT_MAX
# define UNIFORM_MIN INT_MIN
#else
# define UNIFORM_MAX LONG_MAX
# define UNIFORM_MIN LONG_MIN
#endif

static double uniform(long *seed, long start, long end);
static long poisson(long *seed, long mean);

long rtl_dist_poisson(long*seed, long mean)
{
      long i;

      if (mean > 0) {
	    i = poisson(seed,mean);

      } else {
	    vpi_printf("WARNING: Poisson distribution must have "
		       "a positive mean\n");
	    i = 0;
      }

      return 0;
}

/* copied from IEEE1364-2001, with slight midifications for 64bit machines. */
long rtl_dist_uniform(long*seed, long start, long end)
{
	double r;
	long i;

	if (start >= end) return(start);

	if (end != UNIFORM_MAX)
	{
		end++;
		r = uniform( seed, start, end );
		if (r >= 0)
		{
			i = (long) r;
		}
		else
		{
			i = (long) (r-1);
		}
		if (i<start) i = start;
		if (i>=end) i = end-1;
	}
	else if (start!=UNIFORM_MIN)
	{
		start--;
		r = uniform( seed, start, end) + 1.0;
		if (r>=0)
		{
			i = (long) r;
		}
		else
		{
			i = (long) (r-1);
		}
		if (i<=start) i = start+1;
		if (i>end) i = end;
	}
	else
	{
		r = (uniform(seed,start,end)+2147483648.0)/4294967295.0;
		r = r*4294967296.0-2147483648.0;
		if (r>=0)
		{
			i = (long) r;
		}
		else
		{
			i = (long) (r-1);
		}
	}

	return (i);
}

static double uniform(long *seed, long start, long end )
{
	double d = 0.00000011920928955078125;
	double a, b, c;
	unsigned long oldseed, newseed;

	oldseed = *seed;
	if (oldseed == 0)
	      oldseed = 259341593;

	if (start >= end) {
	      a = 0.0;
	      b = 2147483647.0;
	} else {
	      a = (double)start;
	      b = (double)end;
	}

	/* Original routine used signed arithmetic, and the (frequent)
	 * overflows trigger "Undefined Behavior" according to the
	 * C standard (both c89 and c99).  Using unsigned arithmetic
	 * forces a conforming C implementation to get the result
	 * that the IEEE-1364-2001 committee wants.
	 */
	newseed = 69069 * oldseed + 1;

	/* Emulate a 32-bit unsigned long, even if the native machine
	 * uses wider words.
	 */
#if ULONG_MAX > 4294967295UL
	newseed = newseed & 4294967295UL;
#endif
	*seed = newseed;


#if 0
	/* Cadence-donated conversion from unsigned int to double */
	{
		union { float s; unsigned stemp; } u;
		u.stemp = (newseed >> 9) | 0x3f800000;
		c = (double) u.s;
	}
#else
	/* Equivalent conversion without assuming IEEE 32-bit float */
	/* constant is 2^(-23) */
	c = 1.0 + (newseed >> 9) * 0.00000011920928955078125;
#endif


	c = c + (c*d);
	c = ((b - a) * c - 1.0) + a;

	return c;
}

/* copied from IEEE1364-2001, with slight midifications for 64bit machines. */
static long poisson(long*seed, long mean)
{
      long n;
      double p, q;

      n = 0;
      q = -(double)mean;
      p = exp(q);
      q = uniform(seed, 0, 1);
      while (p < q) {
	    n++;
	    q = uniform(seed,0,1) * q;
      }

      return n;
}

static int sys_dist_poisson_calltf(char*name)
{
      s_vpi_value val;
      vpiHandle call_handle;
      vpiHandle argv;
      vpiHandle seed, mean;

      long i_seed, i_mean;

      call_handle = vpi_handle(vpiSysTfCall, 0);
      assert(call_handle);

	/* The presence of correct parameters should be checked at
	   compile time by the compiletf function. */
      argv = vpi_iterate(vpiArgument, call_handle);
      assert(argv);
      seed = vpi_scan(argv);
      assert(seed);
      mean = vpi_scan(argv);
      assert(mean);
      vpi_free_object(argv);

      val.format = vpiIntVal;
      vpi_get_value(seed, &val);
      i_seed = val.value.integer;

      vpi_get_value(mean, &val);
      i_mean = val.value.integer;

      val.format = vpiIntVal;
      val.value.integer = rtl_dist_poisson(&i_seed, i_mean);
      vpi_put_value(call_handle, &val, 0, vpiNoDelay);

      val.format = vpiIntVal;
      val.value.integer = i_seed;
      vpi_put_value(seed, &val, 0, vpiNoDelay);

      return 0;
}

static int sys_dist_poisson_sizetf(char*x)
{
      return 32;
}

static int sys_dist_uniform_calltf(char*name)
{
      s_vpi_value val;
      vpiHandle call_handle;
      vpiHandle argv;
      vpiHandle seed = 0, start, end;

      long i_seed, i_start, i_end;

      call_handle = vpi_handle(vpiSysTfCall, 0);
      assert(call_handle);

      argv = vpi_iterate(vpiArgument, call_handle);
      if (argv == 0) {
	    vpi_printf("ERROR: %s requires parameters "
		       "(seed, start, end)\n", name);
	    return 0;
      }

      seed = vpi_scan(argv);
      assert(seed);
      start = vpi_scan(argv);
      assert(start);
      end = vpi_scan(argv);
      assert(end);

      vpi_free_object(argv);

      val.format = vpiIntVal;
      vpi_get_value(seed, &val);
      i_seed = val.value.integer;

      vpi_get_value(start, &val);
      i_start = val.value.integer;

      vpi_get_value(end, &val);
      i_end = val.value.integer;

      val.format = vpiIntVal;
      val.value.integer = rtl_dist_uniform(&i_seed, i_start, i_end);
      vpi_put_value(call_handle, &val, 0, vpiNoDelay);

      val.format = vpiIntVal;
      val.value.integer = i_seed;
      vpi_put_value(seed, &val, 0, vpiNoDelay);

      return 0;
}

static int sys_dist_uniform_sizetf(char*x)
{
      return 32;
}

static int sys_random_calltf(char*name)
{
      s_vpi_value val;
      vpiHandle call_handle;
      vpiHandle argv;
      vpiHandle seed = 0;
      long i_seed = 0;

      call_handle = vpi_handle(vpiSysTfCall, 0);
      assert(call_handle);

	/* Get the argument list and look for a seed. If it is there,
	   get the value and reseed the random number generator. */
      argv = vpi_iterate(vpiArgument, call_handle);
      if (argv) {
	    seed = vpi_scan(argv);
	    vpi_free_object(argv);

	    val.format = vpiIntVal;
	    vpi_get_value(seed, &val);
	    i_seed = val.value.integer;
      }

      val.format = vpiIntVal;
      val.value.integer = rtl_dist_uniform(&i_seed, INT_MIN, INT_MAX);

      vpi_put_value(call_handle, &val, 0, vpiNoDelay);

        /* Send updated seed back to seed parameter. */
      if (seed) {
	    val.format = vpiIntVal;
	    val.value.integer = i_seed;		
	    vpi_put_value(seed, &val, 0, vpiNoDelay);
      }

      return 0;
}

static int sys_random_sizetf(char*x)
{
      return 32;
}

void sys_random_register()
{
      s_vpi_systf_data tf_data;

      tf_data.type   = vpiSysFunc;
      tf_data.tfname = "$random";
      tf_data.calltf = sys_random_calltf;
      tf_data.compiletf = 0;
      tf_data.sizetf = sys_random_sizetf;
      tf_data.user_data = "$random";
      vpi_register_systf(&tf_data);

      tf_data.type   = vpiSysFunc;
      tf_data.tfname = "$dist_poisson";
      tf_data.calltf = sys_dist_poisson_calltf;
      tf_data.compiletf = 0;
      tf_data.sizetf = sys_dist_poisson_sizetf;
      tf_data.user_data = "$dist_poisson";

      tf_data.type   = vpiSysFunc;
      tf_data.tfname = "$dist_uniform";
      tf_data.calltf = sys_dist_uniform_calltf;
      tf_data.compiletf = 0;
      tf_data.sizetf = sys_dist_uniform_sizetf;
      tf_data.user_data = "$dist_uniform";
      vpi_register_systf(&tf_data);
}

/*
 * $Log: sys_random.c,v $
 * Revision 1.11  2004/06/09 22:14:10  steve
 *  Move Mersenne Twister to $mti_random, and make
 *  the standard $random standard. Also, add $dist_poisson.
 *
 */