1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
5 -- S Y S T E M . R A N D O M _ N U M B E R S --
9 -- Copyright (C) 2007, Free Software Foundation, Inc. --
11 -- GNAT is free software; you can redistribute it and/or modify it under --
12 -- terms of the GNU General Public License as published by the Free Soft- --
13 -- ware Foundation; either version 2, or (at your option) any later ver- --
14 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
15 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
16 -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License --
17 -- for more details. You should have received a copy of the GNU General --
18 -- Public License distributed with GNAT; see file COPYING. If not, write --
19 -- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, --
20 -- Boston, MA 02110-1301, USA. --
22 -- As a special exception, if other files instantiate generics from this --
23 -- unit, or you link this unit with other files to produce an executable, --
24 -- this unit does not by itself cause the resulting executable to be --
25 -- covered by the GNU General Public License. This exception does not --
26 -- however invalidate any other reasons why the executable file might be --
27 -- covered by the GNU Public License. --
29 -- GNAT was originally developed by the GNAT team at New York University. --
30 -- Extensive contributions were provided by Ada Core Technologies Inc. --
32 ------------------------------------------------------------------------------
34 ------------------------------------------------------------------------------
36 -- The implementation here is derived from a C-program for MT19937, with --
37 -- initialization improved 2002/1/26. As required, the following notice is --
38 -- copied from the original program. --
40 -- Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, --
41 -- All rights reserved. --
43 -- Redistribution and use in source and binary forms, with or without --
44 -- modification, are permitted provided that the following conditions --
47 -- 1. Redistributions of source code must retain the above copyright --
48 -- notice, this list of conditions and the following disclaimer. --
50 -- 2. Redistributions in binary form must reproduce the above copyright --
51 -- notice, this list of conditions and the following disclaimer in the --
52 -- documentation and/or other materials provided with the distribution.--
54 -- 3. The names of its contributors may not be used to endorse or promote --
55 -- products derived from this software without specific prior written --
58 -- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS --
59 -- "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT --
60 -- LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR --
61 -- A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT --
62 -- OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, --
63 -- SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED --
64 -- TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR --
65 -- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF --
66 -- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING --
67 -- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS --
68 -- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. --
70 ------------------------------------------------------------------------------
72 ------------------------------------------------------------------------------
74 -- This is an implementation of the Mersenne Twister, twisted generalized --
75 -- feedback shift register of rational normal form, with state-bit --
76 -- reflection and tempering. This version generates 32-bit integers with a --
77 -- period of 2**19937 - 1 (a Mersenne prime, hence the name). For --
78 -- applications requiring more than 32 bits (up to 64), we concatenate two --
81 -- See http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html for --
84 -- In contrast to the original code, we do not generate random numbers in --
85 -- batches of N. Measurement seems to show this has very little if any --
86 -- effect on performance, and it may be marginally better for real-time --
87 -- applications with hard deadlines. --
89 ------------------------------------------------------------------------------
91 with Ada.Calendar; use Ada.Calendar;
92 with Ada.Unchecked_Conversion;
93 with Interfaces; use Interfaces;
97 package body System.Random_Numbers is
99 -------------------------
100 -- Implementation Note --
101 -------------------------
103 -- The design of this spec is very awkward, as a result of Ada 95 not
104 -- permitting in-out parameters for function formals (most naturally,
105 -- Generator values would be passed this way). In pure Ada 95, the only
106 -- solution is to use the heap and pointers, and, to avoid memory leaks,
109 -- This is awfully heavy, so what we do is to use Unrestricted_Access to
110 -- get a pointer to the state in the passed Generator. This works because
111 -- Generator is a limited type and will thus always be passed by reference.
113 Low31_Mask : constant := 2**31-1;
114 Bit31_Mask : constant := 2**31;
116 Matrix_A_X : constant array (State_Val range 0 .. 1) of State_Val :=
119 Y2K : constant Calendar.Time :=
121 (Year => 2000, Month => 1, Day => 1, Seconds => 0.0);
122 -- First Year 2000 day
124 subtype Image_String is String (1 .. Max_Image_Width);
128 procedure Init (Gen : out Generator; Initiator : Unsigned_32);
129 -- Perform a default initialization of the state of Gen. The resulting
130 -- state is identical for identical values of Initiator.
132 procedure Insert_Image
133 (S : in out Image_String;
136 -- Insert image of V into S, in the Index'th 11-character substring
138 function Extract_Value (S : String; Index : Integer) return State_Val;
139 -- Treat S as a sequence of 11-character decimal numerals and return
140 -- the result of converting numeral #Index (numbering from 0)
142 function To_Unsigned is
143 new Unchecked_Conversion (Integer_32, Unsigned_32);
144 function To_Unsigned is
145 new Unchecked_Conversion (Integer_64, Unsigned_64);
151 function Random (Gen : Generator) return Unsigned_32 is
152 G : Generator renames Gen'Unrestricted_Access.all;
160 Y := (G.S (I) and Bit31_Mask) or (G.S (I + 1) and Low31_Mask);
161 Y := G.S (I + M) xor Shift_Right (Y, 1) xor Matrix_A_X (Y and 1);
165 Y := (G.S (I) and Bit31_Mask) or (G.S (I + 1) and Low31_Mask);
166 Y := G.S (I + (M - N))
167 xor Shift_Right (Y, 1)
168 xor Matrix_A_X (Y and 1);
172 Y := (G.S (I) and Bit31_Mask) or (G.S (0) and Low31_Mask);
173 Y := G.S (M - 1) xor Shift_Right (Y, 1) xor Matrix_A_X (Y and 1);
184 Y := Y xor Shift_Right (Y, 11);
185 Y := Y xor (Shift_Left (Y, 7) and 16#9d2c5680#);
186 Y := Y xor (Shift_Left (Y, 15) and 16#efc60000#);
187 Y := Y xor Shift_Right (Y, 18);
192 function Random (Gen : Generator) return Float is
194 -- Note: The application of Float'Machine (...) is necessary to avoid
195 -- returning extra significand bits. Without it, the function's value
196 -- will change if it is spilled, for example, causing
197 -- gratuitous nondeterminism.
199 Result : constant Float :=
201 (Float (Unsigned_32'(Random (Gen))) * 2.0 ** (-32));
206 return Float'Adjacent (1.0, 0.0);
210 function Random (Gen : Generator) return Long_Float is
211 Result : constant Long_Float :=
212 Long_Float'Machine ((Long_Float (Unsigned_32'(Random (Gen)))
214 + (Long_Float (Unsigned_32'(Random (Gen))) * 2.0 ** (-64)));
219 return Long_Float'Adjacent (1.0, 0.0);
223 function Random (Gen : Generator) return Unsigned_64 is
225 return Shift_Left (Unsigned_64 (Unsigned_32'(Random (Gen))), 32)
226 or Unsigned_64 (Unsigned_32'(Random (Gen)));
229 ---------------------
230 -- Random_Discrete --
231 ---------------------
233 function Random_Discrete
235 Min : Result_Subtype := Default_Min;
236 Max : Result_Subtype := Result_Subtype'Last) return Result_Subtype
243 raise Constraint_Error;
245 elsif Result_Subtype'Base'Size > 32 then
247 -- In the 64-bit case, we have to be careful, since not all 64-bit
248 -- unsigned values are representable in GNAT's root_integer type.
249 -- Ignore different-size warnings here; since GNAT's handling
252 pragma Warnings ("Z");
253 function Conv_To_Unsigned is
254 new Unchecked_Conversion (Result_Subtype'Base, Unsigned_64);
255 function Conv_To_Result is
256 new Unchecked_Conversion (Unsigned_64, Result_Subtype'Base);
257 pragma Warnings ("z");
259 N : constant Unsigned_64 :=
260 Conv_To_Unsigned (Max) - Conv_To_Unsigned (Min) + 1;
262 X, Slop : Unsigned_64;
266 return Conv_To_Result (Conv_To_Unsigned (Min) + Random (Gen));
269 Slop := Unsigned_64'Last rem N + 1;
273 exit when Slop = N or else X <= Unsigned_64'Last - Slop;
276 return Conv_To_Result (Conv_To_Unsigned (Min) + X rem N);
280 elsif Result_Subtype'Pos (Max) - Result_Subtype'Pos (Min) =
283 return Result_Subtype'Val
284 (Result_Subtype'Pos (Min) + Unsigned_32'Pos (Random (Gen)));
287 N : constant Unsigned_32 :=
288 Unsigned_32 (Result_Subtype'Pos (Max) -
289 Result_Subtype'Pos (Min) + 1);
290 Slop : constant Unsigned_32 := Unsigned_32'Last rem N + 1;
296 exit when Slop = N or else X <= Unsigned_32'Last - Slop;
301 (Result_Subtype'Pos (Min) + Unsigned_32'Pos (X rem N));
310 function Random_Float (Gen : Generator) return Result_Subtype is
312 if Result_Subtype'Base'Digits > Float'Digits then
313 return Result_Subtype'Machine (Result_Subtype
314 (Long_Float'(Random (Gen))));
316 return Result_Subtype'Machine (Result_Subtype
317 (Float'(Random (Gen))));
325 procedure Reset (Gen : out Generator) is
326 X : constant Unsigned_32 := Unsigned_32 ((Calendar.Clock - Y2K) * 64.0);
331 procedure Reset (Gen : out Generator; Initiator : Integer_32) is
333 Init (Gen, To_Unsigned (Initiator));
336 procedure Reset (Gen : out Generator; Initiator : Unsigned_32) is
338 Init (Gen, Initiator);
341 procedure Reset (Gen : out Generator; Initiator : Integer) is
343 pragma Warnings ("C");
344 -- This is probably an unnecessary precaution against future change, but
345 -- since the test is a static expression, no extra code is involved.
347 if Integer'Size <= 32 then
348 Init (Gen, To_Unsigned (Integer_32 (Initiator)));
352 Initiator1 : constant Unsigned_64 :=
353 To_Unsigned (Integer_64 (Initiator));
354 Init0 : constant Unsigned_32 :=
355 Unsigned_32 (Initiator1 mod 2 ** 32);
356 Init1 : constant Unsigned_32 :=
357 Unsigned_32 (Shift_Right (Initiator1, 32));
359 Reset (Gen, Initialization_Vector'(Init0, Init1));
363 pragma Warnings ("c");
366 procedure Reset (Gen : out Generator; Initiator : Initialization_Vector) is
370 Init (Gen, 19650218);
374 if Initiator'Length > 0 then
375 for K in reverse 1 .. Integer'Max (N, Initiator'Length) loop
378 xor ((Gen.S (I - 1) xor Shift_Right (Gen.S (I - 1), 30))
380 + Initiator (J + Initiator'First) + Unsigned_32 (J);
386 Gen.S (0) := Gen.S (N - 1);
390 if J >= Initiator'Length then
396 for K in reverse 1 .. N - 1 loop
398 (Gen.S (I) xor ((Gen.S (I - 1)
399 xor Shift_Right (Gen.S (I - 1), 30)) * 1566083941))
404 Gen.S (0) := Gen.S (N - 1);
409 Gen.S (0) := Bit31_Mask;
412 procedure Reset (Gen : out Generator; From_State : Generator) is
414 Gen.S := From_State.S;
415 Gen.I := From_State.I;
418 procedure Reset (Gen : out Generator; From_State : State) is
424 procedure Reset (Gen : out Generator; From_Image : String) is
428 for J in 0 .. N - 1 loop
429 Gen.S (J) := Extract_Value (From_Image, J);
437 procedure Save (Gen : Generator; To_State : out State) is
446 To_State (0 .. N - 1 - Gen.I) := Gen.S (Gen.I .. N - 1);
447 To_State (N - Gen.I .. N - 1) := Gen.S (0 .. Gen.I - 1);
455 function Image (Of_State : State) return String is
456 Result : Image_String;
459 Result := (others => ' ');
461 for J in Of_State'Range loop
462 Insert_Image (Result, J, Of_State (J));
468 function Image (Gen : Generator) return String is
469 Result : Image_String;
472 Result := (others => ' ');
474 for J in 0 .. N - 1 loop
475 Insert_Image (Result, J, Gen.S ((J + Gen.I) mod N));
485 function Value (Coded_State : String) return State is
489 Reset (Gen, Coded_State);
498 procedure Init (Gen : out Generator; Initiator : Unsigned_32) is
500 Gen.S (0) := Initiator;
502 for I in 1 .. N - 1 loop
505 * (Gen.S (I - 1) xor Shift_Right (Gen.S (I - 1), 30))
516 procedure Insert_Image
517 (S : in out Image_String;
521 Value : constant String := State_Val'Image (V);
523 S (Index * 11 + 1 .. Index * 11 + Value'Length) := Value;
530 function Extract_Value (S : String; Index : Integer) return State_Val is
532 return State_Val'Value (S (S'First + Index * 11 ..
533 S'First + Index * 11 + 11));
536 end System.Random_Numbers;