]>
Commit | Line | Data |
---|---|---|
41195c94 PH |
1 | ------------------------------------------------------------------------------ |
2 | -- -- | |
3 | -- GNAT RUN-TIME COMPONENTS -- | |
4 | -- -- | |
1c612f29 | 5 | -- G N A T . M B B S _ D I S C R E T E _ R A N D O M -- |
41195c94 PH |
6 | -- -- |
7 | -- B o d y -- | |
8 | -- -- | |
4b490c1e | 9 | -- Copyright (C) 1992-2020, Free Software Foundation, Inc. -- |
41195c94 PH |
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 | with Interfaces; use Interfaces; | |
35 | ||
36 | package body GNAT.MBBS_Discrete_Random is | |
37 | ||
41195c94 PH |
38 | package Calendar renames Ada.Calendar; |
39 | ||
41195c94 PH |
40 | Fits_In_32_Bits : constant Boolean := |
41 | Rst'Size < 31 | |
42 | or else (Rst'Size = 31 | |
43 | and then Rst'Pos (Rst'First) < 0); | |
44 | -- This is set True if we do not need more than 32 bits in the result. If | |
45 | -- we need 64-bits, we will only use the meaningful 48 bits of any 64-bit | |
46 | -- number generated, since if more than 48 bits are required, we split the | |
47 | -- computation into two separate parts, since the algorithm does not behave | |
48 | -- above 48 bits. | |
49 | ||
50 | -- The way this expression works is that obviously if the size is 31 bits, | |
51 | -- it fits in 32 bits. In the 32-bit case, it fits in 32-bit signed if the | |
52 | -- range has negative values. It is too conservative in the case that the | |
53 | -- programmer has set a size greater than the default, e.g. a size of 33 | |
54 | -- for an integer type with a range of 1..10, but an over-conservative | |
55 | -- result is OK. The important thing is that the value is only True if | |
56 | -- we know the result will fit in 32-bits signed. If the value is False | |
57 | -- when it could be True, the behavior will be correct, just a bit less | |
58 | -- efficient than it could have been in some unusual cases. | |
59 | -- | |
60 | -- One might assume that we could get a more accurate result by testing | |
61 | -- the lower and upper bounds of the type Rst against the bounds of 32-bit | |
62 | -- Integer. However, there is no easy way to do that. Why? Because in the | |
308e6f3a | 63 | -- relatively rare case where this expression has to be evaluated at run |
41195c94 PH |
64 | -- time rather than compile time (when the bounds are dynamic), we need a |
65 | -- type to use for the computation. But the possible range of upper bound | |
66 | -- values for Rst (remembering the possibility of 64-bit modular types) is | |
67 | -- from -2**63 to 2**64-1, and no run-time type has a big enough range. | |
68 | ||
69 | ----------------------- | |
70 | -- Local Subprograms -- | |
71 | ----------------------- | |
72 | ||
73 | function Square_Mod_N (X, N : Int) return Int; | |
74 | pragma Inline (Square_Mod_N); | |
75 | -- Computes X**2 mod N avoiding intermediate overflow | |
76 | ||
77 | ----------- | |
78 | -- Image -- | |
79 | ----------- | |
80 | ||
81 | function Image (Of_State : State) return String is | |
82 | begin | |
83 | return Int'Image (Of_State.X1) & | |
84 | ',' & | |
85 | Int'Image (Of_State.X2) & | |
86 | ',' & | |
87 | Int'Image (Of_State.Q); | |
88 | end Image; | |
89 | ||
90 | ------------ | |
91 | -- Random -- | |
92 | ------------ | |
93 | ||
94 | function Random (Gen : Generator) return Rst is | |
9bebf0e9 | 95 | S : State renames Gen.Writable.Self.Gen_State; |
41195c94 PH |
96 | Temp : Int; |
97 | TF : Flt; | |
98 | ||
99 | begin | |
100 | -- Check for flat range here, since we are typically run with checks | |
101 | -- off, note that in practice, this condition will usually be static | |
102 | -- so we will not actually generate any code for the normal case. | |
103 | ||
104 | if Rst'Last < Rst'First then | |
105 | raise Constraint_Error; | |
106 | end if; | |
107 | ||
108 | -- Continue with computation if non-flat range | |
109 | ||
9bebf0e9 AC |
110 | S.X1 := Square_Mod_N (S.X1, S.P); |
111 | S.X2 := Square_Mod_N (S.X2, S.Q); | |
112 | Temp := S.X2 - S.X1; | |
41195c94 | 113 | |
a90bd866 | 114 | -- Following duplication is not an error, it is a loop unwinding |
41195c94 PH |
115 | |
116 | if Temp < 0 then | |
9bebf0e9 | 117 | Temp := Temp + S.Q; |
41195c94 PH |
118 | end if; |
119 | ||
120 | if Temp < 0 then | |
9bebf0e9 | 121 | Temp := Temp + S.Q; |
41195c94 PH |
122 | end if; |
123 | ||
9bebf0e9 | 124 | TF := Offs + (Flt (Temp) * Flt (S.P) + Flt (S.X1)) * S.Scl; |
41195c94 PH |
125 | |
126 | -- Pathological, but there do exist cases where the rounding implicit | |
127 | -- in calculating the scale factor will cause rounding to 'Last + 1. | |
128 | -- In those cases, returning 'First results in the least bias. | |
129 | ||
130 | if TF >= Flt (Rst'Pos (Rst'Last)) + 0.5 then | |
131 | return Rst'First; | |
132 | ||
133 | elsif not Fits_In_32_Bits then | |
134 | return Rst'Val (Interfaces.Integer_64 (TF)); | |
135 | ||
136 | else | |
137 | return Rst'Val (Int (TF)); | |
138 | end if; | |
139 | end Random; | |
140 | ||
141 | ----------- | |
142 | -- Reset -- | |
143 | ----------- | |
144 | ||
145 | procedure Reset (Gen : Generator; Initiator : Integer) is | |
9bebf0e9 | 146 | S : State renames Gen.Writable.Self.Gen_State; |
41195c94 PH |
147 | X1, X2 : Int; |
148 | ||
149 | begin | |
150 | X1 := 2 + Int (Initiator) mod (K1 - 3); | |
151 | X2 := 2 + Int (Initiator) mod (K2 - 3); | |
152 | ||
153 | for J in 1 .. 5 loop | |
154 | X1 := Square_Mod_N (X1, K1); | |
155 | X2 := Square_Mod_N (X2, K2); | |
156 | end loop; | |
157 | ||
158 | -- Eliminate effects of small Initiators | |
159 | ||
9bebf0e9 | 160 | S := |
41195c94 PH |
161 | (X1 => X1, |
162 | X2 => X2, | |
163 | P => K1, | |
164 | Q => K2, | |
165 | FP => K1F, | |
166 | Scl => Scal); | |
167 | end Reset; | |
168 | ||
169 | ----------- | |
170 | -- Reset -- | |
171 | ----------- | |
172 | ||
173 | procedure Reset (Gen : Generator) is | |
9bebf0e9 | 174 | S : State renames Gen.Writable.Self.Gen_State; |
41195c94 PH |
175 | Now : constant Calendar.Time := Calendar.Clock; |
176 | X1 : Int; | |
177 | X2 : Int; | |
178 | ||
179 | begin | |
180 | X1 := Int (Calendar.Year (Now)) * 12 * 31 + | |
b35e5dcb | 181 | Int (Calendar.Month (Now) * 31) + |
41195c94 PH |
182 | Int (Calendar.Day (Now)); |
183 | ||
184 | X2 := Int (Calendar.Seconds (Now) * Duration (1000.0)); | |
185 | ||
186 | X1 := 2 + X1 mod (K1 - 3); | |
187 | X2 := 2 + X2 mod (K2 - 3); | |
188 | ||
189 | -- Eliminate visible effects of same day starts | |
190 | ||
191 | for J in 1 .. 5 loop | |
192 | X1 := Square_Mod_N (X1, K1); | |
193 | X2 := Square_Mod_N (X2, K2); | |
194 | end loop; | |
195 | ||
9bebf0e9 | 196 | S := |
41195c94 PH |
197 | (X1 => X1, |
198 | X2 => X2, | |
199 | P => K1, | |
200 | Q => K2, | |
201 | FP => K1F, | |
202 | Scl => Scal); | |
203 | ||
204 | end Reset; | |
205 | ||
206 | ----------- | |
207 | -- Reset -- | |
208 | ----------- | |
209 | ||
210 | procedure Reset (Gen : Generator; From_State : State) is | |
41195c94 | 211 | begin |
9bebf0e9 | 212 | Gen.Writable.Self.Gen_State := From_State; |
41195c94 PH |
213 | end Reset; |
214 | ||
215 | ---------- | |
216 | -- Save -- | |
217 | ---------- | |
218 | ||
219 | procedure Save (Gen : Generator; To_State : out State) is | |
220 | begin | |
221 | To_State := Gen.Gen_State; | |
222 | end Save; | |
223 | ||
224 | ------------------ | |
225 | -- Square_Mod_N -- | |
226 | ------------------ | |
227 | ||
228 | function Square_Mod_N (X, N : Int) return Int is | |
229 | begin | |
230 | return Int ((Integer_64 (X) ** 2) mod (Integer_64 (N))); | |
231 | end Square_Mod_N; | |
232 | ||
233 | ----------- | |
234 | -- Value -- | |
235 | ----------- | |
236 | ||
237 | function Value (Coded_State : String) return State is | |
238 | Last : constant Natural := Coded_State'Last; | |
239 | Start : Positive := Coded_State'First; | |
240 | Stop : Positive := Coded_State'First; | |
241 | Outs : State; | |
242 | ||
243 | begin | |
244 | while Stop <= Last and then Coded_State (Stop) /= ',' loop | |
245 | Stop := Stop + 1; | |
246 | end loop; | |
247 | ||
248 | if Stop > Last then | |
249 | raise Constraint_Error; | |
250 | end if; | |
251 | ||
252 | Outs.X1 := Int'Value (Coded_State (Start .. Stop - 1)); | |
253 | Start := Stop + 1; | |
254 | ||
255 | loop | |
256 | Stop := Stop + 1; | |
257 | exit when Stop > Last or else Coded_State (Stop) = ','; | |
258 | end loop; | |
259 | ||
260 | if Stop > Last then | |
261 | raise Constraint_Error; | |
262 | end if; | |
263 | ||
264 | Outs.X2 := Int'Value (Coded_State (Start .. Stop - 1)); | |
265 | Outs.Q := Int'Value (Coded_State (Stop + 1 .. Last)); | |
266 | Outs.P := Outs.Q * 2 + 1; | |
267 | Outs.FP := Flt (Outs.P); | |
268 | Outs.Scl := (RstL - RstF + 1.0) / (Flt (Outs.P) * Flt (Outs.Q)); | |
269 | ||
270 | -- Now do *some* sanity checks | |
271 | ||
272 | if Outs.Q < 31 | |
273 | or else Outs.X1 not in 2 .. Outs.P - 1 | |
274 | or else Outs.X2 not in 2 .. Outs.Q - 1 | |
275 | then | |
276 | raise Constraint_Error; | |
277 | end if; | |
278 | ||
279 | return Outs; | |
280 | end Value; | |
281 | ||
282 | end GNAT.MBBS_Discrete_Random; |