|
111
|
1 ------------------------------------------------------------------------------
|
|
|
2 -- --
|
|
|
3 -- GNAT RUN-TIME COMPONENTS --
|
|
|
4 -- --
|
|
|
5 -- G N A T . M B B S _ F L O A T _ R A N D O M --
|
|
|
6 -- --
|
|
|
7 -- B o d y --
|
|
|
8 -- --
|
|
131
|
9 -- Copyright (C) 1992-2018, Free Software Foundation, Inc. --
|
|
111
|
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 with Ada.Calendar;
|
|
|
33
|
|
|
34 package body GNAT.MBBS_Float_Random is
|
|
|
35
|
|
|
36 -------------------------
|
|
|
37 -- Implementation Note --
|
|
|
38 -------------------------
|
|
|
39
|
|
|
40 -- The design of this spec is a bit awkward, as a result of Ada 95 not
|
|
|
41 -- permitting in-out parameters for function formals (most naturally
|
|
|
42 -- Generator values would be passed this way). In pure Ada 95, the only
|
|
|
43 -- solution would be to add a self-referential component to the generator
|
|
|
44 -- allowing access to the generator object from inside the function. This
|
|
|
45 -- would work because the generator is limited, which prevents any copy.
|
|
|
46
|
|
|
47 -- This is a bit heavy, so what we do is to use Unrestricted_Access to
|
|
|
48 -- get a pointer to the state in the passed Generator. This works because
|
|
|
49 -- Generator is a limited type and will thus always be passed by reference.
|
|
|
50
|
|
|
51 package Calendar renames Ada.Calendar;
|
|
|
52
|
|
|
53 type Pointer is access all State;
|
|
|
54
|
|
|
55 -----------------------
|
|
|
56 -- Local Subprograms --
|
|
|
57 -----------------------
|
|
|
58
|
|
|
59 procedure Euclid (P, Q : Int; X, Y : out Int; GCD : out Int);
|
|
|
60
|
|
|
61 function Euclid (P, Q : Int) return Int;
|
|
|
62
|
|
|
63 function Square_Mod_N (X, N : Int) return Int;
|
|
|
64
|
|
|
65 ------------
|
|
|
66 -- Euclid --
|
|
|
67 ------------
|
|
|
68
|
|
|
69 procedure Euclid (P, Q : Int; X, Y : out Int; GCD : out Int) is
|
|
|
70
|
|
|
71 XT : Int := 1;
|
|
|
72 YT : Int := 0;
|
|
|
73
|
|
|
74 procedure Recur
|
|
|
75 (P, Q : Int; -- a (i-1), a (i)
|
|
|
76 X, Y : Int; -- x (i), y (i)
|
|
|
77 XP, YP : in out Int; -- x (i-1), y (i-1)
|
|
|
78 GCD : out Int);
|
|
|
79
|
|
|
80 procedure Recur
|
|
|
81 (P, Q : Int;
|
|
|
82 X, Y : Int;
|
|
|
83 XP, YP : in out Int;
|
|
|
84 GCD : out Int)
|
|
|
85 is
|
|
|
86 Quo : Int := P / Q; -- q <-- |_ a (i-1) / a (i) _|
|
|
|
87 XT : Int := X; -- x (i)
|
|
|
88 YT : Int := Y; -- y (i)
|
|
|
89
|
|
|
90 begin
|
|
|
91 if P rem Q = 0 then -- while does not divide
|
|
|
92 GCD := Q;
|
|
|
93 XP := X;
|
|
|
94 YP := Y;
|
|
|
95 else
|
|
|
96 Recur (Q, P - Q * Quo, XP - Quo * X, YP - Quo * Y, XT, YT, Quo);
|
|
|
97
|
|
|
98 -- a (i) <== a (i)
|
|
|
99 -- a (i+1) <-- a (i-1) - q*a (i)
|
|
|
100 -- x (i+1) <-- x (i-1) - q*x (i)
|
|
|
101 -- y (i+1) <-- y (i-1) - q*y (i)
|
|
|
102 -- x (i) <== x (i)
|
|
|
103 -- y (i) <== y (i)
|
|
|
104
|
|
|
105 XP := XT;
|
|
|
106 YP := YT;
|
|
|
107 GCD := Quo;
|
|
|
108 end if;
|
|
|
109 end Recur;
|
|
|
110
|
|
|
111 -- Start of processing for Euclid
|
|
|
112
|
|
|
113 begin
|
|
|
114 Recur (P, Q, 0, 1, XT, YT, GCD);
|
|
|
115 X := XT;
|
|
|
116 Y := YT;
|
|
|
117 end Euclid;
|
|
|
118
|
|
|
119 function Euclid (P, Q : Int) return Int is
|
|
|
120 X, Y, GCD : Int;
|
|
|
121 pragma Unreferenced (Y, GCD);
|
|
|
122 begin
|
|
|
123 Euclid (P, Q, X, Y, GCD);
|
|
|
124 return X;
|
|
|
125 end Euclid;
|
|
|
126
|
|
|
127 -----------
|
|
|
128 -- Image --
|
|
|
129 -----------
|
|
|
130
|
|
|
131 function Image (Of_State : State) return String is
|
|
|
132 begin
|
|
|
133 return Int'Image (Of_State.X1) & ',' & Int'Image (Of_State.X2)
|
|
|
134 & ',' &
|
|
|
135 Int'Image (Of_State.P) & ',' & Int'Image (Of_State.Q);
|
|
|
136 end Image;
|
|
|
137
|
|
|
138 ------------
|
|
|
139 -- Random --
|
|
|
140 ------------
|
|
|
141
|
|
|
142 function Random (Gen : Generator) return Uniformly_Distributed is
|
|
|
143 Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access;
|
|
|
144
|
|
|
145 begin
|
|
|
146 Genp.X1 := Square_Mod_N (Genp.X1, Genp.P);
|
|
|
147 Genp.X2 := Square_Mod_N (Genp.X2, Genp.Q);
|
|
|
148 return
|
|
|
149 Float ((Flt (((Genp.X2 - Genp.X1) * Genp.X)
|
|
|
150 mod Genp.Q) * Flt (Genp.P)
|
|
|
151 + Flt (Genp.X1)) * Genp.Scl);
|
|
|
152 end Random;
|
|
|
153
|
|
|
154 -----------
|
|
|
155 -- Reset --
|
|
|
156 -----------
|
|
|
157
|
|
|
158 -- Version that works from given initiator value
|
|
|
159
|
|
|
160 procedure Reset (Gen : Generator; Initiator : Integer) is
|
|
|
161 Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access;
|
|
|
162 X1, X2 : Int;
|
|
|
163
|
|
|
164 begin
|
|
|
165 X1 := 2 + Int (Initiator) mod (K1 - 3);
|
|
|
166 X2 := 2 + Int (Initiator) mod (K2 - 3);
|
|
|
167
|
|
|
168 -- Eliminate effects of small initiators
|
|
|
169
|
|
|
170 for J in 1 .. 5 loop
|
|
|
171 X1 := Square_Mod_N (X1, K1);
|
|
|
172 X2 := Square_Mod_N (X2, K2);
|
|
|
173 end loop;
|
|
|
174
|
|
|
175 Genp.all :=
|
|
|
176 (X1 => X1,
|
|
|
177 X2 => X2,
|
|
|
178 P => K1,
|
|
|
179 Q => K2,
|
|
|
180 X => 1,
|
|
|
181 Scl => Scal);
|
|
|
182 end Reset;
|
|
|
183
|
|
|
184 -- Version that works from specific saved state
|
|
|
185
|
|
|
186 procedure Reset (Gen : Generator; From_State : State) is
|
|
|
187 Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access;
|
|
|
188
|
|
|
189 begin
|
|
|
190 Genp.all := From_State;
|
|
|
191 end Reset;
|
|
|
192
|
|
|
193 -- Version that works from calendar
|
|
|
194
|
|
|
195 procedure Reset (Gen : Generator) is
|
|
|
196 Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access;
|
|
|
197 Now : constant Calendar.Time := Calendar.Clock;
|
|
|
198 X1, X2 : Int;
|
|
|
199
|
|
|
200 begin
|
|
|
201 X1 := Int (Calendar.Year (Now)) * 12 * 31 +
|
|
|
202 Int (Calendar.Month (Now)) * 31 +
|
|
|
203 Int (Calendar.Day (Now));
|
|
|
204
|
|
|
205 X2 := Int (Calendar.Seconds (Now) * Duration (1000.0));
|
|
|
206
|
|
|
207 X1 := 2 + X1 mod (K1 - 3);
|
|
|
208 X2 := 2 + X2 mod (K2 - 3);
|
|
|
209
|
|
|
210 -- Eliminate visible effects of same day starts
|
|
|
211
|
|
|
212 for J in 1 .. 5 loop
|
|
|
213 X1 := Square_Mod_N (X1, K1);
|
|
|
214 X2 := Square_Mod_N (X2, K2);
|
|
|
215 end loop;
|
|
|
216
|
|
|
217 Genp.all :=
|
|
|
218 (X1 => X1,
|
|
|
219 X2 => X2,
|
|
|
220 P => K1,
|
|
|
221 Q => K2,
|
|
|
222 X => 1,
|
|
|
223 Scl => Scal);
|
|
|
224
|
|
|
225 end Reset;
|
|
|
226
|
|
|
227 ----------
|
|
|
228 -- Save --
|
|
|
229 ----------
|
|
|
230
|
|
|
231 procedure Save (Gen : Generator; To_State : out State) is
|
|
|
232 begin
|
|
|
233 To_State := Gen.Gen_State;
|
|
|
234 end Save;
|
|
|
235
|
|
|
236 ------------------
|
|
|
237 -- Square_Mod_N --
|
|
|
238 ------------------
|
|
|
239
|
|
|
240 function Square_Mod_N (X, N : Int) return Int is
|
|
|
241 Temp : constant Flt := Flt (X) * Flt (X);
|
|
|
242 Div : Int;
|
|
|
243
|
|
|
244 begin
|
|
|
245 Div := Int (Temp / Flt (N));
|
|
|
246 Div := Int (Temp - Flt (Div) * Flt (N));
|
|
|
247
|
|
|
248 if Div < 0 then
|
|
|
249 return Div + N;
|
|
|
250 else
|
|
|
251 return Div;
|
|
|
252 end if;
|
|
|
253 end Square_Mod_N;
|
|
|
254
|
|
|
255 -----------
|
|
|
256 -- Value --
|
|
|
257 -----------
|
|
|
258
|
|
|
259 function Value (Coded_State : String) return State is
|
|
|
260 Last : constant Natural := Coded_State'Last;
|
|
|
261 Start : Positive := Coded_State'First;
|
|
|
262 Stop : Positive := Coded_State'First;
|
|
|
263 Outs : State;
|
|
|
264
|
|
|
265 begin
|
|
|
266 while Stop <= Last and then Coded_State (Stop) /= ',' loop
|
|
|
267 Stop := Stop + 1;
|
|
|
268 end loop;
|
|
|
269
|
|
|
270 if Stop > Last then
|
|
|
271 raise Constraint_Error;
|
|
|
272 end if;
|
|
|
273
|
|
|
274 Outs.X1 := Int'Value (Coded_State (Start .. Stop - 1));
|
|
|
275 Start := Stop + 1;
|
|
|
276
|
|
|
277 loop
|
|
|
278 Stop := Stop + 1;
|
|
|
279 exit when Stop > Last or else Coded_State (Stop) = ',';
|
|
|
280 end loop;
|
|
|
281
|
|
|
282 if Stop > Last then
|
|
|
283 raise Constraint_Error;
|
|
|
284 end if;
|
|
|
285
|
|
|
286 Outs.X2 := Int'Value (Coded_State (Start .. Stop - 1));
|
|
|
287 Start := Stop + 1;
|
|
|
288
|
|
|
289 loop
|
|
|
290 Stop := Stop + 1;
|
|
|
291 exit when Stop > Last or else Coded_State (Stop) = ',';
|
|
|
292 end loop;
|
|
|
293
|
|
|
294 if Stop > Last then
|
|
|
295 raise Constraint_Error;
|
|
|
296 end if;
|
|
|
297
|
|
|
298 Outs.P := Int'Value (Coded_State (Start .. Stop - 1));
|
|
|
299 Outs.Q := Int'Value (Coded_State (Stop + 1 .. Last));
|
|
|
300 Outs.X := Euclid (Outs.P, Outs.Q);
|
|
|
301 Outs.Scl := 1.0 / (Flt (Outs.P) * Flt (Outs.Q));
|
|
|
302
|
|
|
303 -- Now do *some* sanity checks
|
|
|
304
|
|
|
305 if Outs.Q < 31 or else Outs.P < 31
|
|
|
306 or else Outs.X1 not in 2 .. Outs.P - 1
|
|
|
307 or else Outs.X2 not in 2 .. Outs.Q - 1
|
|
|
308 then
|
|
|
309 raise Constraint_Error;
|
|
|
310 end if;
|
|
|
311
|
|
|
312 return Outs;
|
|
|
313 end Value;
|
|
|
314 end GNAT.MBBS_Float_Random;
|
