| ------------------------------------------------------------------------------ |
| -- -- |
| -- 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, 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 2, 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. See the GNU General Public License -- |
| -- for more details. You should have received a copy of the GNU General -- |
| -- Public License distributed with GNAT; see file COPYING. If not, write -- |
| -- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- |
| -- Boston, MA 02110-1301, USA. -- |
| -- -- |
| -- As a special exception, if other files instantiate generics from this -- |
| -- unit, or you link this unit with other files to produce an executable, -- |
| -- this unit does not by itself cause the resulting executable to be -- |
| -- covered by the GNU General Public License. This exception does not -- |
| -- however invalidate any other reasons why the executable file might be -- |
| -- covered by the GNU Public License. -- |
| -- -- |
| -- 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.Calendar; use Ada.Calendar; |
| with Ada.Unchecked_Conversion; |
| with Interfaces; use Interfaces; |
| |
| use Ada; |
| |
| package body System.Random_Numbers is |
| |
| ------------------------- |
| -- Implementation Note -- |
| ------------------------- |
| |
| -- The design of this spec is very awkward, as a result of Ada 95 not |
| -- permitting in-out parameters for function formals (most naturally, |
| -- Generator values would be passed this way). In pure Ada 95, the only |
| -- solution is to use the heap and pointers, and, to avoid memory leaks, |
| -- controlled types. |
| |
| -- This is awfully heavy, so what we do is to use Unrestricted_Access to |
| -- get a pointer to the state in the passed Generator. This works because |
| -- Generator is a limited type and will thus always be passed by reference. |
| |
| Low31_Mask : constant := 2**31-1; |
| Bit31_Mask : constant := 2**31; |
| |
| Matrix_A_X : constant array (State_Val range 0 .. 1) of State_Val := |
| (0, 16#9908b0df#); |
| |
| Y2K : constant Calendar.Time := |
| Calendar.Time_Of |
| (Year => 2000, Month => 1, Day => 1, Seconds => 0.0); |
| -- First Year 2000 day |
| |
| subtype Image_String is String (1 .. Max_Image_Width); |
| |
| -- Utility functions |
| |
| procedure Init (Gen : out 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'Unrestricted_Access.all; |
| Y : State_Val; |
| I : Integer; |
| |
| begin |
| I := G.I; |
| |
| if I < N - M then |
| Y := (G.S (I) and Bit31_Mask) or (G.S (I + 1) and Low31_Mask); |
| Y := G.S (I + M) xor Shift_Right (Y, 1) xor Matrix_A_X (Y and 1); |
| I := I + 1; |
| |
| elsif I < N - 1 then |
| Y := (G.S (I) and Bit31_Mask) or (G.S (I + 1) and Low31_Mask); |
| Y := G.S (I + (M - N)) |
| xor Shift_Right (Y, 1) |
| xor Matrix_A_X (Y and 1); |
| I := I + 1; |
| |
| elsif I = N - 1 then |
| Y := (G.S (I) and Bit31_Mask) or (G.S (0) and Low31_Mask); |
| Y := G.S (M - 1) xor Shift_Right (Y, 1) xor Matrix_A_X (Y and 1); |
| I := 0; |
| |
| else |
| Init (G, 5489); |
| return Random (Gen); |
| end if; |
| |
| G.S (G.I) := Y; |
| G.I := I; |
| |
| Y := Y xor Shift_Right (Y, 11); |
| Y := Y xor (Shift_Left (Y, 7) and 16#9d2c5680#); |
| Y := Y xor (Shift_Left (Y, 15) and 16#efc60000#); |
| Y := Y xor Shift_Right (Y, 18); |
| |
| return Y; |
| end Random; |
| |
| function Random (Gen : Generator) return Float is |
| |
| -- Note: The application of Float'Machine (...) is necessary to avoid |
| -- returning extra significand bits. Without it, the function's value |
| -- will change if it is spilled, for example, causing |
| -- gratuitous nondeterminism. |
| |
| Result : constant Float := |
| Float'Machine |
| (Float (Unsigned_32'(Random (Gen))) * 2.0 ** (-32)); |
| begin |
| if Result < 1.0 then |
| return Result; |
| else |
| return Float'Adjacent (1.0, 0.0); |
| end if; |
| end Random; |
| |
| function Random (Gen : Generator) return Long_Float is |
| Result : constant Long_Float := |
| Long_Float'Machine ((Long_Float (Unsigned_32'(Random (Gen))) |
| * 2.0 ** (-32)) |
| + (Long_Float (Unsigned_32'(Random (Gen))) * 2.0 ** (-64))); |
| begin |
| if Result < 1.0 then |
| return Result; |
| else |
| return Long_Float'Adjacent (1.0, 0.0); |
| end if; |
| 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 : out Generator) is |
| X : constant Unsigned_32 := Unsigned_32 ((Calendar.Clock - Y2K) * 64.0); |
| begin |
| Init (Gen, X); |
| end Reset; |
| |
| procedure Reset (Gen : out Generator; Initiator : Integer_32) is |
| begin |
| Init (Gen, To_Unsigned (Initiator)); |
| end Reset; |
| |
| procedure Reset (Gen : out Generator; Initiator : Unsigned_32) is |
| begin |
| Init (Gen, Initiator); |
| end Reset; |
| |
| procedure Reset (Gen : out Generator; Initiator : Integer) is |
| begin |
| pragma Warnings ("C"); |
| -- 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 ("c"); |
| end Reset; |
| |
| procedure Reset (Gen : out Generator; Initiator : Initialization_Vector) is |
| I, J : Integer; |
| |
| begin |
| Init (Gen, 19650218); |
| I := 1; |
| J := 0; |
| |
| if Initiator'Length > 0 then |
| for K in reverse 1 .. Integer'Max (N, Initiator'Length) loop |
| Gen.S (I) := |
| (Gen.S (I) |
| xor ((Gen.S (I - 1) xor Shift_Right (Gen.S (I - 1), 30)) |
| * 1664525)) |
| + Initiator (J + Initiator'First) + Unsigned_32 (J); |
| |
| I := I + 1; |
| J := J + 1; |
| |
| if I >= N then |
| Gen.S (0) := Gen.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 |
| Gen.S (I) := |
| (Gen.S (I) xor ((Gen.S (I - 1) |
| xor Shift_Right (Gen.S (I - 1), 30)) * 1566083941)) |
| - Unsigned_32 (I); |
| I := I + 1; |
| |
| if I >= N then |
| Gen.S (0) := Gen.S (N - 1); |
| I := 1; |
| end if; |
| end loop; |
| |
| Gen.S (0) := Bit31_Mask; |
| end Reset; |
| |
| procedure Reset (Gen : out Generator; From_State : Generator) is |
| begin |
| Gen.S := From_State.S; |
| Gen.I := From_State.I; |
| end Reset; |
| |
| procedure Reset (Gen : out Generator; From_State : State) is |
| begin |
| Gen.I := 0; |
| Gen.S := From_State; |
| end Reset; |
| |
| procedure Reset (Gen : out Generator; From_Image : String) is |
| begin |
| Gen.I := 0; |
| |
| for J in 0 .. N - 1 loop |
| Gen.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 : out Generator; Initiator : Unsigned_32) is |
| begin |
| Gen.S (0) := Initiator; |
| |
| for I in 1 .. N - 1 loop |
| Gen.S (I) := |
| 1812433253 |
| * (Gen.S (I - 1) xor Shift_Right (Gen.S (I - 1), 30)) |
| + Unsigned_32 (I); |
| end loop; |
| |
| Gen.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 |
| begin |
| return State_Val'Value (S (S'First + Index * 11 .. |
| S'First + Index * 11 + 11)); |
| end Extract_Value; |
| |
| end System.Random_Numbers; |