Subversion Repositories rng_lib

----------------------------------------------------------------------
----                                                              ----
---- Rand number generator library.                               ----
----                                                              ----
---- This file is part of the Random Number Generator project     ----
---- http://www.opencores.org/cores/rng_lib/                      ----
----                                                              ----
---- Description                                                  ----
---- This library has function for generation random numbers with ----
---- the following distributions:                                 ----
---- - Uniform (continous)                                        ----
---- - Exponential (continous)                                    ----
---- - Gaussian (continous)                                       ----
----                                                              ----
---- Random numbers are produced with a combination of 3          ----
---- Tausworthe generators which gives very good statistical      ----
---- properties.                                                  ----
----                                                              ----
---- NOTE! These functions will NOT synthesize. They are for test ----
----       bench use only!                                        ----
----                                                              ----
---- To Do:                                                       ----
---- -                                                            ----
----                                                              ----
---- Author(s):                                                   ----
---- - Geir Drange, gedra@opencores.org                           ----
----                                                              ----
----------------------------------------------------------------------
----                                                              ----
---- Copyright (C) 2004 Authors and OPENCORES.ORG                 ----
----                                                              ----
---- This source file may be used and distributed without         ----
---- restriction provided that this copyright statement is not    ----
---- removed from the file and that any derivative work contains  ----
---- the original copyright notice and the associated disclaimer. ----
----                                                              ----
---- This source file is free software; you can redistribute it   ----
---- and/or modify it under the terms of the GNU General          ----
---- Public License as published by the Free Software Foundation; ----
---- either version 2.0 of the License, or (at your option) any   ----
---- later version.                                               ----
----                                                              ----
---- This source 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 source; if not, download it   ----
---- from http://www.gnu.org/licenses/gpl.txt                     ----
----                                                              ----
----------------------------------------------------------------------
--
-- CVS Revision History
--
-- $Log: not supported by cvs2svn $
--
--
 
library ieee;
use ieee.std_logic_1164.all;
use ieee.numeric_std.all;
use std.textio.all;
use work.math_lib.all;
 
package rng_lib is
 
  type distribution is (UNIFORM, GAUSSIAN, EXPONENTIAL);
  type rand_var is record  -- random variable record
                     rnd: real;         -- random number
                     rnd_v: unsigned(31 downto 0); -- random number vector
                     dist: distribution;      -- distribution type
                     y, z: real;        -- distribution parameters
                     s1, s2, s3: unsigned(31 downto 0);  -- seeds
                     mask1, mask2, mask3: unsigned(31 downto 0);
                     shft1, shft2, shft3: natural;
                   end record;
 
  function rand (rnd: rand_var) return rand_var;
  function init_uniform(constant a, b, c: natural;
                        constant lo, hi: real) return rand_var;
  function init_gaussian(constant a, b, c: natural;
                         constant mean, stdev: real) return rand_var;
  function init_exponential(constant a, b, c: natural;
                        constant mean: real) return rand_var;
 
  constant q1: natural:=13;
  constant q2: natural:=2;
  constant q3: natural:=3;
  constant p1: natural:=12;
  constant p2: natural:=4;
  constant p3: natural:=17;
 
end rng_lib;
 
package body rng_lib is
 
-- Function to convert 32bit unsigned vector to real
-- Integers only go to 2**31 (VHDL'87), so do it clever
  function unsigned_2_real (constant a: unsigned(31 downto 0)) return real is
    variable r: real;
  begin
    r := 2.0*real(to_integer(a(31 downto 1)));
    if a(0) = '1' then
      r := r + 1.0;
    end if;
    return(r);
  end unsigned_2_real;
 
-- Generate random number using a combination of 3 tausworthe generators 
-- Source: Pierre L'Ecuyer, "Maximally Equidistributed Combined Tausworthe
-- Generators". Mathematics of Computation, vol.65, no.213(1996), pp203--213.
  function rng (rnd: rand_var) return rand_var is
    variable new_rnd : rand_var;
    variable b : unsigned(31 downto 0);
  begin
    new_rnd :=rnd;
    b := ((new_rnd.s1 sll q1) xor new_rnd.s1) srl new_rnd.shft1;
    new_rnd.s1 := ((new_rnd.s1 and new_rnd.mask1) sll p1) xor b;
    b := ((new_rnd.s2 sll q2) xor new_rnd.s2) srl new_rnd.shft2;
    new_rnd.s2 := ((new_rnd.s2 and new_rnd.mask2) sll p2) xor b;
    b := ((new_rnd.s3 sll q3) xor new_rnd.s3) srl new_rnd.shft3;
    new_rnd.s3 := ((new_rnd.s3 and new_rnd.mask3) sll p3) xor b;
    new_rnd.rnd_v := new_rnd.s1 xor new_rnd.s2 xor new_rnd.s3;
    -- normalize to range [0,1)
    new_rnd.rnd := unsigned_2_real(new_rnd.rnd_v) / 65536.0;
    new_rnd.rnd := new_rnd.rnd / 65536.0;
    return (new_rnd);
  end rng;
 
-- rand function generates a random variable with different distributions
  function rand (rnd: rand_var) return rand_var is 
    variable rnd_out : rand_var;
    variable x,y,z : real;
    variable t: real := 0.0;
  begin
    case rnd.dist is
      -- Uniform distribution
      when UNIFORM =>
        rnd_out := rng(rnd);
        rnd_out.rnd := rnd.y + (rnd_out.rnd * (rnd.z - rnd.y));
      -- Gaussian distribution
      when GAUSSIAN =>                  -- Box-Mueller method
        z := 2.0;
        rnd_out := rnd;
        while z > 1.0 or z = 0.0 loop
          -- choose x,y in uniform square (-1,-1) to (+1,+1)
          rnd_out := rng(rnd_out);
          x := -1.0 + 2.0 * rnd_out.rnd;
          rnd_out := rng(rnd_out);
          y := -1.0 + 2.0 * rnd_out.rnd;
          z := (x * x) + (y * y);
        end loop;
        -- Box-Mueller transform
        rnd_out.rnd := rnd_out.y + rnd_out.z * y * sqrt(-2.0 * log(z)/z);
      -- Exponential distribution
      when EXPONENTIAL =>
        rnd_out := rng(rnd);
        rnd_out.rnd := -rnd_out.y * log(1.0 - rnd_out.rnd);
      when others =>
        report "rand() function encountered an error!"
          severity failure;
    end case;
    return (rnd_out);
  end rand;
 
-- Initialize seeds, used by all init_ functions
  function gen_seed (constant a, b, c: natural) return rand_var is
    variable seeded : rand_var;
    variable x : unsigned(31 downto 0):= "11111111111111111111111111111111";
    constant k1: natural:=31;
    constant k2: natural:=29;
    constant k3: natural:=28;
  begin
    seeded.shft1 := k1-p1;
    seeded.shft2 := k2-p2;
    seeded.shft3 := k3-p3;
    seeded.mask1 := x sll (32-k1);
    seeded.mask2 := x sll (32-k2);
    seeded.mask3 := x sll (32-k3);
    seeded.s1 := to_unsigned(390451501, 32); 
    seeded.s2 := to_unsigned(613566701, 32);
    seeded.s3 := to_unsigned(858993401, 32); 
    if to_unsigned(a, 32) > (to_unsigned(1,32) sll (32-k1)) then
      seeded.s1 := to_unsigned(a, 32);
    end if;
    if to_unsigned(b, 32) > (to_unsigned(1,32) sll (32-k2)) then
      seeded.s2 := to_unsigned(b, 32);
    end if;
    if to_unsigned(c, 32) > (to_unsigned(1,32) sll (32-k3)) then
      seeded.s3 := to_unsigned(c, 32);
    end if;
    return(seeded);
  end gen_seed;
 
-- Uniform distribution random variable initialization
-- a,b,c are seeds
-- lo,hi is the range for the uniform distribution
  function init_uniform(constant a, b, c: natural;
                        constant lo, hi: real) return rand_var is
    variable rnd, rout : rand_var;
  begin
    if lo >= hi then
      report "Uniform parameter error: 'hi' must be > 'lo'!"
      severity failure;
    end if;
    rnd := gen_seed(a, b, c);
    rnd.dist := UNIFORM;
    rnd.y := lo;
    rnd.z := hi;
    rout := rand(rnd);
    return(rout);
  end init_uniform;
 
-- Gaussian distribution random variable initialization
-- a,b,c are seeds
-- mean,stdev is mean and standard deviation
  function init_gaussian(constant a, b, c: natural;
                        constant mean, stdev: real) return rand_var is
    variable rnd, rout : rand_var;
  begin
    if stdev = 0.0 then
      report "Gaussian parameter error: 'stdev' must be non-zero!"
      severity failure;
    end if;
    rnd := gen_seed(a, b, c);
    rnd.dist := GAUSSIAN;
    rnd.y := mean;
    rnd.z := stdev;
    rout := rand(rnd);
    return(rout);
  end init_gaussian;
 
-- Exponential distribution random variable initialization
-- a,b,c are seeds
-- mean: mean value
  function init_exponential(constant a, b, c: natural;
                            constant mean: real) return rand_var is
    variable rnd, rout : rand_var;
  begin
    if mean <= 0.0 then
      report "Exponential parameter error: 'mean' must be > 0!"
      severity failure;
    end if;
    rnd := gen_seed(a, b, c);
    rnd.dist := EXPONENTIAL;
    rnd.y := mean;
    rout := rand(rnd);
    return(rout);
  end init_exponential;  
 
end rng_lib;