Subversion Repositories openrisc

------------------------------------------------------------------------------

--                                                                          --

--                         GNAT RUN-TIME COMPONENTS                         --

--                                                                          --

--                S Y S T E M . R A N D O M _ N U M B E R S                 --

--                                                                          --

--                                 B o d y                                  --

--                                                                          --

--          Copyright (C) 2007-2012, Free Software Foundation, Inc.         --

--                                                                          --

-- GNAT is free software;  you can  redistribute it  and/or modify it under --

-- terms of the  GNU General Public License as published  by the Free Soft- --

-- ware  Foundation;  either version 3,  or (at your option) any later ver- --

-- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --

-- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --

-- or FITNESS FOR A PARTICULAR PURPOSE.                                     --

--                                                                          --

-- As a special exception under Section 7 of GPL version 3, you are granted --

-- additional permissions described in the GCC Runtime Library Exception,   --

-- version 3.1, as published by the Free Software Foundation.               --

--                                                                          --

-- You should have received a copy of the GNU General Public License and    --

-- a copy of the GCC Runtime Library Exception along with this program;     --

-- see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see    --

-- <http://www.gnu.org/licenses/>.                                          --

--                                                                          --

-- GNAT was originally developed  by the GNAT team at  New York University. --

-- Extensive contributions were provided by Ada Core Technologies Inc.      --

--                                                                          --

------------------------------------------------------------------------------

------------------------------------------------------------------------------

--                                                                          --

-- The implementation here is derived from a C-program for MT19937, with    --

-- initialization improved 2002/1/26. As required, the following notice is  --

-- copied from the original program.                                        --

--                                                                          --

-- Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,        --

-- All rights reserved.                                                     --

--                                                                          --

-- Redistribution and use in source and binary forms, with or without       --

-- modification, are permitted provided that the following conditions       --

-- are met:                                                                 --

--                                                                          --

--   1. Redistributions of source code must retain the above copyright      --

--      notice, this list of conditions and the following disclaimer.       --

--                                                                          --

--   2. Redistributions in binary form must reproduce the above copyright   --

--      notice, this list of conditions and the following disclaimer in the --

--      documentation and/or other materials provided with the distribution.--

--                                                                          --

--   3. The names of its contributors may not be used to endorse or promote --

--      products derived from this software without specific prior written  --

--      permission.                                                         --

--                                                                          --

-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS      --

-- "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT        --

-- LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR    --

-- A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT    --

-- OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,    --

-- SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED --

-- TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR   --

-- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF   --

-- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING     --

-- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS       --

-- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.             --

--                                                                          --

------------------------------------------------------------------------------

------------------------------------------------------------------------------

--                                                                          --

-- This is an implementation of the Mersenne Twister, twisted generalized   --

-- feedback shift register of rational normal form, with state-bit          --

-- reflection and tempering. This version generates 32-bit integers with a  --

-- period of 2**19937 - 1 (a Mersenne prime, hence the name). For           --

-- applications requiring more than 32 bits (up to 64), we concatenate two  --

-- 32-bit numbers.                                                          --

--                                                                          --

-- See http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html for         --

-- details.                                                                 --

--                                                                          --

-- In contrast to the original code, we do not generate random numbers in   --

-- batches of N. Measurement seems to show this has very little if any      --

-- effect on performance, and it may be marginally better for real-time     --

-- applications with hard deadlines.                                        --

--                                                                          --

------------------------------------------------------------------------------

with Ada.Unchecked_Conversion;

with System.Random_Seed;

with Interfaces; use Interfaces;

use Ada;

package body System.Random_Numbers is

   Image_Numeral_Length : constant := Max_Image_Width / N;

   subtype Image_String is String (1 .. Max_Image_Width);

   ----------------------------

   -- Algorithmic Parameters --

   ----------------------------

   Lower_Mask : constant := 2**31-1;

   Upper_Mask : constant := 2**31;

   Matrix_A   : constant array (State_Val range 0 .. 1) of State_Val

     := (0, 16#9908b0df#);

   --  The twist transformation is represented by a matrix of the form

   --

   --               [  0    I(31) ]

   --               [    _a       ]

   --

   --  where 0 is a 31x31 block of 0s, I(31) is the 31x31 identity matrix and

   --  _a is a particular bit row-vector, represented here by a 32-bit integer.

   --  If integer x represents a row vector of bits (with x(0), the units bit,

   --  last), then

   --           x * A = [0 x(31..1)] xor Matrix_A(x(0)).

   U      : constant := 11;

   S      : constant := 7;

   B_Mask : constant := 16#9d2c5680#;

   T      : constant := 15;

   C_Mask : constant := 16#efc60000#;

   L      : constant := 18;

   --  The tempering shifts and bit masks, in the order applied

   Seed0 : constant := 5489;

   --  Default seed, used to initialize the state vector when Reset not called

   Seed1 : constant := 19650218;

   --  Seed used to initialize the state vector when calling Reset with an

   --  initialization vector.

   Mult0 : constant := 1812433253;

   --  Multiplier for a modified linear congruential generator used to

   --  initialize the state vector when calling Reset with a single integer

   --  seed.

   Mult1 : constant := 1664525;

   Mult2 : constant := 1566083941;

   --  Multipliers for two modified linear congruential generators used to

   --  initialize the state vector when calling Reset with an initialization

   --  vector.

   -----------------------

   -- Local Subprograms --

   -----------------------

   procedure Init (Gen : Generator; Initiator : Unsigned_32);

   --  Perform a default initialization of the state of Gen. The resulting

   --  state is identical for identical values of Initiator.

   procedure Insert_Image

     (S     : in out Image_String;

      Index : Integer;

      V     : State_Val);

   --  Insert image of V into S, in the Index'th 11-character substring

   function Extract_Value (S : String; Index : Integer) return State_Val;

   --  Treat S as a sequence of 11-character decimal numerals and return

   --  the result of converting numeral #Index (numbering from 0)

   function To_Unsigned is

     new Unchecked_Conversion (Integer_32, Unsigned_32);

   function To_Unsigned is

     new Unchecked_Conversion (Integer_64, Unsigned_64);

   ------------

   -- Random --

   ------------

   function Random (Gen : Generator) return Unsigned_32 is

      G : Generator renames Gen.Writable.Self.all;

      Y : State_Val;

      I : Integer;      --  should avoid use of identifier I ???

   begin

      I := G.I;

      if I < N - M then

         Y := (G.S (I) and Upper_Mask) or (G.S (I + 1) and Lower_Mask);

         Y := G.S (I + M) xor Shift_Right (Y, 1) xor Matrix_A (Y and 1);

         I := I + 1;

      elsif I < N - 1 then

         Y := (G.S (I) and Upper_Mask) or (G.S (I + 1) and Lower_Mask);

         Y := G.S (I + (M - N))

                xor Shift_Right (Y, 1)

                xor Matrix_A (Y and 1);

         I := I + 1;

      elsif I = N - 1 then

         Y := (G.S (I) and Upper_Mask) or (G.S (0) and Lower_Mask);

         Y := G.S (M - 1) xor Shift_Right (Y, 1) xor Matrix_A (Y and 1);

         I := 0;

      else

         Init (G, Seed0);

         return Random (Gen);

      end if;

      G.S (G.I) := Y;

      G.I := I;

      Y := Y xor Shift_Right (Y, U);

      Y := Y xor (Shift_Left (Y, S)  and B_Mask);

      Y := Y xor (Shift_Left (Y, T) and C_Mask);

      Y := Y xor Shift_Right (Y, L);

      return Y;

   end Random;

   generic

      type Unsigned is mod <>;

      type Real is digits <>;

      with function Random (G : Generator) return Unsigned is <>;

   function Random_Float_Template (Gen : Generator) return Real;

   pragma Inline (Random_Float_Template);

   --  Template for a random-number generator implementation that delivers

   --  values of type Real in the range [0 .. 1], using values from Gen,

   --  assuming that Unsigned is large enough to hold the bits of a mantissa

   --  for type Real.

   ---------------------------

   -- Random_Float_Template --

   ---------------------------

   function Random_Float_Template (Gen : Generator) return Real is

      pragma Compile_Time_Error

        (Unsigned'Last <= 2**(Real'Machine_Mantissa - 1),

         "insufficiently large modular type used to hold mantissa");

   begin

      --  This code generates random floating-point numbers from unsigned

      --  integers. Assuming that Real'Machine_Radix = 2, it can deliver all

      --  machine values of type Real (as implied by Real'Machine_Mantissa and

      --  Real'Machine_Emin), which is not true of the standard method (to

      --  which we fall back for non-binary radix): computing Real(<random

      --  integer>) / (<max random integer>+1). To do so, we first extract an

      --  (M-1)-bit significand (where M is Real'Machine_Mantissa), and then

      --  decide on a normalized exponent by repeated coin flips, decrementing

      --  from 0 as long as we flip heads (1 bits). This process yields the

      --  proper geometric distribution for the exponent: in a uniformly

      --  distributed set of floating-point numbers, 1/2 of them will be in

      --  (0.5, 1], 1/4 will be in (0.25, 0.5], and so forth. It makes a

      --  further adjustment at binade boundaries (see comments below) to give

      --  the effect of selecting a uniformly distributed real deviate in

      --  [0..1] and then rounding to the nearest representable floating-point

      --  number.  The algorithm attempts to be stingy with random integers. In

      --  the worst case, it can consume roughly -Real'Machine_Emin/32 32-bit

      --  integers, but this case occurs with probability around

      --  2**Machine_Emin, and the expected number of calls to integer-valued

      --  Random is 1.  For another discussion of the issues addressed by this

      --  process, see Allen Downey's unpublished paper at

      --  http://allendowney.com/research/rand/downey07randfloat.pdf.

      if Real'Machine_Radix /= 2 then

         return Real'Machine

           (Real (Unsigned'(Random (Gen))) * 2.0**(-Unsigned'Size));

      else

         declare

            type Bit_Count is range 0 .. 4;

            subtype T is Real'Base;

            Trailing_Ones : constant array (Unsigned_32 range 0 .. 15)

              of Bit_Count :=

                  (2#00000# => 0, 2#00001# => 1, 2#00010# => 0, 2#00011# => 2,

                   2#00100# => 0, 2#00101# => 1, 2#00110# => 0, 2#00111# => 3,

                   2#01000# => 0, 2#01001# => 1, 2#01010# => 0, 2#01011# => 2,

                   2#01100# => 0, 2#01101# => 1, 2#01110# => 0, 2#01111# => 4);

            Pow_Tab : constant array (Bit_Count range 0 .. 3) of Real

              := (0 => 2.0**(0 - T'Machine_Mantissa),

                  1 => 2.0**(-1 - T'Machine_Mantissa),

                  2 => 2.0**(-2 - T'Machine_Mantissa),

                  3 => 2.0**(-3 - T'Machine_Mantissa));

            Extra_Bits : constant Natural :=

                         (Unsigned'Size - T'Machine_Mantissa + 1);

            --  Random bits left over after selecting mantissa

            Mantissa : Unsigned;

            X      : Real;            --  Scaled mantissa

            R      : Unsigned_32;     --  Supply of random bits

            R_Bits : Natural;         --  Number of bits left in R

            K      : Bit_Count;       --  Next decrement to exponent

         begin

            Mantissa := Random (Gen) / 2**Extra_Bits;

            R := Unsigned_32 (Mantissa mod 2**Extra_Bits);

            R_Bits := Extra_Bits;

            X := Real (2**(T'Machine_Mantissa - 1) + Mantissa); -- Exact

            if Extra_Bits < 4 and then R < 2 ** Extra_Bits - 1 then

               --  We got lucky and got a zero in our few extra bits

               K := Trailing_Ones (R);

            else

               Find_Zero : loop

                  --  R has R_Bits unprocessed random bits, a multiple of 4.

                  --  X needs to be halved for each trailing one bit. The

                  --  process stops as soon as a 0 bit is found. If R_Bits

                  --  becomes zero, reload R.

                  --  Process 4 bits at a time for speed: the two iterations

                  --  on average with three tests each was still too slow,

                  --  probably because the branches are not predictable.

                  --  This loop now will only execute once 94% of the cases,

                  --  doing more bits at a time will not help.

                  while R_Bits >= 4 loop

                     K := Trailing_Ones (R mod 16);

                     exit Find_Zero when K < 4; -- Exits 94% of the time

                     R_Bits := R_Bits - 4;

                     X := X / 16.0;

                     R := R / 16;

                  end loop;

                  --  Do not allow us to loop endlessly even in the (very

                  --  unlikely) case that Random (Gen) keeps yielding all ones.

                  exit Find_Zero when X = 0.0;

                  R := Random (Gen);

                  R_Bits := 32;

               end loop Find_Zero;

            end if;

            --  K has the count of trailing ones not reflected yet in X. The

            --  following multiplication takes care of that, as well as the

            --  correction to move the radix point to the left of the mantissa.

            --  Doing it at the end avoids repeated rounding errors in the

            --  exceedingly unlikely case of ever having a subnormal result.

            X := X * Pow_Tab (K);

            --  The smallest value in each binade is rounded to by 0.75 of

            --  the span of real numbers as its next larger neighbor, and

            --  1.0 is rounded to by half of the span of real numbers as its

            --  next smaller neighbor. To account for this, when we encounter

            --  the smallest number in a binade, we substitute the smallest

            --  value in the next larger binade with probability 1/2.

            if Mantissa = 0 and then Unsigned_32'(Random (Gen)) mod 2 = 0 then

               X := 2.0 * X;

            end if;

            return X;

         end;

      end if;

   end Random_Float_Template;

   ------------

   -- Random --

   ------------

   function Random (Gen : Generator) return Float is

      function F is new Random_Float_Template (Unsigned_32, Float);

   begin

      return F (Gen);

   end Random;

   function Random (Gen : Generator) return Long_Float is

      function F is new Random_Float_Template (Unsigned_64, Long_Float);

   begin

      return F (Gen);

   end Random;

   function Random (Gen : Generator) return Unsigned_64 is

   begin

      return Shift_Left (Unsigned_64 (Unsigned_32'(Random (Gen))), 32)

        or Unsigned_64 (Unsigned_32'(Random (Gen)));

   end Random;

   ---------------------

   -- Random_Discrete --

   ---------------------

   function Random_Discrete

     (Gen : Generator;

      Min : Result_Subtype := Default_Min;

      Max : Result_Subtype := Result_Subtype'Last) return Result_Subtype

   is

   begin

      if Max = Min then

         return Max;

      elsif Max < Min then

         raise Constraint_Error;

      elsif Result_Subtype'Base'Size > 32 then

         declare

            --  In the 64-bit case, we have to be careful, since not all 64-bit

            --  unsigned values are representable in GNAT's root_integer type.

            --  Ignore different-size warnings here since GNAT's handling

            --  is correct.

            pragma Warnings ("Z");  -- better to use msg string! ???

            function Conv_To_Unsigned is

               new Unchecked_Conversion (Result_Subtype'Base, Unsigned_64);

            function Conv_To_Result is

               new Unchecked_Conversion (Unsigned_64, Result_Subtype'Base);

            pragma Warnings ("z");

            N : constant Unsigned_64 :=

                  Conv_To_Unsigned (Max) - Conv_To_Unsigned (Min) + 1;

            X, Slop : Unsigned_64;

         begin

            if N = 0 then

               return Conv_To_Result (Conv_To_Unsigned (Min) + Random (Gen));

            else

               Slop := Unsigned_64'Last rem N + 1;

               loop

                  X := Random (Gen);

                  exit when Slop = N or else X <= Unsigned_64'Last - Slop;

               end loop;

               return Conv_To_Result (Conv_To_Unsigned (Min) + X rem N);

            end if;

         end;

      elsif Result_Subtype'Pos (Max) - Result_Subtype'Pos (Min) =

                                                         2 ** 32 - 1

      then

         return Result_Subtype'Val

           (Result_Subtype'Pos (Min) + Unsigned_32'Pos (Random (Gen)));

      else

         declare

            N    : constant Unsigned_32 :=

                     Unsigned_32 (Result_Subtype'Pos (Max) -

                                    Result_Subtype'Pos (Min) + 1);

            Slop : constant Unsigned_32 := Unsigned_32'Last rem N + 1;

            X    : Unsigned_32;

         begin

            loop

               X := Random (Gen);

               exit when Slop = N or else X <= Unsigned_32'Last - Slop;

            end loop;

            return

              Result_Subtype'Val

                (Result_Subtype'Pos (Min) + Unsigned_32'Pos (X rem N));

         end;

      end if;

   end Random_Discrete;

   ------------------

   -- Random_Float --

   ------------------

   function Random_Float (Gen : Generator) return Result_Subtype is

   begin

      if Result_Subtype'Base'Digits > Float'Digits then

         return Result_Subtype'Machine (Result_Subtype

                                         (Long_Float'(Random (Gen))));

      else

         return Result_Subtype'Machine (Result_Subtype

                                         (Float'(Random (Gen))));

      end if;

   end Random_Float;

   -----------

   -- Reset --

   -----------

   procedure Reset (Gen : Generator) is

   begin

      Init (Gen, Unsigned_32'Mod (Random_Seed.Get_Seed));

   end Reset;

   procedure Reset (Gen : Generator; Initiator : Integer_32) is

   begin

      Init (Gen, To_Unsigned (Initiator));

   end Reset;

   procedure Reset (Gen : Generator; Initiator : Unsigned_32) is

   begin

      Init (Gen, Initiator);

   end Reset;

   procedure Reset (Gen : Generator; Initiator : Integer) is

   begin

      pragma Warnings (Off, "condition is always *");

      --  This is probably an unnecessary precaution against future change, but

      --  since the test is a static expression, no extra code is involved.

      if Integer'Size <= 32 then

         Init (Gen, To_Unsigned (Integer_32 (Initiator)));

      else

         declare

            Initiator1 : constant Unsigned_64 :=

                           To_Unsigned (Integer_64 (Initiator));

            Init0      : constant Unsigned_32 :=

                           Unsigned_32 (Initiator1 mod 2 ** 32);

            Init1      : constant Unsigned_32 :=

                           Unsigned_32 (Shift_Right (Initiator1, 32));

         begin

            Reset (Gen, Initialization_Vector'(Init0, Init1));

         end;

      end if;

      pragma Warnings (On, "condition is always *");

   end Reset;

   procedure Reset (Gen : Generator; Initiator : Initialization_Vector) is

      G    : Generator renames Gen.Writable.Self.all;

      I, J : Integer;

   begin

      Init (G, Seed1);

      I := 1;

      J := 0;

      if Initiator'Length > 0 then

         for K in reverse 1 .. Integer'Max (N, Initiator'Length) loop

            G.S (I) :=

              (G.S (I) xor ((G.S (I - 1)

                               xor Shift_Right (G.S (I - 1), 30)) * Mult1))

              + Initiator (J + Initiator'First) + Unsigned_32 (J);

            I := I + 1;

            J := J + 1;

            if I >= N then

               G.S (0) := G.S (N - 1);

               I := 1;

            end if;

            if J >= Initiator'Length then

               J := 0;

            end if;

         end loop;

      end if;

      for K in reverse 1 .. N - 1 loop

         G.S (I) :=

           (G.S (I) xor ((G.S (I - 1)

                            xor Shift_Right (G.S (I - 1), 30)) * Mult2))

           - Unsigned_32 (I);

         I := I + 1;

         if I >= N then

            G.S (0) := G.S (N - 1);

            I := 1;

         end if;

      end loop;

      G.S (0) := Upper_Mask;

   end Reset;

   procedure Reset (Gen : Generator; From_State : Generator) is

      G : Generator renames Gen.Writable.Self.all;

   begin

      G.S := From_State.S;

      G.I := From_State.I;

   end Reset;

   procedure Reset (Gen : Generator; From_State : State) is

      G : Generator renames Gen.Writable.Self.all;

   begin

      G.I := 0;

      G.S := From_State;

   end Reset;

   procedure Reset (Gen : Generator; From_Image : String) is

      G : Generator renames Gen.Writable.Self.all;

   begin

      G.I := 0;

      for J in 0 .. N - 1 loop

         G.S (J) := Extract_Value (From_Image, J);

      end loop;

   end Reset;

   ----------

   -- Save --

   ----------

   procedure Save (Gen : Generator; To_State : out State) is

      Gen2 : Generator;

   begin

      if Gen.I = N then

         Init (Gen2, 5489);

         To_State := Gen2.S;

      else

         To_State (0 .. N - 1 - Gen.I) := Gen.S (Gen.I .. N - 1);

         To_State (N - Gen.I .. N - 1) := Gen.S (0 .. Gen.I - 1);

      end if;

   end Save;

   -----------

   -- Image --

   -----------

   function Image (Of_State : State) return String is

      Result : Image_String;

   begin

      Result := (others => ' ');

      for J in Of_State'Range loop

         Insert_Image (Result, J, Of_State (J));

      end loop;

      return Result;

   end Image;

   function Image (Gen : Generator) return String is

      Result : Image_String;

   begin

      Result := (others => ' ');

      for J in 0 .. N - 1 loop

         Insert_Image (Result, J, Gen.S ((J + Gen.I) mod N));

      end loop;

      return Result;

   end Image;

   -----------

   -- Value --

   -----------

   function Value (Coded_State : String) return State is

      Gen : Generator;

      S   : State;

   begin

      Reset (Gen, Coded_State);

      Save (Gen, S);

      return S;

   end Value;

   ----------

   -- Init --

   ----------

   procedure Init (Gen : Generator; Initiator : Unsigned_32) is

      G : Generator renames Gen.Writable.Self.all;

   begin

      G.S (0) := Initiator;

      for I in 1 .. N - 1 loop

         G.S (I) :=

           (G.S (I - 1) xor Shift_Right (G.S (I - 1), 30)) * Mult0

           + Unsigned_32 (I);

      end loop;

      G.I := 0;

   end Init;

   ------------------

   -- Insert_Image --

   ------------------

   procedure Insert_Image

     (S     : in out Image_String;

      Index : Integer;

      V     : State_Val)

   is

      Value : constant String := State_Val'Image (V);

   begin

      S (Index * 11 + 1 .. Index * 11 + Value'Length) := Value;

   end Insert_Image;

   -------------------

   -- Extract_Value --

   -------------------

   function Extract_Value (S : String; Index : Integer) return State_Val is

      Start : constant Integer := S'First + Index * Image_Numeral_Length;

   begin

      return State_Val'Value (S (Start .. Start + Image_Numeral_Length - 1));

   end Extract_Value;

end System.Random_Numbers;