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-2015, 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.Unchecked_Conversion;
  90
  91 with System.Random_Seed;
  92
  93 with Interfaces; use Interfaces;
  94
  95 use Ada;
  96
  97 package body System.Random_Numbers with
  98   SPARK_Mode => Off
  99 is
 100    Image_Numeral_Length : constant := Max_Image_Width / N;
 101
 102    subtype Image_String is String (1 .. Max_Image_Width);
 103
 104    ----------------------------
 105    -- Algorithmic Parameters --
 106    ----------------------------
 107
 108    Lower_Mask : constant := 2**31 - 1;
 109    Upper_Mask : constant := 2**31;
 110
 111    Matrix_A   : constant array (State_Val range 0 .. 1) of State_Val
 112      := (0, 16#9908b0df#);
 113    --  The twist transformation is represented by a matrix of the form
 114    --
 115    --               [  0    I(31) ]
 116    --               [    _a       ]
 117    --
 118    --  where 0 is a 31x31 block of 0s, I(31) is the 31x31 identity matrix and
 119    --  _a is a particular bit row-vector, represented here by a 32-bit integer.
 120    --  If integer x represents a row vector of bits (with x(0), the units bit,
 121    --  last), then
 122    --           x * A = [0 x(31..1)] xor Matrix_A(x(0)).
 123
 124    U      : constant := 11;
 125    S      : constant := 7;
 126    B_Mask : constant := 16#9d2c5680#;
 127    T      : constant := 15;
 128    C_Mask : constant := 16#efc60000#;
 129    L      : constant := 18;
 130    --  The tempering shifts and bit masks, in the order applied
 131
 132    Seed0 : constant := 5489;
 133    --  Default seed, used to initialize the state vector when Reset not called
 134
 135    Seed1 : constant := 19650218;
 136    --  Seed used to initialize the state vector when calling Reset with an
 137    --  initialization vector.
 138
 139    Mult0 : constant := 1812433253;
 140    --  Multiplier for a modified linear congruential generator used to
 141    --  initialize the state vector when calling Reset with a single integer
 142    --  seed.
 143
 144    Mult1 : constant := 1664525;
 145    Mult2 : constant := 1566083941;
 146    --  Multipliers for two modified linear congruential generators used to
 147    --  initialize the state vector when calling Reset with an initialization
 148    --  vector.
 149
 150    -----------------------
 151    -- Local Subprograms --
 152    -----------------------
 153
 154    procedure Init (Gen : Generator; Initiator : Unsigned_32);
 155    --  Perform a default initialization of the state of Gen. The resulting
 156    --  state is identical for identical values of Initiator.
 157
 158    procedure Insert_Image
 159      (S     : in out Image_String;
 160       Index : Integer;
 161       V     : State_Val);
 162    --  Insert image of V into S, in the Index'th 11-character substring
 163
 164    function Extract_Value (S : String; Index : Integer) return State_Val;
 165    --  Treat S as a sequence of 11-character decimal numerals and return
 166    --  the result of converting numeral #Index (numbering from 0)
 167
 168    function To_Unsigned is
 169      new Unchecked_Conversion (Integer_32, Unsigned_32);
 170    function To_Unsigned is
 171      new Unchecked_Conversion (Integer_64, Unsigned_64);
 172
 173    ------------
 174    -- Random --
 175    ------------
 176
 177    function Random (Gen : Generator) return Unsigned_32 is
 178       G : Generator renames Gen.Writable.Self.all;
 179       Y : State_Val;
 180       I : Integer;      --  should avoid use of identifier I ???
 181
 182    begin
 183       I := G.I;
 184
 185       if I < N - M then
 186          Y := (G.S (I) and Upper_Mask) or (G.S (I + 1) and Lower_Mask);
 187          Y := G.S (I + M) xor Shift_Right (Y, 1) xor Matrix_A (Y and 1);
 188          I := I + 1;
 189
 190       elsif I < N - 1 then
 191          Y := (G.S (I) and Upper_Mask) or (G.S (I + 1) and Lower_Mask);
 192          Y := G.S (I + (M - N))
 193                 xor Shift_Right (Y, 1)
 194                 xor Matrix_A (Y and 1);
 195          I := I + 1;
 196
 197       elsif I = N - 1 then
 198          Y := (G.S (I) and Upper_Mask) or (G.S (0) and Lower_Mask);
 199          Y := G.S (M - 1) xor Shift_Right (Y, 1) xor Matrix_A (Y and 1);
 200          I := 0;
 201
 202       else
 203          Init (G, Seed0);
 204          return Random (Gen);
 205       end if;
 206
 207       G.S (G.I) := Y;
 208       G.I := I;
 209
 210       Y := Y xor Shift_Right (Y, U);
 211       Y := Y xor (Shift_Left (Y, S)  and B_Mask);
 212       Y := Y xor (Shift_Left (Y, T) and C_Mask);
 213       Y := Y xor Shift_Right (Y, L);
 214
 215       return Y;
 216    end Random;
 217
 218    generic
 219       type Unsigned is mod <>;
 220       type Real is digits <>;
 221       with function Random (G : Generator) return Unsigned is <>;
 222    function Random_Float_Template (Gen : Generator) return Real;
 223    pragma Inline (Random_Float_Template);
 224    --  Template for a random-number generator implementation that delivers
 225    --  values of type Real in the range [0 .. 1], using values from Gen,
 226    --  assuming that Unsigned is large enough to hold the bits of a mantissa
 227    --  for type Real.
 228
 229    ---------------------------
 230    -- Random_Float_Template --
 231    ---------------------------
 232
 233    function Random_Float_Template (Gen : Generator) return Real is
 234
 235       pragma Compile_Time_Error
 236         (Unsigned'Last <= 2**(Real'Machine_Mantissa - 1),
 237          "insufficiently large modular type used to hold mantissa");
 238
 239    begin
 240       --  This code generates random floating-point numbers from unsigned
 241       --  integers. Assuming that Real'Machine_Radix = 2, it can deliver all
 242       --  machine values of type Real (as implied by Real'Machine_Mantissa and
 243       --  Real'Machine_Emin), which is not true of the standard method (to
 244       --  which we fall back for non-binary radix): computing Real(<random
 245       --  integer>) / (<max random integer>+1). To do so, we first extract an
 246       --  (M-1)-bit significand (where M is Real'Machine_Mantissa), and then
 247       --  decide on a normalized exponent by repeated coin flips, decrementing
 248       --  from 0 as long as we flip heads (1 bits). This process yields the
 249       --  proper geometric distribution for the exponent: in a uniformly
 250       --  distributed set of floating-point numbers, 1/2 of them will be in
 251       --  (0.5, 1], 1/4 will be in (0.25, 0.5], and so forth. It makes a
 252       --  further adjustment at binade boundaries (see comments below) to give
 253       --  the effect of selecting a uniformly distributed real deviate in
 254       --  [0..1] and then rounding to the nearest representable floating-point
 255       --  number.  The algorithm attempts to be stingy with random integers. In
 256       --  the worst case, it can consume roughly -Real'Machine_Emin/32 32-bit
 257       --  integers, but this case occurs with probability around
 258       --  2**Machine_Emin, and the expected number of calls to integer-valued
 259       --  Random is 1.  For another discussion of the issues addressed by this
 260       --  process, see Allen Downey's unpublished paper at
 261       --  http://allendowney.com/research/rand/downey07randfloat.pdf.
 262
 263       if Real'Machine_Radix /= 2 then
 264          return Real'Machine
 265            (Real (Unsigned'(Random (Gen))) * 2.0**(-Unsigned'Size));
 266
 267       else
 268          declare
 269             type Bit_Count is range 0 .. 4;
 270
 271             subtype T is Real'Base;
 272
 273             Trailing_Ones : constant array (Unsigned_32 range 0 .. 15)
 274               of Bit_Count :=
 275                   (2#00000# => 0, 2#00001# => 1, 2#00010# => 0, 2#00011# => 2,
 276                    2#00100# => 0, 2#00101# => 1, 2#00110# => 0, 2#00111# => 3,
 277                    2#01000# => 0, 2#01001# => 1, 2#01010# => 0, 2#01011# => 2,
 278                    2#01100# => 0, 2#01101# => 1, 2#01110# => 0, 2#01111# => 4);
 279
 280             Pow_Tab : constant array (Bit_Count range 0 .. 3) of Real
 281               := (0 => 2.0**(0 - T'Machine_Mantissa),
 282                   1 => 2.0**(-1 - T'Machine_Mantissa),
 283                   2 => 2.0**(-2 - T'Machine_Mantissa),
 284                   3 => 2.0**(-3 - T'Machine_Mantissa));
 285
 286             Extra_Bits : constant Natural :=
 287                          (Unsigned'Size - T'Machine_Mantissa + 1);
 288             --  Random bits left over after selecting mantissa
 289
 290             Mantissa : Unsigned;
 291
 292             X      : Real;            --  Scaled mantissa
 293             R      : Unsigned_32;     --  Supply of random bits
 294             R_Bits : Natural;         --  Number of bits left in R
 295             K      : Bit_Count;       --  Next decrement to exponent
 296
 297          begin
 298             Mantissa := Random (Gen) / 2**Extra_Bits;
 299             R := Unsigned_32 (Mantissa mod 2**Extra_Bits);
 300             R_Bits := Extra_Bits;
 301             X := Real (2**(T'Machine_Mantissa - 1) + Mantissa); -- Exact
 302
 303             if Extra_Bits < 4 and then R < 2 ** Extra_Bits - 1 then
 304
 305                --  We got lucky and got a zero in our few extra bits
 306
 307                K := Trailing_Ones (R);
 308
 309             else
 310                Find_Zero : loop
 311
 312                   --  R has R_Bits unprocessed random bits, a multiple of 4.
 313                   --  X needs to be halved for each trailing one bit. The
 314                   --  process stops as soon as a 0 bit is found. If R_Bits
 315                   --  becomes zero, reload R.
 316
 317                   --  Process 4 bits at a time for speed: the two iterations
 318                   --  on average with three tests each was still too slow,
 319                   --  probably because the branches are not predictable.
 320                   --  This loop now will only execute once 94% of the cases,
 321                   --  doing more bits at a time will not help.
 322
 323                   while R_Bits >= 4 loop
 324                      K := Trailing_Ones (R mod 16);
 325
 326                      exit Find_Zero when K < 4; -- Exits 94% of the time
 327
 328                      R_Bits := R_Bits - 4;
 329                      X := X / 16.0;
 330                      R := R / 16;
 331                   end loop;
 332
 333                   --  Do not allow us to loop endlessly even in the (very
 334                   --  unlikely) case that Random (Gen) keeps yielding all ones.
 335
 336                   exit Find_Zero when X = 0.0;
 337                   R := Random (Gen);
 338                   R_Bits := 32;
 339                end loop Find_Zero;
 340             end if;
 341
 342             --  K has the count of trailing ones not reflected yet in X. The
 343             --  following multiplication takes care of that, as well as the
 344             --  correction to move the radix point to the left of the mantissa.
 345             --  Doing it at the end avoids repeated rounding errors in the
 346             --  exceedingly unlikely case of ever having a subnormal result.
 347
 348             X := X * Pow_Tab (K);
 349
 350             --  The smallest value in each binade is rounded to by 0.75 of
 351             --  the span of real numbers as its next larger neighbor, and
 352             --  1.0 is rounded to by half of the span of real numbers as its
 353             --  next smaller neighbor. To account for this, when we encounter
 354             --  the smallest number in a binade, we substitute the smallest
 355             --  value in the next larger binade with probability 1/2.
 356
 357             if Mantissa = 0 and then Unsigned_32'(Random (Gen)) mod 2 = 0 then
 358                X := 2.0 * X;
 359             end if;
 360
 361             return X;
 362          end;
 363       end if;
 364    end Random_Float_Template;
 365
 366    ------------
 367    -- Random --
 368    ------------
 369
 370    function Random (Gen : Generator) return Float is
 371       function F is new Random_Float_Template (Unsigned_32, Float);
 372    begin
 373       return F (Gen);
 374    end Random;
 375
 376    function Random (Gen : Generator) return Long_Float is
 377       function F is new Random_Float_Template (Unsigned_64, Long_Float);
 378    begin
 379       return F (Gen);
 380    end Random;
 381
 382    function Random (Gen : Generator) return Unsigned_64 is
 383    begin
 384       return Shift_Left (Unsigned_64 (Unsigned_32'(Random (Gen))), 32)
 385         or Unsigned_64 (Unsigned_32'(Random (Gen)));
 386    end Random;
 387
 388    ---------------------
 389    -- Random_Discrete --
 390    ---------------------
 391
 392    function Random_Discrete
 393      (Gen : Generator;
 394       Min : Result_Subtype := Default_Min;
 395       Max : Result_Subtype := Result_Subtype'Last) return Result_Subtype
 396    is
 397    begin
 398       if Max = Min then
 399          return Max;
 400
 401       elsif Max < Min then
 402          raise Constraint_Error;
 403
 404       elsif Result_Subtype'Base'Size > 32 then
 405          declare
 406             --  In the 64-bit case, we have to be careful, since not all 64-bit
 407             --  unsigned values are representable in GNAT's root_integer type.
 408             --  Ignore different-size warnings here since GNAT's handling
 409             --  is correct.
 410
 411             pragma Warnings ("Z");
 412             function Conv_To_Unsigned is
 413                new Unchecked_Conversion (Result_Subtype'Base, Unsigned_64);
 414             function Conv_To_Result is
 415                new Unchecked_Conversion (Unsigned_64, Result_Subtype'Base);
 416             pragma Warnings ("z");
 417
 418             N : constant Unsigned_64 :=
 419                   Conv_To_Unsigned (Max) - Conv_To_Unsigned (Min) + 1;
 420
 421             X, Slop : Unsigned_64;
 422
 423          begin
 424             if N = 0 then
 425                return Conv_To_Result (Conv_To_Unsigned (Min) + Random (Gen));
 426
 427             else
 428                Slop := Unsigned_64'Last rem N + 1;
 429
 430                loop
 431                   X := Random (Gen);
 432                   exit when Slop = N or else X <= Unsigned_64'Last - Slop;
 433                end loop;
 434
 435                return Conv_To_Result (Conv_To_Unsigned (Min) + X rem N);
 436             end if;
 437          end;
 438
 439       elsif Result_Subtype'Pos (Max) - Result_Subtype'Pos (Min) =
 440                                                          2 ** 32 - 1
 441       then
 442          return Result_Subtype'Val
 443            (Result_Subtype'Pos (Min) + Unsigned_32'Pos (Random (Gen)));
 444       else
 445          declare
 446             N    : constant Unsigned_32 :=
 447                      Unsigned_32 (Result_Subtype'Pos (Max) -
 448                                     Result_Subtype'Pos (Min) + 1);
 449             Slop : constant Unsigned_32 := Unsigned_32'Last rem N + 1;
 450             X    : Unsigned_32;
 451
 452          begin
 453             loop
 454                X := Random (Gen);
 455                exit when Slop = N or else X <= Unsigned_32'Last - Slop;
 456             end loop;
 457
 458             return
 459               Result_Subtype'Val
 460                 (Result_Subtype'Pos (Min) + Unsigned_32'Pos (X rem N));
 461          end;
 462       end if;
 463    end Random_Discrete;
 464
 465    ------------------
 466    -- Random_Float --
 467    ------------------
 468
 469    function Random_Float (Gen : Generator) return Result_Subtype is
 470    begin
 471       if Result_Subtype'Base'Digits > Float'Digits then
 472          return Result_Subtype'Machine (Result_Subtype
 473                                          (Long_Float'(Random (Gen))));
 474       else
 475          return Result_Subtype'Machine (Result_Subtype
 476                                          (Float'(Random (Gen))));
 477       end if;
 478    end Random_Float;
 479
 480    -----------
 481    -- Reset --
 482    -----------
 483
 484    procedure Reset (Gen : Generator) is
 485    begin
 486       Init (Gen, Unsigned_32'Mod (Random_Seed.Get_Seed));
 487    end Reset;
 488
 489    procedure Reset (Gen : Generator; Initiator : Integer_32) is
 490    begin
 491       Init (Gen, To_Unsigned (Initiator));
 492    end Reset;
 493
 494    procedure Reset (Gen : Generator; Initiator : Unsigned_32) is
 495    begin
 496       Init (Gen, Initiator);
 497    end Reset;
 498
 499    procedure Reset (Gen : Generator; Initiator : Integer) is
 500    begin
 501       --  This is probably an unnecessary precaution against future change, but
 502       --  since the test is a static expression, no extra code is involved.
 503
 504       if Integer'Size <= 32 then
 505          Init (Gen, To_Unsigned (Integer_32 (Initiator)));
 506
 507       else
 508          declare
 509             Initiator1 : constant Unsigned_64 :=
 510                            To_Unsigned (Integer_64 (Initiator));
 511             Init0      : constant Unsigned_32 :=
 512                            Unsigned_32 (Initiator1 mod 2 ** 32);
 513             Init1      : constant Unsigned_32 :=
 514                            Unsigned_32 (Shift_Right (Initiator1, 32));
 515          begin
 516             Reset (Gen, Initialization_Vector'(Init0, Init1));
 517          end;
 518       end if;
 519    end Reset;
 520
 521    procedure Reset (Gen : Generator; Initiator : Initialization_Vector) is
 522       G    : Generator renames Gen.Writable.Self.all;
 523       I, J : Integer;
 524
 525    begin
 526       Init (G, Seed1);
 527       I := 1;
 528       J := 0;
 529
 530       if Initiator'Length > 0 then
 531          for K in reverse 1 .. Integer'Max (N, Initiator'Length) loop
 532             G.S (I) :=
 533               (G.S (I) xor ((G.S (I - 1)
 534                                xor Shift_Right (G.S (I - 1), 30)) * Mult1))
 535               + Initiator (J + Initiator'First) + Unsigned_32 (J);
 536
 537             I := I + 1;
 538             J := J + 1;
 539
 540             if I >= N then
 541                G.S (0) := G.S (N - 1);
 542                I := 1;
 543             end if;
 544
 545             if J >= Initiator'Length then
 546                J := 0;
 547             end if;
 548          end loop;
 549       end if;
 550
 551       for K in reverse 1 .. N - 1 loop
 552          G.S (I) :=
 553            (G.S (I) xor ((G.S (I - 1)
 554                             xor Shift_Right (G.S (I - 1), 30)) * Mult2))
 555            - Unsigned_32 (I);
 556          I := I + 1;
 557
 558          if I >= N then
 559             G.S (0) := G.S (N - 1);
 560             I := 1;
 561          end if;
 562       end loop;
 563
 564       G.S (0) := Upper_Mask;
 565    end Reset;
 566
 567    procedure Reset (Gen : Generator; From_State : Generator) is
 568       G : Generator renames Gen.Writable.Self.all;
 569    begin
 570       G.S := From_State.S;
 571       G.I := From_State.I;
 572    end Reset;
 573
 574    procedure Reset (Gen : Generator; From_State : State) is
 575       G : Generator renames Gen.Writable.Self.all;
 576    begin
 577       G.I := 0;
 578       G.S := From_State;
 579    end Reset;
 580
 581    procedure Reset (Gen : Generator; From_Image : String) is
 582       G : Generator renames Gen.Writable.Self.all;
 583    begin
 584       G.I := 0;
 585
 586       for J in 0 .. N - 1 loop
 587          G.S (J) := Extract_Value (From_Image, J);
 588       end loop;
 589    end Reset;
 590
 591    ----------
 592    -- Save --
 593    ----------
 594
 595    procedure Save (Gen : Generator; To_State : out State) is
 596       Gen2 : Generator;
 597
 598    begin
 599       if Gen.I = N then
 600          Init (Gen2, 5489);
 601          To_State := Gen2.S;
 602
 603       else
 604          To_State (0 .. N - 1 - Gen.I) := Gen.S (Gen.I .. N - 1);
 605          To_State (N - Gen.I .. N - 1) := Gen.S (0 .. Gen.I - 1);
 606       end if;
 607    end Save;
 608
 609    -----------
 610    -- Image --
 611    -----------
 612
 613    function Image (Of_State : State) return String is
 614       Result : Image_String;
 615
 616    begin
 617       Result := (others => ' ');
 618
 619       for J in Of_State'Range loop
 620          Insert_Image (Result, J, Of_State (J));
 621       end loop;
 622
 623       return Result;
 624    end Image;
 625
 626    function Image (Gen : Generator) return String is
 627       Result : Image_String;
 628
 629    begin
 630       Result := (others => ' ');
 631       for J in 0 .. N - 1 loop
 632          Insert_Image (Result, J, Gen.S ((J + Gen.I) mod N));
 633       end loop;
 634
 635       return Result;
 636    end Image;
 637
 638    -----------
 639    -- Value --
 640    -----------
 641
 642    function Value (Coded_State : String) return State is
 643       Gen : Generator;
 644       S   : State;
 645    begin
 646       Reset (Gen, Coded_State);
 647       Save (Gen, S);
 648       return S;
 649    end Value;
 650
 651    ----------
 652    -- Init --
 653    ----------
 654
 655    procedure Init (Gen : Generator; Initiator : Unsigned_32) is
 656       G : Generator renames Gen.Writable.Self.all;
 657    begin
 658       G.S (0) := Initiator;
 659
 660       for I in 1 .. N - 1 loop
 661          G.S (I) :=
 662            (G.S (I - 1) xor Shift_Right (G.S (I - 1), 30)) * Mult0
 663            + Unsigned_32 (I);
 664       end loop;
 665
 666       G.I := 0;
 667    end Init;
 668
 669    ------------------
 670    -- Insert_Image --
 671    ------------------
 672
 673    procedure Insert_Image
 674      (S     : in out Image_String;
 675       Index : Integer;
 676       V     : State_Val)
 677    is
 678       Value : constant String := State_Val'Image (V);
 679    begin
 680       S (Index * 11 + 1 .. Index * 11 + Value'Length) := Value;
 681    end Insert_Image;
 682
 683    -------------------
 684    -- Extract_Value --
 685    -------------------
 686
 687    function Extract_Value (S : String; Index : Integer) return State_Val is
 688       Start : constant Integer := S'First + Index * Image_Numeral_Length;
 689    begin
 690       return State_Val'Value (S (Start .. Start + Image_Numeral_Length - 1));
 691    end Extract_Value;
 692
 693 end System.Random_Numbers;