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