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