 /*
  * Khoros: $Id: random.c,v 1.2 1992/03/20 23:42:57 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: random.c,v 1.2 1992/03/20 23:42:57 dkhoros Exp $";
#endif

 /*
  * $Log: random.c,v $
 * Revision 1.2  1992/03/20  23:42:57  dkhoros
 * VirtualPatch5
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1990, University of New Mexico.  All rights reserved.
 * 
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as too the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *----------------------------------------------------------------------
 */

#include "unmcopyright.h"	 /* Copyright 1990 by UNM */
#include "vinclude.h"



/* >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
   >>>>                                                       <<<<
   >>>>                General Math Functions                 <<<<
   >>>>                                                       <<<<
   >>>>			     vfact()                          <<<<
   >>>>			     vimpulse()                       <<<<
   >>>>			     vstep()                          <<<<
   >>>>			     vsign()                          <<<<
   >>>>			     vgnoise()                        <<<<
   >>>>			     vunoise()                        <<<<
   >>>>			     urng()                           <<<<
   >>>>			     vrandom()                        <<<<
   >>>>			     vsrandom()                       <<<<
   >>>>                                                       <<<<
   >>>>>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  */



/************************************************************
*
*  MODULE NAME: vfact
*
*      PURPOSE: computes x factorial (x!)
*
*        INPUT: double precision x
*
*       OUTPUT: Returns x factorial
*               x! = 1 * 2 * 3 * ... * (x-1) * x
*
*    CALLED BY: xvparser & xvexpr & vpoly
*
*   WRITTEN BY: Mark Young
*
*************************************************************/

double vfact(x)
double x;
{
	int i, factor = x;

	for  (i = factor-1; i > 0; i--) 
	   factor *= i;

	return((double) factor);
}


/************************************************************
*
*  MODULE NAME:  vimpulse
*
*      PURPOSE:  impulse function
*
*        INPUT:  double precision x
*
*       OUTPUT:  impulse function of x
*
*    CALLED BY:  xvparser & xvexpr & vpoly
*
*   WRITTEN BY:  Marcelo  Teran
*
*
*************************************************************/

double vimpulse (x)
double x;
{
	static	int	last,actual;

	if ((x > -0.00001) && (x < 0.00001))
	   return (1.0);

	if(x > 0) actual = 1;
	if(x < 0) actual = 0;

	if (last != actual)
	{
	   last = actual;
	   return(1.0);
	}
	else
	{
	   last = actual;
	   return(0.0);
	}
}


/************************************************************
*
*  MODULE NAME:  vstep
*
*      PURPOSE:  gives a step function
*
*        INPUT:  double precision x
*
*       OUTPUT:  gives a step function based on value of x
*
*    CALLED BY:  xvparser & xvexpr & vpoly
*
*   WRITTEN BY:  Marcelo  Teran
*
*************************************************************/

double vstep (x)
double x;
{
	if (x < 0.0)
	   return (0.0);
	else
	   return (1.0);
}

/************************************************************
*
*  MODULE NAME:  vsign
*
*      PURPOSE:  sign function
*
*        INPUT:  double precision x
*
*       OUTPUT:  returns (-1, 0, or 1) based on sign of argument
*
*    CALLED BY:  xvparser & xvexpr & vpoly
*
*   WRITTEN BY:  Tom Sauer
*
*************************************************************/

double vsign (x)
double x;
{
	if (x < 0.0)
	   return (-1.0);
        else if (x > 0.0)
           return(1.0);
        else
          return(0.0);
}


/***********************************************************
*
*       MODULE NAME: vgnoise
*
*	    PURPOSE: generates a gaussian randon number.
*
*	      INPUT: double precision x
*
* 	    OUTPUT : Returns the value of one point of the 
*		     gaussian noise being generated.
*
*        CALLED BY : xvparser & xvexpr & vpoly
*
*        WRITTEN BY: Marcelo  Teran
*
************************************************************/

double vgnoise(x)

double	x;
{
	int	i;
	double	y;
	float   urng();

	for (y = 0.0, i = 0; i < 12; i++)
	    y = y + ((double)urng() - 0.5);

	return(y);
}


/***********************************************************
*
*       MODULE NAME: vunoise
*
*	    PURPOSE: generates a uniform randon number with 
*		     variance 1 and mean 0.
*
*	      INPUT: double precision x
*
* 	    OUTPUT : Returns the value of one point of the 
*		     gaussian noise being generated.
*
*        CALLED BY : xvparser & xvexpr & vpoly
*
*        WRITTEN BY: Marcelo  Teran
*
************************************************************/

double vunoise(x)
double	x;
{
	double	mean,variance;
	float   urng();

	mean = 0.0;	variance = 1.0;
	return(2.0*sqrt(3.0*variance)*((double) urng()-0.5)+mean);
}


/***********************************************************************
*
*   MODULE NAME: urng()
*
*       PURPOSE: generates a uniform random number (float)
*
*         INPUT: none
*
*        OUTPUT: a float random number in the range [0,1] inclusive
*
*        CALLED BY : xvparser & xvexpr & vpoly
*
*    WRITTEN BY: Ramiro Jordan, Marcelo Teran  
*
* MODIFICATIONS: April 23, 1990 - (Jeremy Worley) cleaned it up, 
*                                 fixed up scaling.
*
***********************************************************************/

float urng()
{
  long int i, big, vrandom();
  float f;
  static int iseed = 0;
  int seed;
/*
 * use time as seed 
 */

  if (iseed == 0) {
      vsrandom( _genseed() ); 
      iseed = 1;
  } 

/*
 * generate random number
 */

  i = vrandom(&big);

/*
 * scale it down to range [0,1]
 */

  f = (float)((double)i/(double)big);

  return(f);
}

/***********************************************************************
*
*   MODULE NAME: vsrandom()
*
*       PURPOSE: set the seed for the random number generator
*
*         INPUT: none
*
*        OUTPUT: none
*
*    WRITTEN BY: Mark Young & Tom Sauer
*
* MODIFICATIONS: 
*
***********************************************************************/


vsrandom(seed)

int seed;
{
#ifdef RANDOM
	srandom(seed);
#else
	srand(seed);
#endif
}

/***********************************************************************
*
*   MODULE NAME: vrandom()
*
*       PURPOSE: generates random number
*
*         INPUT: none
*
*        OUTPUT: a random number 
*
*    WRITTEN BY: Mark Young & Tom Sauer
*
* MODIFICATIONS: 
*
***********************************************************************/


#ifdef RANDOM
#define BIG 	2147483647		/* 2^31 -1 */
#else
#define BIG 	     32767		/* 2^15 -1 */
#endif

long vrandom(big)

long *big;
{
#ifdef RANDOM
	long random(), num;
#else
	long rand(), num;
#endif

	if (big != NULL)
	   *big = BIG;

#ifdef RANDOM
	num = random();
#else
	num = rand();
#endif
	return(num);
}


/***********************************************************************
*
*   MODULE NAME: _genseed()
*
*       PURPOSE: generates the seed for the random number generator
*
*         INPUT: none
*
*        OUTPUT: none
*
*    WRITTEN BY: Jeremy Worley 
*
* MODIFICATIONS: 
*
***********************************************************************/

#define PRIME1 4001
#define PRIME2 4091

int _genseed()
{
  int t,p;

  t = (int)time(NULL) % PRIME1;

  p = (int)getpid() % PRIME2; 

  return((int)(t*p + p));
}
