aboutsummaryrefslogtreecommitdiffstats
path: root/gcc-4.6/gcc/ada/s-rannum.adb
blob: d85dd2efacf1a134fb6e876ceb16a41d3dcbea03 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
------------------------------------------------------------------------------
--                                                                          --
--                         GNAT RUN-TIME COMPONENTS                         --
--                                                                          --
--                S Y S T E M . R A N D O M _ N U M B E R S                 --
--                                                                          --
--                                 B o d y                                  --
--                                                                          --
--          Copyright (C) 2007-2010, Free Software Foundation, Inc.         --
--                                                                          --
-- GNAT is free software;  you can  redistribute it  and/or modify it under --
-- terms of the  GNU General Public License as published  by the Free Soft- --
-- ware  Foundation;  either version 3,  or (at your option) any later ver- --
-- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
-- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
-- or FITNESS FOR A PARTICULAR PURPOSE.                                     --
--                                                                          --
-- As a special exception under Section 7 of GPL version 3, you are granted --
-- additional permissions described in the GCC Runtime Library Exception,   --
-- version 3.1, as published by the Free Software Foundation.               --
--                                                                          --
-- You should have received a copy of the GNU General Public License and    --
-- a copy of the GCC Runtime Library Exception along with this program;     --
-- see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see    --
-- <http://www.gnu.org/licenses/>.                                          --
--                                                                          --
-- GNAT was originally developed  by the GNAT team at  New York University. --
-- Extensive contributions were provided by Ada Core Technologies Inc.      --
--                                                                          --
------------------------------------------------------------------------------

------------------------------------------------------------------------------
--                                                                          --
-- The implementation here is derived from a C-program for MT19937, with    --
-- initialization improved 2002/1/26. As required, the following notice is  --
-- copied from the original program.                                        --
--                                                                          --
-- Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,        --
-- All rights reserved.                                                     --
--                                                                          --
-- Redistribution and use in source and binary forms, with or without       --
-- modification, are permitted provided that the following conditions       --
-- are met:                                                                 --
--                                                                          --
--   1. Redistributions of source code must retain the above copyright      --
--      notice, this list of conditions and the following disclaimer.       --
--                                                                          --
--   2. Redistributions in binary form must reproduce the above copyright   --
--      notice, this list of conditions and the following disclaimer in the --
--      documentation and/or other materials provided with the distribution.--
--                                                                          --
--   3. The names of its contributors may not be used to endorse or promote --
--      products derived from this software without specific prior written  --
--      permission.                                                         --
--                                                                          --
-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS      --
-- "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT        --
-- LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR    --
-- A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT    --
-- OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,    --
-- SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED --
-- TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR   --
-- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF   --
-- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING     --
-- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS       --
-- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.             --
--                                                                          --
------------------------------------------------------------------------------

------------------------------------------------------------------------------
--                                                                          --
-- This is an implementation of the Mersenne Twister, twisted generalized   --
-- feedback shift register of rational normal form, with state-bit          --
-- reflection and tempering. This version generates 32-bit integers with a  --
-- period of 2**19937 - 1 (a Mersenne prime, hence the name). For           --
-- applications requiring more than 32 bits (up to 64), we concatenate two  --
-- 32-bit numbers.                                                          --
--                                                                          --
-- See http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html for         --
-- details.                                                                 --
--                                                                          --
-- In contrast to the original code, we do not generate random numbers in   --
-- batches of N. Measurement seems to show this has very little if any      --
-- effect on performance, and it may be marginally better for real-time     --
-- applications with hard deadlines.                                        --
--                                                                          --
------------------------------------------------------------------------------

with Ada.Calendar;             use Ada.Calendar;
with Ada.Unchecked_Conversion;

with Interfaces; use Interfaces;

use Ada;

package body System.Random_Numbers is

   Y2K : constant Calendar.Time :=
           Calendar.Time_Of
             (Year => 2000, Month => 1, Day => 1, Seconds => 0.0);
   --  First day of Year 2000 (what is this for???)

   Image_Numeral_Length : constant := Max_Image_Width / N;
   subtype Image_String is String (1 .. Max_Image_Width);

   ----------------------------
   -- Algorithmic Parameters --
   ----------------------------

   Lower_Mask : constant := 2**31-1;
   Upper_Mask : constant := 2**31;

   Matrix_A   : constant array (State_Val range 0 .. 1) of State_Val
     := (0, 16#9908b0df#);
   --  The twist transformation is represented by a matrix of the form
   --
   --               [  0    I(31) ]
   --               [    _a       ]
   --
   --  where 0 is a 31x31 block of 0s, I(31) is the 31x31 identity matrix and
   --  _a is a particular bit row-vector, represented here by a 32-bit integer.
   --  If integer x represents a row vector of bits (with x(0), the units bit,
   --  last), then
   --           x * A = [0 x(31..1)] xor Matrix_A(x(0)).

   U      : constant := 11;
   S      : constant := 7;
   B_Mask : constant := 16#9d2c5680#;
   T      : constant := 15;
   C_Mask : constant := 16#efc60000#;
   L      : constant := 18;
   --  The tempering shifts and bit masks, in the order applied

   Seed0 : constant := 5489;
   --  Default seed, used to initialize the state vector when Reset not called

   Seed1 : constant := 19650218;
   --  Seed used to initialize the state vector when calling Reset with an
   --  initialization vector.

   Mult0 : constant := 1812433253;
   --  Multiplier for a modified linear congruential generator used to
   --  initialize the state vector when calling Reset with a single integer
   --  seed.

   Mult1 : constant := 1664525;
   Mult2 : constant := 1566083941;
   --  Multipliers for two modified linear congruential generators used to
   --  initialize the state vector when calling Reset with an initialization
   --  vector.

   -----------------------
   -- Local Subprograms --
   -----------------------

   procedure Init (Gen : Generator; Initiator : Unsigned_32);
   --  Perform a default initialization of the state of Gen. The resulting
   --  state is identical for identical values of Initiator.

   procedure Insert_Image
     (S     : in out Image_String;
      Index : Integer;
      V     : State_Val);
   --  Insert image of V into S, in the Index'th 11-character substring

   function Extract_Value (S : String; Index : Integer) return State_Val;
   --  Treat S as a sequence of 11-character decimal numerals and return
   --  the result of converting numeral #Index (numbering from 0)

   function To_Unsigned is
     new Unchecked_Conversion (Integer_32, Unsigned_32);
   function To_Unsigned is
     new Unchecked_Conversion (Integer_64, Unsigned_64);

   ------------
   -- Random --
   ------------

   function Random (Gen : Generator) return Unsigned_32 is
      G : Generator renames Gen.Writable.Self.all;
      Y : State_Val;
      I : Integer;      --  should avoid use of identifier I ???

   begin
      I := G.I;

      if I < N - M then
         Y := (G.S (I) and Upper_Mask) or (G.S (I + 1) and Lower_Mask);
         Y := G.S (I + M) xor Shift_Right (Y, 1) xor Matrix_A (Y and 1);
         I := I + 1;

      elsif I < N - 1 then
         Y := (G.S (I) and Upper_Mask) or (G.S (I + 1) and Lower_Mask);
         Y := G.S (I + (M - N))
                xor Shift_Right (Y, 1)
                xor Matrix_A (Y and 1);
         I := I + 1;

      elsif I = N - 1 then
         Y := (G.S (I) and Upper_Mask) or (G.S (0) and Lower_Mask);
         Y := G.S (M - 1) xor Shift_Right (Y, 1) xor Matrix_A (Y and 1);
         I := 0;

      else
         Init (G, Seed0);
         return Random (Gen);
      end if;

      G.S (G.I) := Y;
      G.I := I;

      Y := Y xor Shift_Right (Y, U);
      Y := Y xor (Shift_Left (Y, S)  and B_Mask);
      Y := Y xor (Shift_Left (Y, T) and C_Mask);
      Y := Y xor Shift_Right (Y, L);

      return Y;
   end Random;

   generic
      type Unsigned is mod <>;
      type Real is digits <>;
      with function Random (G : Generator) return Unsigned is <>;
   function Random_Float_Template (Gen : Generator) return Real;
   pragma Inline (Random_Float_Template);
   --  Template for a random-number generator implementation that delivers
   --  values of type Real in the range [0 .. 1], using values from Gen,
   --  assuming that Unsigned is large enough to hold the bits of a mantissa
   --  for type Real.

   ---------------------------
   -- Random_Float_Template --
   ---------------------------

   function Random_Float_Template (Gen : Generator) return Real is

      pragma Compile_Time_Error
        (Unsigned'Last <= 2**(Real'Machine_Mantissa - 1),
         "insufficiently large modular type used to hold mantissa");

   begin
      --  This code generates random floating-point numbers from unsigned
      --  integers. Assuming that Real'Machine_Radix = 2, it can deliver all
      --  machine values of type Real (as implied by Real'Machine_Mantissa and
      --  Real'Machine_Emin), which is not true of the standard method (to
      --  which we fall back for non-binary radix): computing Real(<random
      --  integer>) / (<max random integer>+1). To do so, we first extract an
      --  (M-1)-bit significand (where M is Real'Machine_Mantissa), and then
      --  decide on a normalized exponent by repeated coin flips, decrementing
      --  from 0 as long as we flip heads (1 bits). This process yields the
      --  proper geometric distribution for the exponent: in a uniformly
      --  distributed set of floating-point numbers, 1/2 of them will be in
      --  (0.5, 1], 1/4 will be in (0.25, 0.5], and so forth. It makes a
      --  further adjustment at binade boundaries (see comments below) to give
      --  the effect of selecting a uniformly distributed real deviate in
      --  [0..1] and then rounding to the nearest representable floating-point
      --  number.  The algorithm attempts to be stingy with random integers. In
      --  the worst case, it can consume roughly -Real'Machine_Emin/32 32-bit
      --  integers, but this case occurs with probability around
      --  2**Machine_Emin, and the expected number of calls to integer-valued
      --  Random is 1.  For another discussion of the issues addressed by this
      --  process, see Allen Downey's unpublished paper at
      --  http://allendowney.com/research/rand/downey07randfloat.pdf.

      if Real'Machine_Radix /= 2 then
         return Real'Machine
           (Real (Unsigned'(Random (Gen))) * 2.0**(-Unsigned'Size));

      else
         declare
            type Bit_Count is range 0 .. 4;

            subtype T is Real'Base;

            Trailing_Ones : constant array (Unsigned_32 range 0 .. 15)
              of Bit_Count :=
                  (2#00000# => 0, 2#00001# => 1, 2#00010# => 0, 2#00011# => 2,
                   2#00100# => 0, 2#00101# => 1, 2#00110# => 0, 2#00111# => 3,
                   2#01000# => 0, 2#01001# => 1, 2#01010# => 0, 2#01011# => 2,
                   2#01100# => 0, 2#01101# => 1, 2#01110# => 0, 2#01111# => 4);

            Pow_Tab : constant array (Bit_Count range 0 .. 3) of Real
              := (0 => 2.0**(0 - T'Machine_Mantissa),
                  1 => 2.0**(-1 - T'Machine_Mantissa),
                  2 => 2.0**(-2 - T'Machine_Mantissa),
                  3 => 2.0**(-3 - T'Machine_Mantissa));

            Extra_Bits : constant Natural :=
                         (Unsigned'Size - T'Machine_Mantissa + 1);
            --  Random bits left over after selecting mantissa

            Mantissa : Unsigned;

            X      : Real;            --  Scaled mantissa
            R      : Unsigned_32;     --  Supply of random bits
            R_Bits : Natural;         --  Number of bits left in R
            K      : Bit_Count;       --  Next decrement to exponent

         begin
            Mantissa := Random (Gen) / 2**Extra_Bits;
            R := Unsigned_32 (Mantissa mod 2**Extra_Bits);
            R_Bits := Extra_Bits;
            X := Real (2**(T'Machine_Mantissa - 1) + Mantissa); -- Exact

            if Extra_Bits < 4 and then R < 2 ** Extra_Bits - 1 then

               --  We got lucky and got a zero in our few extra bits

               K := Trailing_Ones (R);

            else
               Find_Zero : loop

                  --  R has R_Bits unprocessed random bits, a multiple of 4.
                  --  X needs to be halved for each trailing one bit. The
                  --  process stops as soon as a 0 bit is found. If R_Bits
                  --  becomes zero, reload R.

                  --  Process 4 bits at a time for speed: the two iterations
                  --  on average with three tests each was still too slow,
                  --  probably because the branches are not predictable.
                  --  This loop now will only execute once 94% of the cases,
                  --  doing more bits at a time will not help.

                  while R_Bits >= 4 loop
                     K := Trailing_Ones (R mod 16);

                     exit Find_Zero when K < 4; -- Exits 94% of the time

                     R_Bits := R_Bits - 4;
                     X := X / 16.0;
                     R := R / 16;
                  end loop;

                  --  Do not allow us to loop endlessly even in the (very
                  --  unlikely) case that Random (Gen) keeps yielding all ones.

                  exit Find_Zero when X = 0.0;
                  R := Random (Gen);
                  R_Bits := 32;
               end loop Find_Zero;
            end if;

            --  K has the count of trailing ones not reflected yet in X. The
            --  following multiplication takes care of that, as well as the
            --  correction to move the radix point to the left of the mantissa.
            --  Doing it at the end avoids repeated rounding errors in the
            --  exceedingly unlikely case of ever having a subnormal result.

            X := X * Pow_Tab (K);

            --  The smallest value in each binade is rounded to by 0.75 of
            --  the span of real numbers as its next larger neighbor, and
            --  1.0 is rounded to by half of the span of real numbers as its
            --  next smaller neighbor. To account for this, when we encounter
            --  the smallest number in a binade, we substitute the smallest
            --  value in the next larger binade with probability 1/2.

            if Mantissa = 0 and then Unsigned_32'(Random (Gen)) mod 2 = 0 then
               X := 2.0 * X;
            end if;

            return X;
         end;
      end if;
   end Random_Float_Template;

   ------------
   -- Random --
   ------------

   function Random (Gen : Generator) return Float is
      function F is new Random_Float_Template (Unsigned_32, Float);
   begin
      return F (Gen);
   end Random;

   function Random (Gen : Generator) return Long_Float is
      function F is new Random_Float_Template (Unsigned_64, Long_Float);
   begin
      return F (Gen);
   end Random;

   function Random (Gen : Generator) return Unsigned_64 is
   begin
      return Shift_Left (Unsigned_64 (Unsigned_32'(Random (Gen))), 32)
        or Unsigned_64 (Unsigned_32'(Random (Gen)));
   end Random;

   ---------------------
   -- Random_Discrete --
   ---------------------

   function Random_Discrete
     (Gen : Generator;
      Min : Result_Subtype := Default_Min;
      Max : Result_Subtype := Result_Subtype'Last) return Result_Subtype
   is
   begin
      if Max = Min then
         return Max;

      elsif Max < Min then
         raise Constraint_Error;

      elsif Result_Subtype'Base'Size > 32 then
         declare
            --  In the 64-bit case, we have to be careful, since not all 64-bit
            --  unsigned values are representable in GNAT's root_integer type.
            --  Ignore different-size warnings here since GNAT's handling
            --  is correct.

            pragma Warnings ("Z");  -- better to use msg string! ???
            function Conv_To_Unsigned is
               new Unchecked_Conversion (Result_Subtype'Base, Unsigned_64);
            function Conv_To_Result is
               new Unchecked_Conversion (Unsigned_64, Result_Subtype'Base);
            pragma Warnings ("z");

            N : constant Unsigned_64 :=
                  Conv_To_Unsigned (Max) - Conv_To_Unsigned (Min) + 1;

            X, Slop : Unsigned_64;

         begin
            if N = 0 then
               return Conv_To_Result (Conv_To_Unsigned (Min) + Random (Gen));

            else
               Slop := Unsigned_64'Last rem N + 1;

               loop
                  X := Random (Gen);
                  exit when Slop = N or else X <= Unsigned_64'Last - Slop;
               end loop;

               return Conv_To_Result (Conv_To_Unsigned (Min) + X rem N);
            end if;
         end;

      elsif Result_Subtype'Pos (Max) - Result_Subtype'Pos (Min) =
                                                         2 ** 32 - 1
      then
         return Result_Subtype'Val
           (Result_Subtype'Pos (Min) + Unsigned_32'Pos (Random (Gen)));
      else
         declare
            N    : constant Unsigned_32 :=
                     Unsigned_32 (Result_Subtype'Pos (Max) -
                                    Result_Subtype'Pos (Min) + 1);
            Slop : constant Unsigned_32 := Unsigned_32'Last rem N + 1;
            X    : Unsigned_32;

         begin
            loop
               X := Random (Gen);
               exit when Slop = N or else X <= Unsigned_32'Last - Slop;
            end loop;

            return
              Result_Subtype'Val
                (Result_Subtype'Pos (Min) + Unsigned_32'Pos (X rem N));
         end;
      end if;
   end Random_Discrete;

   ------------------
   -- Random_Float --
   ------------------

   function Random_Float (Gen : Generator) return Result_Subtype is
   begin
      if Result_Subtype'Base'Digits > Float'Digits then
         return Result_Subtype'Machine (Result_Subtype
                                         (Long_Float'(Random (Gen))));
      else
         return Result_Subtype'Machine (Result_Subtype
                                         (Float'(Random (Gen))));
      end if;
   end Random_Float;

   -----------
   -- Reset --
   -----------

   procedure Reset (Gen : Generator) is
      Clock              : constant Time := Calendar.Clock;
      Duration_Since_Y2K : constant Duration := Clock - Y2K;

      X : constant Unsigned_32 :=
            Unsigned_32'Mod (Unsigned_64 (Duration_Since_Y2K) * 64);

   begin
      Init (Gen, X);
   end Reset;

   procedure Reset (Gen : Generator; Initiator : Integer_32) is
   begin
      Init (Gen, To_Unsigned (Initiator));
   end Reset;

   procedure Reset (Gen : Generator; Initiator : Unsigned_32) is
   begin
      Init (Gen, Initiator);
   end Reset;

   procedure Reset (Gen : Generator; Initiator : Integer) is
   begin
      pragma Warnings (Off, "condition is always *");
      --  This is probably an unnecessary precaution against future change, but
      --  since the test is a static expression, no extra code is involved.

      if Integer'Size <= 32 then
         Init (Gen, To_Unsigned (Integer_32 (Initiator)));

      else
         declare
            Initiator1 : constant Unsigned_64 :=
                           To_Unsigned (Integer_64 (Initiator));
            Init0      : constant Unsigned_32 :=
                           Unsigned_32 (Initiator1 mod 2 ** 32);
            Init1      : constant Unsigned_32 :=
                           Unsigned_32 (Shift_Right (Initiator1, 32));
         begin
            Reset (Gen, Initialization_Vector'(Init0, Init1));
         end;
      end if;

      pragma Warnings (On, "condition is always *");
   end Reset;

   procedure Reset (Gen : Generator; Initiator : Initialization_Vector) is
      G    : Generator renames Gen.Writable.Self.all;
      I, J : Integer;

   begin
      Init (G, Seed1);
      I := 1;
      J := 0;

      if Initiator'Length > 0 then
         for K in reverse 1 .. Integer'Max (N, Initiator'Length) loop
            G.S (I) :=
              (G.S (I) xor ((G.S (I - 1)
                               xor Shift_Right (G.S (I - 1), 30)) * Mult1))
              + Initiator (J + Initiator'First) + Unsigned_32 (J);

            I := I + 1;
            J := J + 1;

            if I >= N then
               G.S (0) := G.S (N - 1);
               I := 1;
            end if;

            if J >= Initiator'Length then
               J := 0;
            end if;
         end loop;
      end if;

      for K in reverse 1 .. N - 1 loop
         G.S (I) :=
           (G.S (I) xor ((G.S (I - 1)
                            xor Shift_Right (G.S (I - 1), 30)) * Mult2))
           - Unsigned_32 (I);
         I := I + 1;

         if I >= N then
            G.S (0) := G.S (N - 1);
            I := 1;
         end if;
      end loop;

      G.S (0) := Upper_Mask;
   end Reset;

   procedure Reset (Gen : Generator; From_State : Generator) is
      G : Generator renames Gen.Writable.Self.all;
   begin
      G.S := From_State.S;
      G.I := From_State.I;
   end Reset;

   procedure Reset (Gen : Generator; From_State : State) is
      G : Generator renames Gen.Writable.Self.all;
   begin
      G.I := 0;
      G.S := From_State;
   end Reset;

   procedure Reset (Gen : Generator; From_Image : String) is
      G : Generator renames Gen.Writable.Self.all;
   begin
      G.I := 0;

      for J in 0 .. N - 1 loop
         G.S (J) := Extract_Value (From_Image, J);
      end loop;
   end Reset;

   ----------
   -- Save --
   ----------

   procedure Save (Gen : Generator; To_State : out State) is
      Gen2 : Generator;

   begin
      if Gen.I = N then
         Init (Gen2, 5489);
         To_State := Gen2.S;

      else
         To_State (0 .. N - 1 - Gen.I) := Gen.S (Gen.I .. N - 1);
         To_State (N - Gen.I .. N - 1) := Gen.S (0 .. Gen.I - 1);
      end if;
   end Save;

   -----------
   -- Image --
   -----------

   function Image (Of_State : State) return String is
      Result : Image_String;

   begin
      Result := (others => ' ');

      for J in Of_State'Range loop
         Insert_Image (Result, J, Of_State (J));
      end loop;

      return Result;
   end Image;

   function Image (Gen : Generator) return String is
      Result : Image_String;

   begin
      Result := (others => ' ');
      for J in 0 .. N - 1 loop
         Insert_Image (Result, J, Gen.S ((J + Gen.I) mod N));
      end loop;

      return Result;
   end Image;

   -----------
   -- Value --
   -----------

   function Value (Coded_State : String) return State is
      Gen : Generator;
      S   : State;
   begin
      Reset (Gen, Coded_State);
      Save (Gen, S);
      return S;
   end Value;

   ----------
   -- Init --
   ----------

   procedure Init (Gen : Generator; Initiator : Unsigned_32) is
      G : Generator renames Gen.Writable.Self.all;
   begin
      G.S (0) := Initiator;

      for I in 1 .. N - 1 loop
         G.S (I) :=
           (G.S (I - 1) xor Shift_Right (G.S (I - 1), 30)) * Mult0
           + Unsigned_32 (I);
      end loop;

      G.I := 0;
   end Init;

   ------------------
   -- Insert_Image --
   ------------------

   procedure Insert_Image
     (S     : in out Image_String;
      Index : Integer;
      V     : State_Val)
   is
      Value : constant String := State_Val'Image (V);
   begin
      S (Index * 11 + 1 .. Index * 11 + Value'Length) := Value;
   end Insert_Image;

   -------------------
   -- Extract_Value --
   -------------------

   function Extract_Value (S : String; Index : Integer) return State_Val is
      Start : constant Integer := S'First + Index * Image_Numeral_Length;
   begin
      return State_Val'Value (S (Start .. Start + Image_Numeral_Length - 1));
   end Extract_Value;
end System.Random_Numbers;