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, 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 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.                                              --
  21 --                                                                          --
  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.                                      --
  28 --                                                                          --
  29 -- GNAT was originally developed  by the GNAT team at  New York University. --
  30 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  31 --                                                                          --
  32 ------------------------------------------------------------------------------
  33
  34 ------------------------------------------------------------------------------
  35 --                                                                          --
  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.                                        --
  39 --                                                                          --
  40 -- Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,        --
  41 -- All rights reserved.                                                     --
  42 --                                                                          --
  43 -- Redistribution and use in source and binary forms, with or without       --
  44 -- modification, are permitted provided that the following conditions       --
  45 -- are met:                                                                 --
  46 --                                                                          --
  47 --   1. Redistributions of source code must retain the above copyright      --
  48 --      notice, this list of conditions and the following disclaimer.       --
  49 --                                                                          --
  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.--
  53 --                                                                          --
  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  --
  56 --      permission.                                                         --
  57 --                                                                          --
  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.             --
  69 --                                                                          --
  70 ------------------------------------------------------------------------------
  71
  72 ------------------------------------------------------------------------------
  73 --                                                                          --
  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  --
  79 -- 32-bit numbers.                                                          --
  80 --                                                                          --
  81 -- See http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html for         --
  82 -- details.                                                                 --
  83 --                                                                          --
  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.                                        --
  88 --                                                                          --
  89 ------------------------------------------------------------------------------
  90
  91 with Ada.Calendar;              use Ada.Calendar;
  92 with Ada.Unchecked_Conversion;
  93 with Interfaces;                use Interfaces;
  94
  95 use Ada;
  96
  97 package body System.Random_Numbers is
  98
  99    -------------------------
 100    -- Implementation Note --
 101    -------------------------
 102
 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,
 107    --  controlled types.
 108
 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.
 112
 113    Low31_Mask : constant := 2**31-1;
 114    Bit31_Mask : constant := 2**31;
 115
 116    Matrix_A_X : constant array (State_Val range 0 .. 1) of State_Val :=
 117                   (0, 16#9908b0df#);
 118
 119    Y2K : constant Calendar.Time :=
 120            Calendar.Time_Of
 121              (Year => 2000, Month => 1, Day => 1, Seconds => 0.0);
 122    --  First Year 2000 day
 123
 124    subtype Image_String is String (1 .. Max_Image_Width);
 125
 126    --  Utility functions
 127
 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.
 131
 132    procedure Insert_Image
 133      (S     : in out Image_String;
 134       Index : Integer;
 135       V     : State_Val);
 136    --  Insert image of V into S, in the Index'th 11-character substring
 137
 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)
 141
 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);
 146
 147    ------------
 148    -- Random --
 149    ------------
 150
 151    function Random (Gen : Generator) return Unsigned_32 is
 152       G : Generator renames Gen'Unrestricted_Access.all;
 153       Y : State_Val;
 154       I : Integer;
 155
 156    begin
 157       I := G.I;
 158
 159       if I < N - M then
 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);
 162          I := I + 1;
 163
 164       elsif I < N - 1 then
 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);
 169          I := I + 1;
 170
 171       elsif I = N - 1 then
 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);
 174          I := 0;
 175
 176       else
 177          Init (G, 5489);
 178          return Random (Gen);
 179       end if;
 180
 181       G.S (G.I) := Y;
 182       G.I := I;
 183
 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);
 188
 189       return Y;
 190    end Random;
 191
 192    function Random (Gen : Generator) return Float is
 193
 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.
 198
 199       Result : constant Float :=
 200                  Float'Machine
 201                    (Float (Unsigned_32'(Random (Gen))) * 2.0 ** (-32));
 202    begin
 203       if Result < 1.0 then
 204          return Result;
 205       else
 206          return Float'Adjacent (1.0, 0.0);
 207       end if;
 208    end Random;
 209
 210    function Random (Gen : Generator) return Long_Float is
 211       Result : constant Long_Float :=
 212                  Long_Float'Machine ((Long_Float (Unsigned_32'(Random (Gen)))
 213                    * 2.0 ** (-32))
 214                    + (Long_Float (Unsigned_32'(Random (Gen))) * 2.0 ** (-64)));
 215    begin
 216       if Result < 1.0 then
 217          return Result;
 218       else
 219          return Long_Float'Adjacent (1.0, 0.0);
 220       end if;
 221    end Random;
 222
 223    function Random (Gen : Generator) return Unsigned_64 is
 224    begin
 225       return Shift_Left (Unsigned_64 (Unsigned_32'(Random (Gen))), 32)
 226         or Unsigned_64 (Unsigned_32'(Random (Gen)));
 227    end Random;
 228
 229    ---------------------
 230    -- Random_Discrete --
 231    ---------------------
 232
 233    function Random_Discrete
 234      (Gen : Generator;
 235       Min : Result_Subtype := Default_Min;
 236       Max : Result_Subtype := Result_Subtype'Last) return Result_Subtype
 237    is
 238    begin
 239       if Max = Min then
 240          return Max;
 241
 242       elsif Max < Min then
 243          raise Constraint_Error;
 244
 245       elsif Result_Subtype'Base'Size > 32 then
 246          declare
 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
 250             --  is correct.
 251
 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");
 258
 259             N : constant Unsigned_64 :=
 260                   Conv_To_Unsigned (Max) - Conv_To_Unsigned (Min) + 1;
 261
 262             X, Slop : Unsigned_64;
 263
 264          begin
 265             if N = 0 then
 266                return Conv_To_Result (Conv_To_Unsigned (Min) + Random (Gen));
 267
 268             else
 269                Slop := Unsigned_64'Last rem N + 1;
 270
 271                loop
 272                   X := Random (Gen);
 273                   exit when Slop = N or else X <= Unsigned_64'Last - Slop;
 274                end loop;
 275
 276                return Conv_To_Result (Conv_To_Unsigned (Min) + X rem N);
 277             end if;
 278          end;
 279
 280       elsif Result_Subtype'Pos (Max) - Result_Subtype'Pos (Min) =
 281                                                          2 ** 32 - 1
 282       then
 283          return Result_Subtype'Val
 284            (Result_Subtype'Pos (Min) + Unsigned_32'Pos (Random (Gen)));
 285       else
 286          declare
 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;
 291             X    : Unsigned_32;
 292
 293          begin
 294             loop
 295                X := Random (Gen);
 296                exit when Slop = N or else X <= Unsigned_32'Last - Slop;
 297             end loop;
 298
 299             return
 300               Result_Subtype'Val
 301                 (Result_Subtype'Pos (Min) + Unsigned_32'Pos (X rem N));
 302          end;
 303       end if;
 304    end Random_Discrete;
 305
 306    ------------------
 307    -- Random_Float --
 308    ------------------
 309
 310    function Random_Float (Gen : Generator) return Result_Subtype is
 311    begin
 312       if Result_Subtype'Base'Digits > Float'Digits then
 313          return Result_Subtype'Machine (Result_Subtype
 314                                          (Long_Float'(Random (Gen))));
 315       else
 316          return Result_Subtype'Machine (Result_Subtype
 317                                          (Float'(Random (Gen))));
 318       end if;
 319    end Random_Float;
 320
 321    -----------
 322    -- Reset --
 323    -----------
 324
 325    procedure Reset (Gen : out Generator) is
 326       X : constant Unsigned_32 := Unsigned_32 ((Calendar.Clock - Y2K) * 64.0);
 327    begin
 328       Init (Gen, X);
 329    end Reset;
 330
 331    procedure Reset (Gen : out Generator; Initiator : Integer_32) is
 332    begin
 333       Init (Gen, To_Unsigned (Initiator));
 334    end Reset;
 335
 336    procedure Reset (Gen : out Generator; Initiator : Unsigned_32) is
 337    begin
 338       Init (Gen, Initiator);
 339    end Reset;
 340
 341    procedure Reset (Gen : out Generator; Initiator : Integer) is
 342    begin
 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.
 346
 347       if Integer'Size <= 32 then
 348          Init (Gen, To_Unsigned (Integer_32 (Initiator)));
 349
 350       else
 351          declare
 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));
 358          begin
 359             Reset (Gen, Initialization_Vector'(Init0, Init1));
 360          end;
 361       end if;
 362
 363       pragma Warnings ("c");
 364    end Reset;
 365
 366    procedure Reset (Gen : out Generator; Initiator : Initialization_Vector) is
 367       I, J : Integer;
 368
 369    begin
 370       Init (Gen, 19650218);
 371       I := 1;
 372       J := 0;
 373
 374       if Initiator'Length > 0 then
 375          for K in reverse 1 .. Integer'Max (N, Initiator'Length) loop
 376             Gen.S (I) :=
 377               (Gen.S (I)
 378                  xor ((Gen.S (I - 1) xor Shift_Right (Gen.S (I - 1), 30))
 379                                                                  * 1664525))
 380               + Initiator (J + Initiator'First) + Unsigned_32 (J);
 381
 382             I := I + 1;
 383             J := J + 1;
 384
 385             if I >= N then
 386                Gen.S (0) := Gen.S (N - 1);
 387                I := 1;
 388             end if;
 389
 390             if J >= Initiator'Length then
 391                J := 0;
 392             end if;
 393          end loop;
 394       end if;
 395
 396       for K in reverse 1 .. N - 1 loop
 397          Gen.S (I) :=
 398            (Gen.S (I) xor ((Gen.S (I - 1)
 399                             xor Shift_Right (Gen.S (I - 1), 30)) * 1566083941))
 400            - Unsigned_32 (I);
 401          I := I + 1;
 402
 403          if I >= N then
 404             Gen.S (0) := Gen.S (N - 1);
 405             I := 1;
 406          end if;
 407       end loop;
 408
 409       Gen.S (0) := Bit31_Mask;
 410    end Reset;
 411
 412    procedure Reset (Gen : out Generator; From_State : Generator) is
 413    begin
 414       Gen.S := From_State.S;
 415       Gen.I := From_State.I;
 416    end Reset;
 417
 418    procedure Reset (Gen : out Generator; From_State : State) is
 419    begin
 420       Gen.I := 0;
 421       Gen.S := From_State;
 422    end Reset;
 423
 424    procedure Reset (Gen : out Generator; From_Image : String) is
 425    begin
 426       Gen.I := 0;
 427
 428       for J in 0 .. N - 1 loop
 429          Gen.S (J) := Extract_Value (From_Image, J);
 430       end loop;
 431    end Reset;
 432
 433    ----------
 434    -- Save --
 435    ----------
 436
 437    procedure Save (Gen : Generator; To_State : out State) is
 438       Gen2 : Generator;
 439
 440    begin
 441       if Gen.I = N then
 442          Init (Gen2, 5489);
 443          To_State := Gen2.S;
 444
 445       else
 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);
 448       end if;
 449    end Save;
 450
 451    -----------
 452    -- Image --
 453    -----------
 454
 455    function Image (Of_State : State) return String is
 456       Result : Image_String;
 457
 458    begin
 459       Result := (others => ' ');
 460
 461       for J in Of_State'Range loop
 462          Insert_Image (Result, J, Of_State (J));
 463       end loop;
 464
 465       return Result;
 466    end Image;
 467
 468    function Image (Gen : Generator) return String is
 469       Result : Image_String;
 470
 471    begin
 472       Result := (others => ' ');
 473
 474       for J in 0 .. N - 1 loop
 475          Insert_Image (Result, J, Gen.S ((J + Gen.I) mod N));
 476       end loop;
 477
 478       return Result;
 479    end Image;
 480
 481    -----------
 482    -- Value --
 483    -----------
 484
 485    function Value (Coded_State : String) return State is
 486       Gen : Generator;
 487       S   : State;
 488    begin
 489       Reset (Gen, Coded_State);
 490       Save (Gen, S);
 491       return S;
 492    end Value;
 493
 494    ----------
 495    -- Init --
 496    ----------
 497
 498    procedure Init (Gen : out Generator; Initiator : Unsigned_32) is
 499    begin
 500       Gen.S (0) := Initiator;
 501
 502       for I in 1 .. N - 1 loop
 503          Gen.S (I) :=
 504            1812433253
 505              * (Gen.S (I - 1) xor Shift_Right (Gen.S (I - 1), 30))
 506            + Unsigned_32 (I);
 507       end loop;
 508
 509       Gen.I := 0;
 510    end Init;
 511
 512    ------------------
 513    -- Insert_Image --
 514    ------------------
 515
 516    procedure Insert_Image
 517      (S     : in out Image_String;
 518       Index : Integer;
 519       V     : State_Val)
 520    is
 521       Value : constant String := State_Val'Image (V);
 522    begin
 523       S (Index * 11 + 1 .. Index * 11 + Value'Length) := Value;
 524    end Insert_Image;
 525
 526    -------------------
 527    -- Extract_Value --
 528    -------------------
 529
 530    function Extract_Value (S : String; Index : Integer) return State_Val is
 531    begin
 532       return State_Val'Value (S (S'First + Index * 11 ..
 533                                  S'First + Index * 11 + 11));
 534    end Extract_Value;
 535
 536 end System.Random_Numbers;