From qiclab!mdcgwy.mdc.com!HANNA%VAXEC1.decnet Tue Nov 23 14:23:01 1993
Date: 23 Nov 93 13:02:00 CST
From: "VAXEC1::HANNA" <HANNA%VAXEC1.decnet@mdcgwy.mdc.com>
Subject: Random Number?Random Vector Generator! HAPPY THANKSGIVING!!!!
To: "vhdlsynth" <vhdlsynth@vhdl.org>
Content-Length: 56574
Date: 11/23/1993
Subj: Free Distribution of VHDL_Based Rand Number/Random
Vector Generation (RNG/RVG) S/W
From: Bill Hanna
McDonnell Douglas Aerospace- Defense & Electronics Systems
MDA-D&ES
MC 500 4224
P.O. Box 426
St. Charles, MO 63302, USA
1. The Following is a random number generator package and a test bench
to generate random vectors: 1s and 0s. This code was developed for
internal use within MDC (McDonnell Douglas Corporation). I am making
this code available to memebers of IEEE-DASC with permission from
MDC without liability to MDC.
2. MDC and the writer do not assume any responsibility for the
correctness and/or accuracy of this S/W package. Use is according to
restrictions in the first paragraph of the S/W Package.
3. We developed the S/W originally for use with VHDL tools operating under
UNIX. The attached copy is a modified version which eliminates all UNIX Math.
Library Calls by encoding math. functions in VHDL. This might not be as
precise mathematically as a UNIX Production Quality Math. Library, but it
eliminates the need for Library Calls; Self contained VHDL. We ran this
package on several VHDL applications with positive results.
3. Please feel free to use it. Also, we would appreciate some feedback
from users. We will try to respond to questions as soon as we can, but
do not expect a vendor type of service since we are not in the S/W
business. Good Luck!
Bill Hanna
hanna%vaxec1.decnet@mdcgwy.mdc.com
--Subject: Rnd testbench
-- Here's the main testbench I used to check out the routines. It generates
-- 100 Test Vectors (18-bit wide) based on 1800
-- random values and writes them to the file 'rnd.out' and to Output. It also
-- calculates the mean and standard deviation.
-- To change the distribution, just change the initial values of the fields of
-- the variable rnd_rec and recompile.
--
-- To Seed the RNG, Use a File "XYZ.IN" Provided or your own seed
--
---------------------------------------- CUT HERE ----------------------------
-- **********************************************************************
-- * This Software Package was developed for internal use within *
-- * McDonnell Douglas Corporation (MDC). MDC and the authors or *
-- * distributers of this code are not responsible for correctnes or *
-- * accuracy of the procedures/functions in this package. This code *
-- * may not be distributed for commercial purposes. IN NO EVENT SHALL *
-- * MDC BE LIABLE FOR ANY INCIDENTAL, INDIRECT , SPECIAL, OR *
-- * CONSEQUENTIAL DAMAGES WHATSOEVER ARISING OUT OF OR RELATING TO *
-- * THE USE OF THE INFORMATION/SOFTWARE PROVIDED IN THIS PACKAGE. *
-- **********************************************************************
Use Std.TextIO.All;
Library Utility2;
-- Use Utility.c_math.all; -- This Line Calls Unix Math Utility.
Use Utility2.rnd2.All; -- Utility2 Is our Own VHDL Coded Math.
Entity rnd_test2 Is
End rnd_test2;
Architecture a2 of rnd_test2 Is
Constant imax : integer:= 100; -- Number of Vectors
Constant jmax : integer:= 18; -- Number of Bits Per Vector
Signal ATGV: bit_vector(0 to jmax-1); -- Random ATG Vector
File f: Text Is Out "rnd.out";
File file_f: Text Is In "XYZ.IN";
Begin
Process
Constant str1: String(1 to 4) := "RN: "; -- Random Number
Constant str2: String(1 to 11) := "Seed(3-0): ";
Constant str3: String(1 to 5) := "Avg: ";
Constant str4: String(1 to 12) := ", Std Dev: ";
Constant str6: String(1 to 13) :="TEST_VECTOR: ";
Constant bound_h: Real :=10.0;
Constant bound_l: Real := -10.0;
Variable rnd_rec: rnd_rec_t := (
distribution => rnd_empirical_c,
seed => (0, 0, 0, 123),
--seed => (367, 23, 4, 191),
--seed => (250, 1843, 3687, 991),
bound_l => -10.0,
bound_h => 10.0,
trials => 10,
p_success => 0.7,
mean => 2.0,
std_dev => 1.0,
others => 0.0
);
Variable emp: rnd_empirical_vector_t(0 to 5) := (
(0.0, 0.0),
(1.0, 0.1),
(2.0, 0.3),
(3.0, 0.8),
(4.0, 0.8),
(5.0, 1.0)
);
Variable l,fl: Line; -- Peculiar Way of Writing Twice
Variable m,fm: LIne; -- Special Line for Test_Vector
Variable i: Integer;
Variable avg, sd, sum_of_squares: Real;
Type TV1553T Is array (1 to imax, 0 to jmax-1) of Bit;
Type TEMPVT Is array (0 to jmax-1) of Bit;
Variable TV1553: TV1553T;
Variable TEMP_V: bit_Vector (0 to jmax-1);
Begin
avg := 0.0;
sd := 0.0;
sum_of_squares := 0.0;
Rnd_read_Record (rnd_rec); -- Testing Read Seed From File
-- Comment the Previous Line if you do not want to use File for Seeding!
-- Seed File Name is "XYZ.IN"
For i In 1 to imax Loop -- Outer Loop
Wait For 10 fs; -- attach time to test vectors
Write(l, str1);
For j in 0 to jmax-1 Loop -- Inner Loop Takes care of One Vector
Rnd_Random(rnd_rec, emp);
--Write(l, rnd_rec.rnd);
--WriteLine(Output, l);
--Write(l, integer(rnd_rec.rnd), field => 10);
Write(l, rnd_rec.rnd);
Write(l, ' '); -- Space
-- WriteLine(f, l);
--Write(l, integer(rnd_rec.rnd), field => 10);
--Write(l, rnd_rec.rnd);
--WriteLine(Output, l);
--Write(l, str2);
--Write(l, rnd_rec.seed(3), field => 5);
--Write(l, rnd_rec.seed(2), field => 5);
--Write(l, rnd_rec.seed(1), field => 5);
--Write(l, rnd_rec.seed(0), field => 5);
--WriteLine(Output, l);
avg := avg + rnd_rec.rnd;
sum_of_squares := sum_of_squares + rnd_rec.rnd * rnd_rec.rnd;
-- Create One 1553 Vector
if rnd_rec.rnd < 2.5 then TV1553(i,j):='1';
else TV1553(i,j):='0';
end if;
-- Set a bit at random!
TEMP_V (j) := TV1553(i,j); -- Copy It to TEMP_V
End Loop; -- Ends Inner Loop
fl:=l; WriteLine(Output, l);
WriteLine(f, fl);
write(m, str6); -- Write Label "TEST_VECTOR"
write(m, TEMP_V); -- Write Test Vector
ATGV <= TEMP_V; -- TEST VECTOR!!!!
fm:=m; writeLine(Output, m); -- Print test vector
writeLine(f, fm); -- Write test vector to File F
End Loop;
sd := SqRt((sum_of_squares * Real(imax * jmax) - avg * avg) /
(Real(imax *jmax * (imax *jmax- 1))));
avg := avg / Real(imax * jmax);
Write(l, str3);
Write(l, avg);
Write(l, str4);
Write(l, sd);
fl:=l; WriteLine(Output, l);
writeLine(f, fl);
End Process;
End a2;
-- **********************************************************************
-- * This Software Package was developed for internal use within *
-- * McDonnell Douglas Corporation (MDC). MDC and the authors or *
-- * distributers of this code are not responsible for correctnes or *
-- * accuracy of the procedures/functions in this package. This code *
-- * may not be distributed for commercial purposes. IN NO EVENT SHALL *
-- * MDC BE LIABLE FOR ANY INCIDENTAL, INDIRECT , SPECIAL, OR *
-- * CONSEQUENTIAL DAMAGES WHATSOEVER ARISING OUT OF OR RELATING TO *
-- * THE USE OF THE INFORMATION/SOFTWARE PROVIDED IN THIS PACKAGE. *
-- **********************************************************************
Package rnd2 Is
--------------------------------------------------------------------------
-- --
-- Random Number Generation Package --
-- --
-- Authors: --
-- John A. Breen, Ken Christensen --
-- McDonnell Douglas Missile Systems Company, July 1989 --
-- Dept E435, Bldg 111/2/N-3, (314)234-4341 --
-- Updates/Distribution: --
-- William A. Hanna --
-- McDonnell Douglas Aerospace-Defense & Electronic Systems, Nov. 1993 --
-- MC 500 4224 --
-- MDA-D&ES --
-- P.O.Box 426 --
-- St. Charles, MO 63302 --
-- --
-- --
-- Date: 12-July-1989 --
-- --
-- This package contains routines for generating random numbers from --
-- 14 different discrete and uniform distributions. Each distribution --
-- is accessed through a single procedure call, Rnd_Random. All --
-- parameters to this procedure are passed via a record of type --
-- rnd_rec_t, with the exception of the empirical array, which is --
-- optional (We would have liked to include the array in the record, --
-- but VHDL does not allow record fields to be unconstrained arrays). --
-- It would have been nice if we could have made Rnd_Random a function --
-- instead of a procedure, but the subprogram has to modify the seed, --
-- and VHDL doesn't allow functions to modify their parameters. --
-- --
-- In addition, there is a procedure for loading a rnd_rec_t record --
-- from an ASCII file; this routine is named Rnd_Read_Record, and takes--
-- two parameters: the file variable, and the record to be loaded. --
-- --
-- The basic random number generator, from which all of the various --
-- distributions are derived, is a mixed linear congruential generator --
-- with a 48-bit seed (in fact, it emulates the UNIX System V function --
-- erand48). To make it easier for the user to choose independent --
-- initial seeds, a constant array, rnd_seeds, is provided, which the --
-- user can specify directly in his models (although Rnd_Read_Rec --
-- isn't smart enough to recognize it in an ASCII file, at least not --
-- yet). --
-- --
-- Dependencies: Removed in this version! --
-- Package TEXTIO --
-- A math package (in this implementation called C_MATH) which --
-- contains the following functions: --
-- Log(x) (natural log) --
-- Ceil(x) (smallest integral value >= x) --
-- Floor(x) (largest integral value <= x) --
-- SqRt(x) (exponentiation with Real exponent) --
-- "**"(x, y) (exponentiation with Real exponent) --
-- Each of these functions takes Real arguments and returns a --
-- Real. We cheated a little by "hooking" in the UNIX math --
-- library, but there was no easy way around it. --
-- --
-- This package was developed for MDC internal use only. Public Use --
-- Is Subject to Resrictions as Stated In the First Paragraph of this --
-- Document. --
--------------------------------------------------------------------------
-- Use Std.textIO.Text; -- Some Tools Require This Line
Use STD.TEXTIO.ALL; -- Other VHDL Tolls Require This Instead
Type rnd_distribution_t Is (
--
-- Discrete distributions
--
rnd_constant, -- Constant value (always returns the mean)
rnd_uniform_d, -- Integer uniform, [bound_l, bound_h)
rnd_poisson,
rnd_binomial,
rnd_geometric,
rnd_empirical_d, -- User-specified array (cumulative distribution)
--
-- Continuous distributions
--
rnd_uniform_c, -- Real uniform, [bound_l, bound_h)
rnd_normal,
rnd_triangular,
rnd_exponential,
rnd_hyperexponential,
rnd_gamma,
rnd_erlang,
rnd_empirical_c -- User-specified array (cumulative distribution)
);
Type rnd_seed_t Is Array (3 downto 0) of Integer;
Type rnd_seed_vector_t Is Array (Positive Range <>) of rnd_seed_t;
Type rnd_empirical_t Is Record
x: Real; -- Value of random variable
p: Real; -- cumulative probability of value x
End Record;
Type rnd_empirical_vector_t Is Array (Natural Range <>) of rnd_empirical_t;
Type rnd_rec_t Is Record
rnd: Real;
distribution: rnd_distribution_t;
seed: rnd_seed_t;
mean, std_dev: Real;
bound_l, bound_h: Real;
trials: Integer;
p_success: Real;
End Record;
File file_f: Text Is in "XYZ.IN"; -- Special Variable for file_f
-- ------------------------------------------------------------------------
-- VHDL CODING OF MATH FUNCTIONS.
-- WE RECOMMEND THE USE OF A MATH LIBRARY IF ONE IS AVAILABLE
-- SQRT Is A Square Root Polynomial Type Function
Function SQRT (X:Real) return Real;
-- Natural Logarithm
Function Log (X:Real) return Real;
-- Exponential Function
Function Exp (X: Real) return Real;
-- Exponential Function #2
Function "**" (X,Y: Real) return Real;
-- Floor Function
Function Floor (X:Real) return Real;
-- Ceil Function
Function Ceil (X:Real) return Real;
-- ---------------------------------------------------------------------
-- Procedure Rnd_Read_Record(file_f: Text; rnd_rec: inout rnd_rec_t);
Procedure Rnd_Read_Record(rnd_rec: inout rnd_rec_t);
Procedure Rnd_Random(
rnd_rec: InOut rnd_rec_t;
empirical_data: In rnd_empirical_vector_t);
-- V-Systems Does Not Like Initial Values
-- := (0 => (0.0, 0.0), 1 => (1.0, 1.0)));
Constant rnd_seeds: rnd_seed_vector_t(1 to 50) := (
( 80, 2466, 2681, 2459), ( 62, 493, 536, 699), (1430, 2639, 3931, 3035),
( 153, 2415, 2756, 1275), (3114, 3364, 3281, 3611), (2389, 3195, 3588, 1851),
(3616, 2485, 1756, 91), (4015, 1093, 4056, 2427), (4022, 2774, 378, 667),
(1343, 3394, 1088, 3003), (2234, 3739, 172, 1243), (2914, 1828, 3900, 3579),
( 494, 1624, 2162, 1819), (2711, 1250, 1229, 59), (1641, 1698, 3276, 2395),
(3559, 3246, 2289, 635), (2303, 1870, 442, 2971), ( 586, 2045, 4009, 1211),
(3491, 878, 2876, 3547), (2567, 900, 3317, 1787), (1738, 2392, 3411, 27),
(1050, 3895, 1237, 2363), (1746, 2720, 3069, 603), (1947, 3657, 2889, 2939),
(1379, 3098, 2875, 1179), ( 797, 1840, 1105, 3515), (2261, 1498, 3853, 1755),
( 20, 927, 1005, 4091), ( 914, 821, 2931, 2331), (2618, 138, 3614, 571),
(1404, 703, 1133, 2907), (2801, 3623, 1762, 1147), ( 336, 1611, 3579, 3483),
(2219, 4073, 570, 1723), (1547, 751, 3101, 4059), (2633, 3204, 1062, 2299),
(3148, 255, 724, 539), (3713, 3565, 166, 2875), (1977, 3086, 1566, 1115),
(1860, 343, 3002, 3451), (3174, 753, 2556, 1691), (2333, 1848, 2402, 4027),
(2098, 1981, 622, 2267), (3740, 835, 3487, 507), ( 50, 4037, 884, 2843),
(4019, 3181, 3183, 1083), ( 531, 926, 270, 3419), (3360, 3306, 2515, 1659),
(2425, 3867, 3900, 3995), (4045, 558, 2507, 2235)
);
End rnd2;
Library Utility2;
-- Use Utility.C_Math.All;
Use Utility2.all;
Package Body rnd2 Is
Constant rnd_lcg_a: rnd_seed_t := (16#0#, 16#5DE#, 16#ECE#, 16#66D#);
Constant rnd_lcg_c: rnd_seed_t := (16#0#, 16#0#, 16#0#, 16#B#);
Constant rnd_lcg_m: Integer := 2**12;
-- ----------------------------------------------------------------------
-- Body of Math Functions within VHDL!
-- VHDL CODING OF MATH FUNCTIONS. WE RECOMMEND THE USE OF A MATH LIBRARY
-- IF ONE IS AVAILABLE
-- SQRT Is A Square Root Polynomial Type Function
Function SQRT (X:Real) return Real Is
Variable A: Real; -- TEMP for Evaluation
begin
A:=0.25497+X*(0.11722-X*(0.67206+X*(0.30750-X*(0.62662))));
return A;
end;
-- -------------------------------------------------------------------
-- Natural Logarithm
Function Log (X:REal) return Real Is
Variable A: Real := 0.0; -- TEMP for Evaluation
begin
A :=0.2+X*(0.66666+X*(0.40020+X*(0.28233+X*(0.25095
-X*(0.55393+X*(0.41272 ))))));
return A;
end;
-- ---------------------------------------------------------------------
-- Exponential Function
Function Exp (X: Real) return Real Is
Variable A: Real :=1.0; -- TEMP for Evaluation
begin
A:= 1.0+X*(1.0+X*(1.0/2.0+X*(1.0/6.0+X*(1.0/24.0))));
return A;
end;
-- ---------------------------------------------------------------------
-- Exponential Function #2
Function "**" (X,Y: Real) return Real Is
Variable A: Real:=1.0; -- TEMP for Evaluation
begin
A := exp ( Y * Log (X));
return A;
end;
-- ---------------------------------------------------------------------
-- Floor Function
Function Floor (X:Real) return Real Is
Variable A: Real:=0.0; -- TEMP for Evaluation
begin
for i in -32768.0 to 32767.0 loop
if X>Real(i) then A:=Real(i)-1.0; return A;
else A:=Real(i);
end if;
end loop;
return A;
end;
-- -------------------------------------------------------------------
-- Ceil Function
Function Ceil (X:Real) return Real Is
Variable A: Real:=0.0; -- TEMP for Evaluation
begin
for i in -32768.0 to 32767.0 loop
if X>Real(i) then A:=Real(i)+1.0; return A;
else A:=Real(i);
end if;
end loop;
return A;
end;
-- -------------------------------------------------------------------
Procedure rnd48(rnd: Out Real; seed: InOut rnd_seed_t) Is
------------------------------------------------------------------
--Rnd48 : Receives a 48 bit number, the seed, and returns a --
-- random number between 0 and 1. --
-- --
--Reference : Domain/IX Programmer's Reference for System V. --
-- Apollo Computer Inc. --
-- --
--Used Fields : rnd, seed. --
-- --
--Changed Fields : rnd, seed. --
------------------------------------------------------------------
Variable i : integer;
Function Mult(rnd_lcg_b : rnd_seed_t) return rnd_seed_t is
Variable temp : rnd_seed_t;
Variable i, j : integer;
Begin
for i in 0 to 3 loop
temp(i) := 0;
end loop;
for i in 0 to 3 loop
for j in 0 to (3 - i) loop
temp(i + j) := temp(i + j) + (rnd_lcg_a(i) * rnd_lcg_b(j));
end loop;
end loop;
for i in 0 to 2 loop
if temp(i) >= rnd_lcg_m then
temp(i + 1) := temp(i + 1) + ( temp(i) / rnd_lcg_m);
temp(i) := temp(i) mod rnd_lcg_m;
end if;
end loop;
temp(3) := temp(3) mod rnd_lcg_m;
return temp;
End Mult;
Function Add(x : rnd_seed_t) return rnd_seed_t is
Variable tempa : rnd_seed_t;
Variable i : integer;
Begin
for i in 0 to 3 loop
tempa(i) := x(i) + rnd_lcg_c(i);
end loop;
for i in 0 to 2 loop
if tempa(i) >= rnd_lcg_m then
tempa(i +1) := tempa(i +1) + (tempa(i) / rnd_lcg_m);
tempa(i) := tempa(i) mod rnd_lcg_m;
end if;
end loop;
tempa(3) := x(3) mod rnd_lcg_m;
return tempa;
End Add;
Begin
For i in 0 to 3 loop
If seed(i) >= rnd_lcg_m then
Assert False
Report "The seed element must be less than the 4096"
Severity Warning;
seed(i) := seed(i) mod rnd_lcg_m;
Elsif
Seed(i) < 0 then
Assert False
Report "The seed must be a positive number"
Severity Warning;
seed(i) := (-seed(i)) mod rnd_lcg_m;
End if;
End loop;
seed := add(mult(seed));
rnd := ((((Real(seed(0)) / Real(rnd_lcg_m)
+ Real(seed(1))) / Real(rnd_lcg_m) + Real(seed(2))) /
Real(rnd_lcg_m) + Real(seed(3))) / Real(rnd_lcg_m));
End rnd48;
Procedure uniform_d(rnd_rec: InOut rnd_rec_t) Is
-----------------------------------------------------------------
--uniform_d : Generate a uniformly distributed integer --
-- random number in the range [bound_l, bound_h). --
-- --
--Used Fields : seed, bound_l, bound_h. --
-- --
--Changed Fields : rnd, seed. --
-----------------------------------------------------------------
Variable tmp : real;
Begin
If (rnd_rec.bound_l > rnd_rec.bound_h) then
Assert False
Report "Lower bound is greater than the upper bound."
Severity Warning;
rnd_rec.rnd := rnd_rec.bound_h;
Elsif (rnd_rec.bound_l /= real(integer(rnd_rec.bound_l))) or
(rnd_rec.bound_h /= real(integer(rnd_rec.bound_h))) then
Assert False
Report "The lower and upper bounds are not integer numbers"
Severity Warning;
rnd_rec.rnd := floor(rnd_rec.bound_h);
Else
rnd48(tmp, rnd_rec.seed);
rnd_rec.rnd := rnd_rec.bound_l + floor
((rnd_rec.bound_h + 1.0 - rnd_rec.bound_l) * tmp);
End if;
End uniform_d;
Procedure poisson(rnd_rec: InOut rnd_rec_t) Is
----------------------------------------------------------------
--poisson : Generate a random number from the Poisson --
-- distribution. --
-- --
--Reference : Lavenberg, computer performance modeling hand- --
-- book, 1983, section 5.2.1, algorithm 5.4. --
-- --
--Used Fields : mean, seed. --
-- --
--Changed Fields : rnd, seed. --
----------------------------------------------------------------
Variable temp, term, sum, x : real;
Begin
If (rnd_rec.mean > 0.0) then
x := 0.0;
rnd48(temp, rnd_rec.seed);
term := Exp(-rnd_rec.mean);
If (term > 0.0) then
sum := term;
While temp > sum loop
x := x + 1.0;
term := rnd_rec.mean / x * term;
sum := sum + term;
End loop;
rnd_rec.rnd := x;
Else
Assert False
Report " The mean is too large"
Severity Warning;
rnd_rec.rnd := rnd_rec.mean;
End if;
Else
Assert False
Report "The mean must be greater than zero"
Severity Warning;
rnd_rec.rnd := rnd_rec.mean;
End if;
End poisson;
Procedure binomial(rnd_rec: InOut rnd_rec_t) Is
----------------------------------------------------------------
--binomial : Generate a random number from the Binomial --
-- Distribution with the trials having a probability--
-- for success, p_success. --
-- --
--Reference : Lavenburg, computer performance modeling hand- --
-- book, 1983, section 5.2.1, algorithm 5.7. --
-- --
--Used Fields : trials, p_success, seed. --
-- --
--Changed Fields : rnd, seed. --
----------------------------------------------------------------
Variable ngood,temp : real;
Variable i :integer;
Begin
If ((rnd_rec.p_success >= 0.0) and (rnd_rec.p_success <= 1.0)) then
If rnd_rec.trials >= 1 then
ngood := 0.0;
For i in 1 to rnd_rec.trials loop
rnd48(temp, rnd_rec.seed);
If (temp < rnd_rec.p_success) then
ngood := ngood + 1.0;
End if;
rnd_rec.rnd := ngood;
End loop;
Else
Assert False
Report "The number of trials must be a positive integer"
Severity Warning;
rnd_rec.rnd := 0.0;
End if;
Else
Assert False
Report "The probability of success must be between zero and one"
Severity Warning;
rnd_rec.rnd := 0.0;
End if;
End binomial;
Procedure geometric(rnd_rec: InOut rnd_rec_t) Is
----------------------------------------------------------------
--geometric : Generate a random number from the Geometric --
-- Distribution with probability of success, --
-- p_success. --
-- --
--Reference : Lavenberg, computer performance modeling hand- --
-- book, 1983, section 5.2.1, --
-- --
--Used Fields : seed, p_success. --
-- --
--Changed Fields : rnd, seed. --
----------------------------------------------------------------
Variable rgeom, temp : real;
Begin
If ((rnd_rec.p_success > 0.0) and (rnd_rec.p_success <= 1.0)) then
If (rnd_rec.p_success < 1.0) then
rnd48(temp, rnd_rec.seed);
rnd_rec.rnd := ceil(log(temp) / log(1.0 - rnd_rec.p_success));
Else
rnd_rec.rnd := 1.0;
End if;
Else
Assert False
Report "The probability of success must be between zero and one"
Severity Warning;
rnd_rec.rnd := 0.0;
End if;
End geometric;
Procedure empirical_d(rnd_rec: InOut rnd_rec_t;
emp_data: In rnd_empirical_vector_t) Is
----------------------------------------------------------------
--empirical_d : Generate a random number from the cumulative --
-- distribution specified in emp_data. The last --
-- element of this array must have a probability --
-- of 1.0; the first element must have a --
-- probability GREATER than 0. --
-- --
--Used Fields : seed --
-- --
--Changed Fields : rnd, seed --
----------------------------------------------------------------
Variable u: Real;
Variable i: Integer;
Begin
If ((emp_data'Length = 0) Or (emp_data'Left > emp_data'Right)) Then
Assert False
Report "An ascending array of empirical data must be given"
Severity Warning;
rnd_rec.rnd := 0.0;
Elsif (emp_data(emp_data'Left).p <= 0.0) Then
Assert False
Report "The first data point's probability must be greater than 0.0"
Severity Warning;
rnd_rec.rnd := 0.0;
Elsif (emp_data(emp_data'Right).p /= 1.0) Then
Assert False
Report "The last data point's probability must be 1.0"
Severity Warning;
rnd_rec.rnd := 0.0;
Else
rnd48(u, rnd_rec.seed);
i := emp_data'Left;
While (emp_data(i).p < u) Loop
i := i + 1;
End Loop;
rnd_rec.rnd := emp_data(i).x;
End If;
End empirical_d;
Procedure uniform_c(rnd_rec: InOut rnd_rec_t) Is
--------------------------------------------------------------------------
-- uniform_c: Generate a uniformly distributed random number --
-- in the range [bound_l, bound_h) --
-- --
-- Used fields: seed, bound_l, bound_h --
-- --
-- Changed fields: rnd, seed --
--------------------------------------------------------------------------
Variable tmp_rnd: Real;
Begin
If (rnd_rec.bound_l > rnd_rec.bound_h) Then
Assert False
Report "Lower bound is greater than upper bound."
Severity Warning;
rnd_rec.rnd := rnd_rec.bound_h;
Else
rnd48(tmp_rnd, rnd_rec.seed);
rnd_rec.rnd := tmp_rnd * (rnd_rec.bound_h - rnd_rec.bound_l)
+ rnd_rec.bound_l;
End If;
End uniform_c;
Procedure normal(rnd_rec: InOut rnd_rec_t) Is
----------------------------------------------------------------
--normal : Generate a random number from the Normal --
-- Distribution(polar method). --
-- --
--Reference : Lavenberg, computer performance modeling hand- --
-- book, 1983, polar method. --
-- --
--Used Fields : std_dev, seed, mean. --
-- --
--Changed Fields : rnd, seed. --
----------------------------------------------------------------
Variable temp1, temp2, v1, v2, s : real;
Begin
s := 1.0;
While s >= 1.0 loop
rnd48(temp1, rnd_rec.seed);
rnd48(temp2, rnd_rec.seed);
v1 := 2.0 * temp1 - 1.0;
v2 := 2.0 * temp2 - 1.0;
s := v1 * v1 + v2 * v2;
end loop;
If rnd_rec.std_dev < 0.0 then
Assert False
Report "The standard deviation must be a positive number"
Severity Warning;
rnd_rec.rnd := rnd_rec.mean;
Else
rnd_rec.rnd := rnd_rec.std_dev * v1 * sqrt((-2.0 * log(s)) / s)
+ rnd_rec.mean;
End if;
End normal;
Procedure triangular(rnd_rec: InOut rnd_rec_t) Is
----------------------------------------------------------------
--triangle : Generate a random number from the Triangular --
-- Distribution. --
-- --
--Reference : Pritsker, introduction to simulation and --
-- Slam II, 1984, appendix g. --
-- --
--Used Fields : mean, bound_l, bound_h, seed. --
-- --
--Changed Fields : rnd, seed. --
----------------------------------------------------------------
Variable temp : real;
Begin
If (rnd_rec.bound_l < rnd_rec.bound_h) and
(rnd_rec.bound_l < rnd_rec.mean) and
(rnd_rec.mean < rnd_rec.bound_h) then
rnd48(temp, rnd_rec.seed);
If (temp > ((rnd_rec.mean - rnd_rec.bound_l) /
(rnd_rec.bound_h - rnd_rec.bound_l))) then
rnd_rec.rnd := rnd_rec.bound_h - sqrt
((rnd_rec.bound_h - rnd_rec.mean) *
(rnd_rec.bound_h - rnd_rec.bound_l) * (1.0 - temp));
Else
rnd_rec.rnd := rnd_rec.bound_l + sqrt((rnd_rec.mean
- rnd_rec.bound_l) * (rnd_rec.bound_h
- rnd_rec.bound_l) * temp);
End if;
Else
Assert False
Report "Improper lower bound, upper bound, or mean"
Severity Warning;
rnd_rec.rnd := rnd_rec.mean;
End if;
End triangular;
Procedure hyperexponential(rnd_rec: InOut rnd_rec_t) Is
----------------------------------------------------------------
--hyperexponential : Generate a hyperexponentially --
-- Distributed random number. --
-- --
--Used Fields : mean, std_dev, seed. --
-- --
--Changed Fields : rnd, seed. --
----------------------------------------------------------------
Variable coeff, cosq, p, fac, temp, temp1 : real;
Begin
If ((rnd_rec.mean <= 0.0) or (rnd_rec.std_dev <= 0.0)) then
Assert False
Report "The mean and standard deviation must be positive"
Severity Warning;
rnd_rec.rnd := 1.0;
Else
coeff := rnd_rec.std_dev / rnd_rec.mean;
If coeff < 1.0 then
Assert False Report
"The ratio of the deviation to the mean must be greater than 1"
Severity Warning;
rnd_rec.rnd := rnd_rec.mean;
Else
cosq := coeff * coeff;
p := (1.0 + (cosq) - sqrt(cosq * cosq - 1.0)) /
(2.0 * (1.0 + cosq));
If ((p <= 0.0) or (p > 1.0)) then
Assert False
Report "Improper mean and standard deviation"
Severity Warning;
rnd_rec.rnd := rnd_rec.mean;
Else
fac := 1.0 / (2.0 * p);
rnd48(temp, rnd_rec.seed);
If (temp > p) then
fac := 1.0 / (2.0 * (1.0 - p));
End if;
rnd48(temp1, rnd_Rec.seed);
rnd_rec.rnd := (-rnd_rec.mean * fac) * log(temp1);
End if;
End if;
End if;
End hyperexponential;
Procedure gamma(rnd_rec: InOut rnd_rec_t) Is
----------------------------------------------------------------
--gamma : Generate a random number from the Gamma --
-- Distibution. --
-- --
--Reference : Pritsker, simulation and Slam II, 1984 --
-- Cheng, "The generation of gamma variables with --
-- nonintegeral shape parameters," applied --
-- statistics, vol 26, no 1, 1977, pp 71-75. --
-- --
--Used Fields : mean, std_dev, seed. --
-- --
--Changed Fields : rnd, seed. --
----------------------------------------------------------------
Variable accept : boolean;
Variable temp1, temp2, alpha, ra, a, b : real;
Variable rima, beta, x, y, w, c, v, z, r: real;
Begin
If rnd_rec.mean <= 0.0 then
Assert False
Report "The mean must be greater than zero"
Severity Warning;
rnd_rec.rnd := 1.0;
Elsif
rnd_rec.std_dev <= 0.0 then
Assert False
Report " The standard deviation must be greater than zero"
Severity Warning;
rnd_rec.rnd := rnd_rec.mean;
Else
alpha := (rnd_rec.mean * rnd_rec.mean) /
(rnd_rec.std_dev * rnd_rec.std_dev);
If alpha < 1.0 then
ra := 1.0 / alpha;
rima := 1.0 / (1.0 - alpha);
beta := rnd_rec.mean / alpha;
x := 1.0;
y := 1.0;
While x + y > 1.0 loop
rnd48(temp1, rnd_rec.seed);
rnd48(temp2, rnd_rec.seed);
x := temp1 ** ra;
y := temp2 ** rima;
End loop;
w := x / (x + y);
rnd48(temp1, rnd_rec.seed);
rnd_rec.rnd := w * (-log(temp1)) * beta;
Elsif
alpha <= 1.0 then
rnd48(temp1, rnd_rec.seed);
rnd_rec.rnd := -rnd_rec.mean * log(temp1);
Else
a := 1.0 / sqrt(2.0 * alpha - 1.0);
b := alpha - log(4.0);
c := Alpha + 1.0 / a;
accept := false;
While not accept loop
rnd48(temp1, rnd_rec.seed);
rnd48(temp2, rnd_rec.seed);
v := a * log(temp1 / (1.0 - temp1));
x := alpha * exp(v);
z := temp1 * temp1 * temp2;
r := b + c * v - x;
accept := ((r + 2.5040774 - 4.5 * z) >= 0.0
or (r > log(z)));
End loop;
rnd_rec.rnd := rnd_rec.std_dev * rnd_rec.std_dev
/ rnd_rec.mean * x;
End if;
End if;
End gamma;
Procedure erlang(rnd_rec: InOut rnd_rec_t) Is
----------------------------------------------------------------
--erlang : Generate a random number from a generalized --
-- Erlang distribution. --
-- --
--Used Fields : mean, std_dev, seed. --
-- --
--Changed Fields : rnd, seed. --
----------------------------------------------------------------
Variable coeff, fn, cfsq, p, tranx : real;
Variable temp1, temp2, stgmn, cumlog: real;
Variable n, i : integer;
Constant rnd_min : real := 1.0E-25;
Begin
If(rnd_rec.mean <= 0.0) or (rnd_rec.std_dev <= 0.0) then
Assert False
Report "The mean and the standard deviation must be positive"
Severity Warning;
rnd_rec.rnd := 1.0;
Return;
End if;
coeff := rnd_rec.std_dev / rnd_rec.mean;
If coeff > 1.0 then
Assert False Report
"The ratio of standard deviation to the mean must be less than one"
Severity Warning;
rnd_rec.rnd := rnd_rec.mean;
Return;
End if;
fn := 1.0 / (coeff * coeff);
n := Integer(floor(fn));
If fn /= real(n) then n := n + 1;
If n = 1 then n := 2;
End if;
fn := real(n);
cfsq := coeff * coeff;
p := ((2.0 * fn * cfsq) + fn - 2.0 - sqrt((fn * fn) + 4.0
- (4.0 * fn * cfsq))) / (2.0 * (cfsq + 1.0) * (fn - 1.0));
If (p < 0.0) or (p > 1.0) then
Assert False
Report "Improper mean and standard deviation"
Severity Warning;
rnd_rec.rnd := rnd_rec.mean;
Return;
End if;
Else
p := 0.0;
End if;
stgmn := rnd_Rec.mean / (fn - p * (fn - 1.0));
rnd48(temp1, rnd_rec.seed);
cumlog := 0.0;
If p = 0.0 then tranx := 0.0;
Else
rnd48(tranx, rnd_rec.seed);
End if;
If (tranx >= p) and (n > 1) then
For i in 2 to n loop
If temp1 <= rnd_min then
cumlog := cumlog + log(temp1);
temp1 := 1.0;
End if;
rnd48(temp2, rnd_rec.seed);
temp1 := temp1 * temp2;
End loop;
End if;
rnd_rec.rnd := (-stgmn) * (cumlog + log(temp1));
End erlang;
Procedure exponential(rnd_rec: Inout rnd_rec_t) Is
----------------------------------------------------------------
--exponential : Generate a random number exponentially --
-- distributed. --
-- --
--Used Fields : mean seed. --
-- --
--Changed Fields : rnd, seed. --
----------------------------------------------------------------
Variable temp : real;
Begin
If rnd_rec.mean <= 0.0 then
Assert False
Report "The mean must be a positive number"
Severity Warning;
End if;
temp := 0.0;
While temp = 0.0 loop
rnd48(temp, rnd_rec.seed);
End loop;
rnd_rec.rnd := -rnd_rec.mean * log(temp);
End exponential;
Procedure empirical_c(rnd_rec: InOut rnd_rec_t;
emp_data: In rnd_empirical_vector_t) Is
----------------------------------------------------------------
--empirical_c : Generate a random number from the cumulative --
-- distribution specified in emp_data. The last --
-- element of this array must have a probability --
-- of 1.0; the first element must have a --
-- probability of 0.0. --
-- --
--Used Fields : seed --
-- --
--Changed Fields : rnd, seed --
----------------------------------------------------------------
Variable u: Real;
Variable i: Integer;
Begin
If ((emp_data'Length = 0) Or (emp_data'Left > emp_data'Right)) Then
Assert False
Report "An ascending array of empirical data must be given"
Severity Warning;
rnd_rec.rnd := 0.0;
Elsif (emp_data(emp_data'Left).p /= 0.0) Then
Assert False
Report "The first data point's probability must be 0.0"
Severity Warning;
rnd_rec.rnd := 0.0;
Elsif (emp_data(emp_data'Right).p /= 1.0) Then
Assert False
Report "The last data point's probability must be 1.0"
Severity Warning;
rnd_rec.rnd := 0.0;
Else
rnd48(u, rnd_rec.seed);
i := emp_data'Left + 1;
While (emp_data(i).p < u) Loop
i := i + 1;
End Loop;
rnd_rec.rnd := (u - emp_data(i - 1).p) / (emp_data(i).p - emp_data(i - 1).p)
* (emp_data(i).x - emp_data(i - 1).x) + emp_data(i - 1).x;
End If;
End empirical_c;
-------------------------------------------------------------------------------
-- For some reason, this code causes some VHDL compilers to crash;
-- since it's rarely (if ever) used, we originally left it out.
-- I reinistated it as a test case for V-System. B. Hanna 11/18/93
-- It seems to work now, the collbrate is passing FILE Variables!!!
-- Procedure Rnd_Read_Record(file_f: Text; rnd_rec: inOut rnd_rec_t) Is
Procedure Rnd_Read_Record(rnd_rec: inOut rnd_rec_t) Is
-- ----------------------------------------------------------------
-- --rnd_read_rec : This procedure is called to read data from a --
-- -- file that the user calls. This procedure then --
-- -- calls the procedures rnd_reead_character and --
-- -- random_dist1. The data is then checked against--
-- -- the fields of the record rnd_rec_t and stored --
-- -- as those field. --
-- ----------------------------------------------------------------
Use Std.TextIO.All;
Type rnd_rec1_t is (mean, distribution, std_dev, trials, seed,
rnd, bound_l, bound_h, p_success);
Type rnd_lead_array is array(rnd_rec1_t) of string(1 to 21);
Constant rnd_lead : rnd_lead_array :=
( "mean ",
"distribution ",
"std_dev ",
"trials ",
"seed ",
"rnd ",
"bound_l ",
"bound_h ",
"p_success " );
Variable read_str, read_str1, read_str2, read_str3, read_str4 : string(1 to 21);
Variable l, lpr: Line;
Variable check : boolean;
File f : Text is out "rnd.out";
Variable i, int : integer;
Procedure rnd_read_character(l: inout line; read_str: inout string(1 to 21)) is
-- ----------------------------------------------------------------
-- --rnd_read_character : Called to read the data file and store --
-- -- the data. While reading the file, blank --
-- -- spaces are thrown out. --
-- ----------------------------------------------------------------
Variable i, j : integer;
Variable space: character;
Variable check: boolean;
variable ld : line;
Begin
If (l'Length = 0) Then
For j in 1 to 21 loop
read_str(j) := ' ';
End loop;
Else
read(l, space, check);
If check = true then
While (space = ' ') And (l'Length /= 0) loop
read(l, space);
End loop;
read_str(1) := space;
i := 2;
While i <= 21 loop
If (l'Length = 0) Then
Exit;
End If;
Read(l, read_str(i));
If read_str(i) = ' ' then
Exit;
End if;
i := i + 1;
End loop;
For j in i to 21 loop
read_str(j) := ' ';
End loop;
Else
Assert False
Report "Error in reading the data file"
Severity Warning;
End if;
End if;
End rnd_read_character;
--
Procedure random_dist1(rnd_rec : inout rnd_rec_t; l : inout line) Is
-- ----------------------------------------------------------------
-- --random_dist1 : The data read is checked against all of the --
-- -- distributions and are stored in the record --
-- -- field of rnd_rec.distribution --
-- ----------------------------------------------------------------
Type Dist_str_array is array(rnd_distribution_t) of string(1 to 21);
Constant dist_str : Dist_str_array :=
( "rnd_constant ",
"rnd_uniform_d ",
"rnd_poisson ",
"rnd_binomial ",
"rnd_geometric ",
"rnd_empirical_d ",
"rnd_uniform_c ",
"rnd_normal ",
"rnd_triangular ",
"rnd_exponential ",
"rnd_hyperexponential ",
"rnd_gamma ",
"rnd_erlang ",
"rnd_empirical_c " );
Variable dist : String(1 to 21);
variable lpr : line;
Variable check: boolean;
variable i : rnd_distribution_t;
Begin
rnd_read_character(l, dist);
For i in rnd_distribution_t loop
If dist = (dist_str(i)) then
rnd_rec.distribution := i;
exit;
End if;
Assert i /= rnd_distribution_t'High
Report "Unrecognized distribution found"
Severity Warning;
End loop;
End random_dist1;
Begin
While not Endfile(file_f) loop
readline(file_f, l);
rnd_read_character(l, read_str);
For i in rnd_rec1_t loop
If read_str = rnd_lead(i) then
If rnd_lead(i) = rnd_lead(seed) then
read(l, int);
rnd_rec.seed(0) := int;
read(l, int);
rnd_rec.seed(1) := int;
read(l, int);
rnd_rec.seed(2) := int;
read(l, int);
rnd_rec.seed(3) := int;
Elsif (rnd_lead(i) = rnd_lead(distribution)) then
random_dist1(rnd_rec, l);
Else
If (read_str = rnd_lead(mean)) Then
read(l, rnd_rec.mean);
Elsif (read_str = rnd_lead(std_dev)) Then
read(l, rnd_rec.std_dev);
Elsif (read_str = rnd_lead(trials)) Then
read(l, rnd_rec.trials);
Elsif (read_str = rnd_lead(rnd)) Then
read(l, rnd_rec.rnd);
Elsif (read_str = rnd_lead(bound_l)) Then
read(l, rnd_rec.bound_l);
Elsif (read_str = rnd_lead(bound_h)) Then
read(l, rnd_rec.bound_h);
Elsif (read_str = rnd_lead(p_success)) Then
read(l, rnd_rec.p_success);
Else
Assert False
report "Invalid field name found"
Severity warning;
End If;
End if;
End if;
End loop;
End loop;
End Rnd_Read_Record;
-- -------------------------------------------------------------------------
Procedure Rnd_Random (rnd_rec: InOut rnd_rec_t; empirical_data:
In rnd_empirical_vector_t) Is
-- V-Systems Does Not Like Initial Values
-- := (0 => (0.0, 0.0),1 => (1.0, 1.0))
Begin
Case rnd_rec.distribution Is
When rnd_constant => rnd_rec.rnd := rnd_rec.mean;
When rnd_uniform_d => uniform_d(rnd_rec);
When rnd_poisson => poisson(rnd_rec);
When rnd_binomial => binomial(rnd_rec);
When rnd_geometric => geometric(rnd_rec);
When rnd_uniform_c => uniform_c(rnd_rec);
When rnd_normal => normal(rnd_rec);
When rnd_triangular => triangular(rnd_rec);
When rnd_hyperexponential => hyperexponential(rnd_rec);
When rnd_gamma => gamma(rnd_rec);
When rnd_erlang => erlang(rnd_rec);
When rnd_empirical_d => empirical_d(rnd_rec, empirical_data);
When rnd_empirical_c => empirical_c(rnd_rec, empirical_data);
When rnd_exponential => exponential(rnd_rec);
End Case;
End Rnd_Random;
End rnd2;
----------------------------- Cut Here --------------------------------------
The following Is File "XYZ.IN" !!!
Ignore it if you do not call Rnd_read_record ........
1.10185E-1 4.550113 1.438351 2.684018 1.590097 2.695197 4.359265 4.191798
4.954218 1.311093 2.333256 1.181042 2.948619 2.086714 8.816002E-1 2.973476
2.203053 4.05536
2.929811 3.159666E-1 4.214381 2.610496 2.989147E-1 1.388424 2.746835
2.475202 2.294238 2.224455 6.181239E-1 2.641986 4.204303 2.617279 2.558914
2.772824 2.828106 2.867225
2.892393E-1 2.331572 4.872393 2.730826 1.165035 1.911058 2.42206 1.770379
2.974291 2.599989 1.861426E-1 2.86104 2.453901 1.374658 2.232845 1.706726
2.21552 2.137208
1.316083 2.279989 3.401187E-1 2.170335 1.417249 4.985046 2.505599 1.097503
1.417657 2.098764 4.84458 4.569949 7.809391E-1 2.283342 2.838585 2.556404
2.494201 2.293147E-1
2.449824 2.159662 2.213958 4.300379 2.841128 2.486782 1.846618 6.957985E-1
2.480724 4.250188 2.722524 2.236344 4.217936 4.003117 2.280461 2.630268
2.550001E-1 1.228973
2.432141 1.665156 3.641811E-1 2.181681 4.620957 2.044269 2.570266 2.833069
2.439487 1.008646 2.252072 1.581691E-1 2.589334 2.094404 2.81995 4.517756
2.766231 2.353488
2.288909 2.228878 2.615364 2.520471 1.887305 4.244486E-1 2.072604 4.56562
4.248774 1.713845E-2 2.681152 4.703273 2.404612 2.180274 2.51785 6.619296E-2
4.843644 4.633335
2.853335 2.025907 7.740762E-2 2.00919 1.547245 2.984337 2.518233 4.878979
2.003666 2.374277 2.700533 2.308045 8.026191E-1 2.757151 2.983902 2.251004
2.09121 2.258416
1.619651 8.716932E-1 4.936877 1.789543 4.941626 2.993652E-1 4.72275E-1
2.005302 1.472767 2.51507 2.734172 2.894203 4.549603 4.099712 1.781163
9.400108E-1 1.349926 2.123261
1.483356 2.188843 2.597879 1.921498 8.743594E-1 6.380821E-1 4.496956
4.625475 1.910821 4.574652 4.953065 4.031111 2.283712 2.297944 1.563671
1.75171 1.317444 4.533722
2.960496 4.680359 4.389577 2.766362 2.322055 2.051096 7.56996E-1 1.719276
2.033643 1.477548 1.606924 4.49642 2.797199 4.435889 1.374035 2.682927
2.632812 2.65763
2.475743 2.661489 6.052993E-1 2.85374 2.575287 4.731117 4.70633 1.719285
2.616551 2.27777 2.463141 1.691975 2.954899 1.582471 2.948478 4.091076
4.227309 2.716923
2.554578 2.023056 2.860676 2.053273 2.266132 1.925981 2.739061 2.552348
4.569175 4.985281 2.525038 2.709172 1.561703E-1 6.001482E-1 2.838769 1.175169
2.963118E-1 2.58157
RN: 2.342977 2.546691 1.443457 4.682113 2.581236 1.553273 2.411718 4.010745
1.661809 2.145081 4.804929 2.362327 4.491977 1.927268 2.650728 4.511765
4.079527 4.546187
<div align="center"><br /><script type="text/javascript"><!--
google_ad_client = "pub-7293844627074885";
//468x60, Created at 07. 11. 25
google_ad_slot = "8619794253";
google_ad_width = 468;
google_ad_height = 60;
//--></script>
<script type="text/javascript" src="http://pagead2.googlesyndication.com/pagead/show_ads.js">
</script><br /> </div>