gcc/ada/s-rannum.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUN-TIME COMPONENTS                         --
   4 --                                                                          --
   5 --                S Y S T E M . R A N D O M _ N U M B E R S                 --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                                                                          --
   9 --          Copyright (C) 2007,2009  Free Software Foundation, Inc.         --
  10 --                                                                          --
  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 3,  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.                                     --
  17 --                                                                          --
  18 -- As a special exception under Section 7 of GPL version 3, you are granted --
  19 -- additional permissions described in the GCC Runtime Library Exception,   --
  20 -- version 3.1, as published by the Free Software Foundation.               --
  21 --                                                                          --
  22 -- You should have received a copy of the GNU General Public License and    --
  23 -- a copy of the GCC Runtime Library Exception along with this program;     --
  24 -- see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see    --
  25 -- <http://www.gnu.org/licenses/>.                                          --
  26 --                                                                          --
  27 -- GNAT was originally developed  by the GNAT team at  New York University. --
  28 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  29 --                                                                          --
  30 ------------------------------------------------------------------------------
  31
  32 ------------------------------------------------------------------------------
  33 --                                                                          --
  34 -- The implementation here is derived from a C-program for MT19937, with    --
  35 -- initialization improved 2002/1/26. As required, the following notice is  --
  36 -- copied from the original program.                                        --
  37 --                                                                          --
  38 -- Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,        --
  39 -- All rights reserved.                                                     --
  40 --                                                                          --
  41 -- Redistribution and use in source and binary forms, with or without       --
  42 -- modification, are permitted provided that the following conditions       --
  43 -- are met:                                                                 --
  44 --                                                                          --
  45 --   1. Redistributions of source code must retain the above copyright      --
  46 --      notice, this list of conditions and the following disclaimer.       --
  47 --                                                                          --
  48 --   2. Redistributions in binary form must reproduce the above copyright   --
  49 --      notice, this list of conditions and the following disclaimer in the --
  50 --      documentation and/or other materials provided with the distribution.--
  51 --                                                                          --
  52 --   3. The names of its contributors may not be used to endorse or promote --
  53 --      products derived from this software without specific prior written  --
  54 --      permission.                                                         --
  55 --                                                                          --
  56 -- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS      --
  57 -- "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT        --
  58 -- LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR    --
  59 -- A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT    --
  60 -- OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,    --
  61 -- SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED --
  62 -- TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR   --
  63 -- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF   --
  64 -- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING     --
  65 -- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS       --
  66 -- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.             --
  67 --                                                                          --
  68 ------------------------------------------------------------------------------
  69
  70 ------------------------------------------------------------------------------
  71 --                                                                          --
  72 -- This is an implementation of the Mersenne Twister, twisted generalized   --
  73 -- feedback shift register of rational normal form, with state-bit          --
  74 -- reflection and tempering. This version generates 32-bit integers with a  --
  75 -- period of 2**19937 - 1 (a Mersenne prime, hence the name). For           --
  76 -- applications requiring more than 32 bits (up to 64), we concatenate two  --
  77 -- 32-bit numbers.                                                          --
  78 --                                                                          --
  79 -- See http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html for         --
  80 -- details.                                                                 --
  81 --                                                                          --
  82 -- In contrast to the original code, we do not generate random numbers in   --
  83 -- batches of N. Measurement seems to show this has very little if any      --
  84 -- effect on performance, and it may be marginally better for real-time     --
  85 -- applications with hard deadlines.                                        --
  86 --                                                                          --
  87 ------------------------------------------------------------------------------
  88
  89 with Ada.Calendar;              use Ada.Calendar;
  90 with Ada.Unchecked_Conversion;
  91 with Interfaces;                use Interfaces;
  92
  93 use Ada;
  94
  95 package body System.Random_Numbers is
  96
  97    -------------------------
  98    -- Implementation Note --
  99    -------------------------
 100
 101    --  The design of this spec is very awkward, as a result of Ada 95 not
 102    --  permitting in-out parameters for function formals (most naturally,
 103    --  Generator values would be passed this way). In pure Ada 95, the only
 104    --  solution is to use the heap and pointers, and, to avoid memory leaks,
 105    --  controlled types.
 106
 107    --  This is awfully heavy, so what we do is to use Unrestricted_Access to
 108    --  get a pointer to the state in the passed Generator. This works because
 109    --  Generator is a limited type and will thus always be passed by reference.
 110
 111    Low31_Mask : constant := 2**31-1;
 112    Bit31_Mask : constant := 2**31;
 113
 114    Matrix_A_X : constant array (State_Val range 0 .. 1) of State_Val :=
 115                   (0, 16#9908b0df#);
 116
 117    Y2K : constant Calendar.Time :=
 118            Calendar.Time_Of
 119              (Year => 2000, Month => 1, Day => 1, Seconds => 0.0);
 120    --  First Year 2000 day
 121
 122    subtype Image_String is String (1 .. Max_Image_Width);
 123
 124    --  Utility functions
 125
 126    procedure Init (Gen : out Generator; Initiator : Unsigned_32);
 127    --  Perform a default initialization of the state of Gen. The resulting
 128    --  state is identical for identical values of Initiator.
 129
 130    procedure Insert_Image
 131      (S     : in out Image_String;
 132       Index : Integer;
 133       V     : State_Val);
 134    --  Insert image of V into S, in the Index'th 11-character substring
 135
 136    function Extract_Value (S : String; Index : Integer) return State_Val;
 137    --  Treat S as a sequence of 11-character decimal numerals and return
 138    --  the result of converting numeral #Index (numbering from 0)
 139
 140    function To_Unsigned is
 141      new Unchecked_Conversion (Integer_32, Unsigned_32);
 142    function To_Unsigned is
 143      new Unchecked_Conversion (Integer_64, Unsigned_64);
 144
 145    ------------
 146    -- Random --
 147    ------------
 148
 149    function Random (Gen : Generator) return Unsigned_32 is
 150       G : Generator renames Gen'Unrestricted_Access.all;
 151       Y : State_Val;
 152       I : Integer;
 153
 154    begin
 155       I := G.I;
 156
 157       if I < N - M then
 158          Y := (G.S (I) and Bit31_Mask) or (G.S (I + 1) and Low31_Mask);
 159          Y := G.S (I + M) xor Shift_Right (Y, 1) xor Matrix_A_X (Y and 1);
 160          I := I + 1;
 161
 162       elsif I < N - 1 then
 163          Y := (G.S (I) and Bit31_Mask) or (G.S (I + 1) and Low31_Mask);
 164          Y := G.S (I + (M - N))
 165                 xor Shift_Right (Y, 1)
 166                 xor Matrix_A_X (Y and 1);
 167          I := I + 1;
 168
 169       elsif I = N - 1 then
 170          Y := (G.S (I) and Bit31_Mask) or (G.S (0) and Low31_Mask);
 171          Y := G.S (M - 1) xor Shift_Right (Y, 1) xor Matrix_A_X (Y and 1);
 172          I := 0;
 173
 174       else
 175          Init (G, 5489);
 176          return Random (Gen);
 177       end if;
 178
 179       G.S (G.I) := Y;
 180       G.I := I;
 181
 182       Y := Y xor Shift_Right (Y, 11);
 183       Y := Y xor (Shift_Left (Y, 7)  and 16#9d2c5680#);
 184       Y := Y xor (Shift_Left (Y, 15) and 16#efc60000#);
 185       Y := Y xor Shift_Right (Y, 18);
 186
 187       return Y;
 188    end Random;
 189
 190    function Random (Gen : Generator) return Float is
 191
 192       --  Note: The application of Float'Machine (...) is necessary to avoid
 193       --  returning extra significand bits. Without it, the function's value
 194       --  will change if it is spilled, for example, causing
 195       --  gratuitous nondeterminism.
 196
 197       Result : constant Float :=
 198                  Float'Machine
 199                    (Float (Unsigned_32'(Random (Gen))) * 2.0 ** (-32));
 200    begin
 201       if Result < 1.0 then
 202          return Result;
 203       else
 204          return Float'Adjacent (1.0, 0.0);
 205       end if;
 206    end Random;
 207
 208    function Random (Gen : Generator) return Long_Float is
 209       Result : constant Long_Float :=
 210                  Long_Float'Machine ((Long_Float (Unsigned_32'(Random (Gen)))
 211                    * 2.0 ** (-32))
 212                    + (Long_Float (Unsigned_32'(Random (Gen))) * 2.0 ** (-64)));
 213    begin
 214       if Result < 1.0 then
 215          return Result;
 216       else
 217          return Long_Float'Adjacent (1.0, 0.0);
 218       end if;
 219    end Random;
 220
 221    function Random (Gen : Generator) return Unsigned_64 is
 222    begin
 223       return Shift_Left (Unsigned_64 (Unsigned_32'(Random (Gen))), 32)
 224         or Unsigned_64 (Unsigned_32'(Random (Gen)));
 225    end Random;
 226
 227    ---------------------
 228    -- Random_Discrete --
 229    ---------------------
 230
 231    function Random_Discrete
 232      (Gen : Generator;
 233       Min : Result_Subtype := Default_Min;
 234       Max : Result_Subtype := Result_Subtype'Last) return Result_Subtype
 235    is
 236    begin
 237       if Max = Min then
 238          return Max;
 239
 240       elsif Max < Min then
 241          raise Constraint_Error;
 242
 243       elsif Result_Subtype'Base'Size > 32 then
 244          declare
 245             --  In the 64-bit case, we have to be careful, since not all 64-bit
 246             --  unsigned values are representable in GNAT's root_integer type.
 247             --  Ignore different-size warnings here; since GNAT's handling
 248             --  is correct.
 249
 250             pragma Warnings ("Z");
 251             function Conv_To_Unsigned is
 252                new Unchecked_Conversion (Result_Subtype'Base, Unsigned_64);
 253             function Conv_To_Result is
 254                new Unchecked_Conversion (Unsigned_64, Result_Subtype'Base);
 255             pragma Warnings ("z");
 256
 257             N : constant Unsigned_64 :=
 258                   Conv_To_Unsigned (Max) - Conv_To_Unsigned (Min) + 1;
 259
 260             X, Slop : Unsigned_64;
 261
 262          begin
 263             if N = 0 then
 264                return Conv_To_Result (Conv_To_Unsigned (Min) + Random (Gen));
 265
 266             else
 267                Slop := Unsigned_64'Last rem N + 1;
 268
 269                loop
 270                   X := Random (Gen);
 271                   exit when Slop = N or else X <= Unsigned_64'Last - Slop;
 272                end loop;
 273
 274                return Conv_To_Result (Conv_To_Unsigned (Min) + X rem N);
 275             end if;
 276          end;
 277
 278       elsif Result_Subtype'Pos (Max) - Result_Subtype'Pos (Min) =
 279                                                          2 ** 32 - 1
 280       then
 281          return Result_Subtype'Val
 282            (Result_Subtype'Pos (Min) + Unsigned_32'Pos (Random (Gen)));
 283       else
 284          declare
 285             N    : constant Unsigned_32 :=
 286                      Unsigned_32 (Result_Subtype'Pos (Max) -
 287                                     Result_Subtype'Pos (Min) + 1);
 288             Slop : constant Unsigned_32 := Unsigned_32'Last rem N + 1;
 289             X    : Unsigned_32;
 290
 291          begin
 292             loop
 293                X := Random (Gen);
 294                exit when Slop = N or else X <= Unsigned_32'Last - Slop;
 295             end loop;
 296
 297             return
 298               Result_Subtype'Val
 299                 (Result_Subtype'Pos (Min) + Unsigned_32'Pos (X rem N));
 300          end;
 301       end if;
 302    end Random_Discrete;
 303
 304    ------------------
 305    -- Random_Float --
 306    ------------------
 307
 308    function Random_Float (Gen : Generator) return Result_Subtype is
 309    begin
 310       if Result_Subtype'Base'Digits > Float'Digits then
 311          return Result_Subtype'Machine (Result_Subtype
 312                                          (Long_Float'(Random (Gen))));
 313       else
 314          return Result_Subtype'Machine (Result_Subtype
 315                                          (Float'(Random (Gen))));
 316       end if;
 317    end Random_Float;
 318
 319    -----------
 320    -- Reset --
 321    -----------
 322
 323    procedure Reset (Gen : out Generator) is
 324       X : constant Unsigned_32 := Unsigned_32 ((Calendar.Clock - Y2K) * 64.0);
 325    begin
 326       Init (Gen, X);
 327    end Reset;
 328
 329    procedure Reset (Gen : out Generator; Initiator : Integer_32) is
 330    begin
 331       Init (Gen, To_Unsigned (Initiator));
 332    end Reset;
 333
 334    procedure Reset (Gen : out Generator; Initiator : Unsigned_32) is
 335    begin
 336       Init (Gen, Initiator);
 337    end Reset;
 338
 339    procedure Reset (Gen : out Generator; Initiator : Integer) is
 340    begin
 341       pragma Warnings ("C");
 342       --  This is probably an unnecessary precaution against future change, but
 343       --  since the test is a static expression, no extra code is involved.
 344
 345       if Integer'Size <= 32 then
 346          Init (Gen, To_Unsigned (Integer_32 (Initiator)));
 347
 348       else
 349          declare
 350             Initiator1 : constant Unsigned_64 :=
 351                            To_Unsigned (Integer_64 (Initiator));
 352             Init0      : constant Unsigned_32 :=
 353                            Unsigned_32 (Initiator1 mod 2 ** 32);
 354             Init1      : constant Unsigned_32 :=
 355                            Unsigned_32 (Shift_Right (Initiator1, 32));
 356          begin
 357             Reset (Gen, Initialization_Vector'(Init0, Init1));
 358          end;
 359       end if;
 360
 361       pragma Warnings ("c");
 362    end Reset;
 363
 364    procedure Reset (Gen : out Generator; Initiator : Initialization_Vector) is
 365       I, J : Integer;
 366
 367    begin
 368       Init (Gen, 19650218);
 369       I := 1;
 370       J := 0;
 371
 372       if Initiator'Length > 0 then
 373          for K in reverse 1 .. Integer'Max (N, Initiator'Length) loop
 374             Gen.S (I) :=
 375               (Gen.S (I)
 376                  xor ((Gen.S (I - 1) xor Shift_Right (Gen.S (I - 1), 30))
 377                                                                  * 1664525))
 378               + Initiator (J + Initiator'First) + Unsigned_32 (J);
 379
 380             I := I + 1;
 381             J := J + 1;
 382
 383             if I >= N then
 384                Gen.S (0) := Gen.S (N - 1);
 385                I := 1;
 386             end if;
 387
 388             if J >= Initiator'Length then
 389                J := 0;
 390             end if;
 391          end loop;
 392       end if;
 393
 394       for K in reverse 1 .. N - 1 loop
 395          Gen.S (I) :=
 396            (Gen.S (I) xor ((Gen.S (I - 1)
 397                             xor Shift_Right (Gen.S (I - 1), 30)) * 1566083941))
 398            - Unsigned_32 (I);
 399          I := I + 1;
 400
 401          if I >= N then
 402             Gen.S (0) := Gen.S (N - 1);
 403             I := 1;
 404          end if;
 405       end loop;
 406
 407       Gen.S (0) := Bit31_Mask;
 408    end Reset;
 409
 410    procedure Reset (Gen : out Generator; From_State : Generator) is
 411    begin
 412       Gen.S := From_State.S;
 413       Gen.I := From_State.I;
 414    end Reset;
 415
 416    procedure Reset (Gen : out Generator; From_State : State) is
 417    begin
 418       Gen.I := 0;
 419       Gen.S := From_State;
 420    end Reset;
 421
 422    procedure Reset (Gen : out Generator; From_Image : String) is
 423    begin
 424       Gen.I := 0;
 425
 426       for J in 0 .. N - 1 loop
 427          Gen.S (J) := Extract_Value (From_Image, J);
 428       end loop;
 429    end Reset;
 430
 431    ----------
 432    -- Save --
 433    ----------
 434
 435    procedure Save (Gen : Generator; To_State : out State) is
 436       Gen2 : Generator;
 437
 438    begin
 439       if Gen.I = N then
 440          Init (Gen2, 5489);
 441          To_State := Gen2.S;
 442
 443       else
 444          To_State (0 .. N - 1 - Gen.I) := Gen.S (Gen.I .. N - 1);
 445          To_State (N - Gen.I .. N - 1) := Gen.S (0 .. Gen.I - 1);
 446       end if;
 447    end Save;
 448
 449    -----------
 450    -- Image --
 451    -----------
 452
 453    function Image (Of_State : State) return String is
 454       Result : Image_String;
 455
 456    begin
 457       Result := (others => ' ');
 458
 459       for J in Of_State'Range loop
 460          Insert_Image (Result, J, Of_State (J));
 461       end loop;
 462
 463       return Result;
 464    end Image;
 465
 466    function Image (Gen : Generator) return String is
 467       Result : Image_String;
 468
 469    begin
 470       Result := (others => ' ');
 471
 472       for J in 0 .. N - 1 loop
 473          Insert_Image (Result, J, Gen.S ((J + Gen.I) mod N));
 474       end loop;
 475
 476       return Result;
 477    end Image;
 478
 479    -----------
 480    -- Value --
 481    -----------
 482
 483    function Value (Coded_State : String) return State is
 484       Gen : Generator;
 485       S   : State;
 486    begin
 487       Reset (Gen, Coded_State);
 488       Save (Gen, S);
 489       return S;
 490    end Value;
 491
 492    ----------
 493    -- Init --
 494    ----------
 495
 496    procedure Init (Gen : out Generator; Initiator : Unsigned_32) is
 497    begin
 498       Gen.S (0) := Initiator;
 499
 500       for I in 1 .. N - 1 loop
 501          Gen.S (I) :=
 502            1812433253
 503              * (Gen.S (I - 1) xor Shift_Right (Gen.S (I - 1), 30))
 504            + Unsigned_32 (I);
 505       end loop;
 506
 507       Gen.I := 0;
 508    end Init;
 509
 510    ------------------
 511    -- Insert_Image --
 512    ------------------
 513
 514    procedure Insert_Image
 515      (S     : in out Image_String;
 516       Index : Integer;
 517       V     : State_Val)
 518    is
 519       Value : constant String := State_Val'Image (V);
 520    begin
 521       S (Index * 11 + 1 .. Index * 11 + Value'Length) := Value;
 522    end Insert_Image;
 523
 524    -------------------
 525    -- Extract_Value --
 526    -------------------
 527
 528    function Extract_Value (S : String; Index : Integer) return State_Val is
 529    begin
 530       return State_Val'Value (S (S'First + Index * 11 ..
 531                                  S'First + Index * 11 + 11));
 532    end Extract_Value;
 533
 534 end System.Random_Numbers;