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-- -- 9-- Copyright (C) 1992-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 32with Ada.Calendar; 33 34package 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; 314end GNAT.MBBS_Float_Random; 315