From: Adam Krolnik (adamk@cyrix.com)
Date: Tue May 26 1998 - 21:59:25 PDT
Behavioral Task Force - Enhancement Request
Assigned Enhancement Request Number: B22
Enhancement Name (Description): Define an algorithm for $random()
Date Submitted: 970418
Requestor: adamk@cyrix.com (Adam Krolnik)
<p>I would like to propose either drand48() which may be a public domain algorithm, or
BSD random() which has the BSD Copyright requirements on it.
The specification of a standard algorithm for $random ($dist_uniform and the other
random distribution functions) allows more portability between simulators for
users choosing to use multiple simulators. It also reduces the risk of a simulator
using an algorithm that has less desirable random properties.
Tom Fitspatrick has volunteered Cadence's algorithm.
Proposal:
[I have removed the DEBUG and LINT preprocessor sections.
What should the copyright say?]
Add section 14.10.3 "Algorithm for Probabilistic distribution functions"
"The algorithm for these functions shall be defined by the C code below.
The verilog function is listed with the corresponding C function.
$dist_uniform rtl_dist_uniform
$dist_normal rtl_dist_normal
$dist_exponential rtl_dist_exponential
$dist_poisson rtl_dist_poisson
$dist_chi_square rtl_dist_chi_square
$dist_t rtl_dist_t
$dist_erlang rtl_dist_erlang
$random rtl_dist_uniform( seed, LONG_MIN, LONG_MAX )
<p><p>/*
* (c) Copyright 1995, Cadence Design Systems, Inc.
* All Rights Reserved.
*/
#include <limits.h>
<p>static double uniform(ARG3( long *seed, long start, long end ));
static double normal(ARG3( long *seed, long mean, long deviation));
static double exponential(ARG2( long *seed, long mean));
static long poisson(ARG2( long *seed, long mean));
static double chi_square(ARG2( long *seed, long deg_of_free));
static double t(ARG2( long *seed, long deg_of_free));
static double erlangian(ARG3( long *seed, long k, long mean));
<p><p>long
rtl_dist_chi_square( seed, df )
long *seed;
long df;
{
double r;
long i;
if(df>0)
{
r=chi_square(seed,df);
if(r>=0)
{
i=(long)(r+0.5);
}
else
{
r = -r;
i=(long)(r+0.5);
i = -i;
}
}
else
{
/*("WARNING [CSDPDF]: Chi_square distribution must have positive degree of freedom\n"); */
i=0;
}
return (i);
}
long
rtl_dist_erlang( seed, k, mean )
long *seed;
long k, mean;
{
double r;
long i;
if(k>0)
{
r=erlangian(seed,k,mean);
if(r>=0)
{
i=(long)(r+0.5);
}
else
{
r = -r;
i=(long)(r+0.5);
i = -i;
}
}
else
{
/*("WARNING [KEDPK]: k-stage erlangian distribution must have positive k\n");*/
i=0;
}
return (i);
}
long
rtl_dist_exponential( seed, mean )
long *seed;
long mean;
{
double r;
long i;
if(mean>0)
{
r=exponential(seed,mean);
if(r>=0)
{
i=(long)(r+0.5);
}
else
{
r = -r;
i=(long)(r+0.5);
i = -i;
}
}
else
{
/*("WARNING [EXDMPM]: Exponential distribution must have a positive mean\n"); */
i=0;
}
return (i);
}
long
rtl_dist_normal( seed, mean, sd )
long *seed;
long mean, sd;
{
double r;
long i;
r=normal(seed,mean,sd);
if(r>=0)
{
i=(long)(r+0.5);
}
else
{
r = -r;
i=(long)(r+0.5);
i = -i;
}
return (i);
}
long
rtl_dist_poisson( seed, mean )
long *seed;
long mean;
{
long i;
if(mean>0)
{
i=poisson(seed,mean);
}
else
{
/*("WARNING [PDMPM]: Poisson distribution must have a positive mean\n");*/
i=0;
}
return (i);
}
long
rtl_dist_t( seed, df )
long *seed;
long df;
{
double r;
long i;
if(df>0)
{
r=t(seed,df);
if(r>=0)
{
i=(long)(r+0.5);
}
else
{
r = -r;
i=(long)(r+0.5);
i = -i;
}
}
else
{
/* ("WARNING [TFPDF]: t distribution must have positive degree of freedom\n");*/
i=0;
}
return (i);
}
long
rtl_dist_uniform(seed, start, end)
long *seed;
long start, end;
{
double r;
long i;
if (start >= end) return(start);
if (end != LONG_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!=LONG_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);
}
}
<p> return (i);
}
tatic double
uniform( seed, start, end )
long *seed, start, end;
{
union u_s
{
float s;
unsigned stemp;
} u;
double d = 0.00000011920928955078125;
double a,b,c;
if ((*seed) == 0)
*seed = 259341593;
if (start >= end)
{
a = 0.0;
b = 2147483647.0;
}
else
{
a = (double) start;
b = (double) end;
}
*seed = 69069 * (*seed) + 1;
u.stemp = *seed;
#if defined(HP9000) || defined(SUN4) || defined(SUN4V)
u.stemp = (u.stemp >> 9) | 0x3f800000;
#endif /* SUN ... */
/* Does RS6000 actually use a different floating point format or was
* this intended for some other IBM machine?
*/
#ifdef RS6000
u.stemp = (u.stemp >> 8) | 0x40000000;
#endif /* RS6000 */
c = (double) u.s;
#ifdef RS6000
c = c + 1.0;
#endif
c = c+(c*d);
c = ((b - a) * (c - 1.0)) + a;
return (c);
}
<p>static double
normal(seed,mean,deviation)
long *seed,mean,deviation;
{
double v1,v2,s;
double log(), sqrt();
s = 1.0;
while((s >= 1.0) || (s == 0.0))
{
v1 = uniform(seed,-1,1);
v2 = uniform(seed,-1,1);
s = v1 * v1 + v2 * v2;
}
s = v1 * sqrt(-2.0 * log(s) / s);
v1 = (double) deviation;
v2 = (double) mean;
return(s * v1 + v2);
}
tatic double
exponential(seed,mean)
long *seed,mean;
{
double log(),n;
n = uniform(seed,0,1);
if(n != 0)
{
n = -log(n) * mean;
}
return(n);
}
tatic long
poisson(seed,mean)
long *seed,mean;
{
long n;
double p,q;
double exp();
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);
}
tatic double
chi_square(seed,deg_of_free)
long *seed,deg_of_free;
{
double x;
long k;
if(deg_of_free % 2)
{
x = normal(seed,0,1);
x = x * x;
}
else
{
x = 0.0;
}
for(k = 2;k <= deg_of_free;k = k + 2)
{
x = x + 2 * exponential(seed,1);
}
return(x);
}
tatic double
t(seed,deg_of_free)
long *seed,deg_of_free;
{
double sqrt(),x;
double chi2 = chi_square(seed,deg_of_free);
double div = chi2 / (double)deg_of_free;
double root = sqrt(div);
x = normal(seed,0,1) / root;
return(x);
}
tatic double
erlangian(seed,k,mean)
long *seed,k,mean;
{
double x,log(),a,b;
long i;
x=1.0;
for(i=1;i<=k;i++)
{
x = x * uniform(seed,0,1);
}
a=(double)mean;
b=(double)k;
x= -a*log(x)/b;
return(x);
}
This archive was generated by hypermail 2.1.4
: Mon Jul 08 2002 - 12:52:53 PDT
and
sponsored by Boyd Technology, Inc.