/*
 * 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.14 2004/10/04 01:10:58 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);

	  /* NOTE: The cast of r to i can overflow and generate
	     strange values, so cast to unsigned long
	     first. This eliminates the underflow and gets the
	     twos complement value. That in turn can be cast
	     to the long value that is expected. */

	if (end != UNIFORM_MAX)
	{
		end++;
		r = uniform( seed, start, end );
		if (r >= 0)
		{
			i = (unsigned long) r;
		}
		else
		{
			i = (unsigned 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 = (unsigned long) r;
		}
		else
		{
			i = (unsigned 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 = (unsigned long) r;
		}
		else
		{
			i = (unsigned 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;
      static 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.14  2004/10/04 01:10:58  steve
 *  Clean up spurious trailing white space.
 *
 * Revision 1.13  2004/06/17 14:44:01  steve
 *  Save seed in static variable, in case user doesnt pass it.
 *
 * Revision 1.12  2004/06/10 02:14:42  steve
 *  Fix transcription error scaling c in uniform range.
 *
 * 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.
 *
 */

