------------------------------------------------------------------------------ -- -- -- 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-2018, 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 -- -- . -- -- -- -- 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 with SPARK_Mode => Off 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 nonbinary radix): computing Real() / (+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"); 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 -- 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; 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;