1------------------------------------------------------------------------------
2--                                                                          --
3--                         GNAT RUN-TIME COMPONENTS                         --
4--                                                                          --
5--                  G N A T . R A N D O M _ N U M B E R S                   --
6--                                                                          --
7--                                 B o d y                                  --
8--                                                                          --
9--      Copyright (C) 2007-2009  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
32with Ada.Numerics.Long_Elementary_Functions;
33use Ada.Numerics.Long_Elementary_Functions;
34with Ada.Unchecked_Conversion;
35with System.Random_Numbers; use System.Random_Numbers;
36
37package body GNAT.Random_Numbers is
38
39   Sys_Max_Image_Width : constant := System.Random_Numbers.Max_Image_Width;
40
41   subtype Image_String is String (1 .. Max_Image_Width);
42
43   --  Utility function declarations
44
45   procedure Insert_Image
46     (S     : in out Image_String;
47      Index : Integer;
48      V     : Integer_64);
49   --  Insert string representation of V in S starting at position Index
50
51   ---------------
52   -- To_Signed --
53   ---------------
54
55   function To_Signed is
56     new Ada.Unchecked_Conversion (Unsigned_32, Integer_32);
57   function To_Signed is
58     new Ada.Unchecked_Conversion (Unsigned_64, Integer_64);
59
60   ------------------
61   -- Insert_Image --
62   ------------------
63
64   procedure Insert_Image
65     (S     : in out Image_String;
66      Index : Integer;
67      V     : Integer_64)
68   is
69      Image : constant String := Integer_64'Image (V);
70   begin
71      S (Index .. Index + Image'Length - 1) := Image;
72   end Insert_Image;
73
74   ---------------------
75   -- Random_Discrete --
76   ---------------------
77
78   function Random_Discrete
79     (Gen   : Generator;
80      Min   : Result_Subtype := Default_Min;
81      Max   : Result_Subtype := Result_Subtype'Last) return Result_Subtype
82   is
83      function F is
84        new System.Random_Numbers.Random_Discrete
85              (Result_Subtype, Default_Min);
86   begin
87      return F (Gen.Rep, Min, Max);
88   end Random_Discrete;
89
90   ------------
91   -- Random --
92   ------------
93
94   function Random (Gen : Generator) return Float is
95   begin
96      return Random (Gen.Rep);
97   end Random;
98
99   function Random (Gen : Generator) return Long_Float is
100   begin
101      return Random (Gen.Rep);
102   end Random;
103
104   function Random (Gen : Generator) return Interfaces.Unsigned_32 is
105   begin
106      return Random (Gen.Rep);
107   end Random;
108
109   function Random (Gen : Generator) return Interfaces.Unsigned_64 is
110   begin
111      return Random (Gen.Rep);
112   end Random;
113
114   function Random (Gen : Generator) return Integer_64 is
115   begin
116      return To_Signed (Unsigned_64'(Random (Gen)));
117   end Random;
118
119   function Random (Gen : Generator) return Integer_32 is
120   begin
121      return To_Signed (Unsigned_32'(Random (Gen)));
122   end Random;
123
124   function Random (Gen : Generator) return Long_Integer is
125      function Random_Long_Integer is new Random_Discrete (Long_Integer);
126   begin
127      return Random_Long_Integer (Gen);
128   end Random;
129
130   function Random (Gen : Generator) return Integer is
131      function Random_Integer is new Random_Discrete (Integer);
132   begin
133      return Random_Integer (Gen);
134   end Random;
135
136   ------------------
137   -- Random_Float --
138   ------------------
139
140   function Random_Float (Gen   : Generator) return Result_Subtype is
141      function F is new System.Random_Numbers.Random_Float (Result_Subtype);
142   begin
143      return F (Gen.Rep);
144   end Random_Float;
145
146   ---------------------
147   -- Random_Gaussian --
148   ---------------------
149
150   --  Generates pairs of normally distributed values using the polar method of
151   --  G. E. P. Box, M. E. Muller, and G. Marsaglia. See Donald E. Knuth, The
152   --  Art of Computer Programming, Vol 2: Seminumerical Algorithms, section
153   --  3.4.1, subsection C, algorithm P. Returns half of the pair on each call,
154   --  using the Next_Gaussian field of Gen to hold the second member on
155   --  even-numbered calls.
156
157   function Random_Gaussian (Gen : Generator) return Long_Float is
158      G : Generator renames Gen'Unrestricted_Access.all;
159
160      V1, V2, Rad2, Mult : Long_Float;
161
162   begin
163      if G.Have_Gaussian then
164         G.Have_Gaussian := False;
165         return G.Next_Gaussian;
166
167      else
168         loop
169            V1 := 2.0 * Random (G) - 1.0;
170            V2 := 2.0 * Random (G) - 1.0;
171            Rad2 := V1 ** 2 + V2 ** 2;
172            exit when Rad2 < 1.0 and then Rad2 /= 0.0;
173         end loop;
174
175         --  Now V1 and V2 are coordinates in the unit circle
176
177         Mult := Sqrt (-2.0 * Log (Rad2) / Rad2);
178         G.Next_Gaussian := V2 * Mult;
179         G.Have_Gaussian := True;
180         return Long_Float'Machine (V1 * Mult);
181      end if;
182   end Random_Gaussian;
183
184   function Random_Gaussian (Gen : Generator) return Float is
185      V : constant Long_Float := Random_Gaussian (Gen);
186   begin
187      return Float'Machine (Float (V));
188   end Random_Gaussian;
189
190   -----------
191   -- Reset --
192   -----------
193
194   procedure Reset (Gen : out Generator) is
195   begin
196      Reset (Gen.Rep);
197      Gen.Have_Gaussian := False;
198   end Reset;
199
200   procedure Reset
201     (Gen       : out Generator;
202      Initiator : Initialization_Vector)
203   is
204   begin
205      Reset (Gen.Rep, Initiator);
206      Gen.Have_Gaussian := False;
207   end Reset;
208
209   procedure Reset
210     (Gen       : out Generator;
211      Initiator : Interfaces.Integer_32)
212   is
213   begin
214      Reset (Gen.Rep, Initiator);
215      Gen.Have_Gaussian := False;
216   end Reset;
217
218   procedure Reset
219     (Gen       : out Generator;
220      Initiator : Interfaces.Unsigned_32)
221   is
222   begin
223      Reset (Gen.Rep, Initiator);
224      Gen.Have_Gaussian := False;
225   end Reset;
226
227   procedure Reset
228     (Gen       : out Generator;
229      Initiator : Integer)
230   is
231   begin
232      Reset (Gen.Rep, Initiator);
233      Gen.Have_Gaussian := False;
234   end Reset;
235
236   procedure Reset
237     (Gen        : out Generator;
238      From_State : Generator)
239   is
240   begin
241      Reset (Gen.Rep, From_State.Rep);
242      Gen.Have_Gaussian := From_State.Have_Gaussian;
243      Gen.Next_Gaussian := From_State.Next_Gaussian;
244   end Reset;
245
246   Frac_Scale : constant Long_Float :=
247                  Long_Float
248                    (Long_Float'Machine_Radix) ** Long_Float'Machine_Mantissa;
249
250   function Val64 (Image : String) return Integer_64;
251   --  Renames Integer64'Value
252   --  We cannot use a 'renames Integer64'Value' since for some strange
253   --  reason, this requires a dependency on s-auxdec.ads which not all
254   --  run-times support ???
255
256   function Val64 (Image : String) return Integer_64 is
257   begin
258      return Integer_64'Value (Image);
259   end Val64;
260
261   procedure Reset
262     (Gen        : out Generator;
263      From_Image : String)
264   is
265      F0 : constant Integer := From_Image'First;
266      T0 : constant Integer := From_Image'First + Sys_Max_Image_Width;
267
268   begin
269      Reset (Gen.Rep, From_Image (F0 .. F0 + Sys_Max_Image_Width));
270
271      if From_Image (T0 + 1) = '1' then
272         Gen.Have_Gaussian := True;
273         Gen.Next_Gaussian :=
274           Long_Float (Val64 (From_Image (T0 + 3 .. T0 + 23))) / Frac_Scale
275           * Long_Float (Long_Float'Machine_Radix)
276           ** Integer (Val64 (From_Image (T0 + 25 .. From_Image'Last)));
277      else
278         Gen.Have_Gaussian := False;
279      end if;
280   end Reset;
281
282   -----------
283   -- Image --
284   -----------
285
286   function Image (Gen : Generator) return String is
287      Result : Image_String;
288
289   begin
290      Result := (others => ' ');
291      Result (1 .. Sys_Max_Image_Width) := Image (Gen.Rep);
292
293      if Gen.Have_Gaussian then
294         Result (Sys_Max_Image_Width + 2) := '1';
295         Insert_Image (Result, Sys_Max_Image_Width + 4,
296                       Integer_64 (Long_Float'Fraction (Gen.Next_Gaussian)
297                                   * Frac_Scale));
298         Insert_Image (Result, Sys_Max_Image_Width + 24,
299                       Integer_64 (Long_Float'Exponent (Gen.Next_Gaussian)));
300
301      else
302         Result (Sys_Max_Image_Width + 2) := '0';
303      end if;
304
305      return Result;
306   end Image;
307
308end GNAT.Random_Numbers;
309