Geek Gadgets 1

home *** CD-ROM | disk | FTP | other *** search

/ Geek Gadgets 1 / ADE-1.bin / ade-bin / lib / gcc-lib / m68k-amigaos / 2.7.2.1 / adainclude / a-nuflra.adb < prev next >

Wrap

Text File | 1995-10-17 | 9.0 KB | 299 lines

------------------------------------------------------------------------------ -- -- -- GNAT RUNTIME COMPONENTS -- -- -- -- A D A . N U M E R I C S . F L O A T _ R A N D O M -- -- -- -- B o d y -- -- -- -- $Revision: 1.11 $ -- -- -- -- Copyright (c) 1992,1993,1994 NYU, All Rights Reserved -- -- -- -- The GNAT library is free software; you can redistribute it and/or modify -- -- it under terms of the GNU Library General Public License as published by -- -- the Free Software Foundation; either version 2, or (at your option) any -- -- later version. The GNAT library is distributed in the hope that it will -- -- be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -- -- of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -- -- Library General Public License for more details. You should have -- -- received a copy of the GNU Library General Public License along with -- -- the GNAT library; see the file COPYING.LIB. If not, write to the Free -- -- Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -- -- -- ------------------------------------------------------------------------------ with Ada.Calendar; package body Ada.Numerics.Float_Random is ----------------------- -- Local Subprograms -- ----------------------- procedure Euclid (P, Q : in Int_32; X, Y : out Int_32; GCD : out Int_32); function Euclid (P, Q : Int_32) return Int_32; function Square_Mod_N (X, N : Int_32) return Int_32; pragma Inline (Square_Mod_N); procedure Test_Math (N : in Int_32); ------------ -- Create -- ------------ function Create return Pointer is begin return new State' (Initial_State); end Create; ------------ -- Euclid -- ------------ procedure Euclid (P, Q : in Int_32; X, Y : out Int_32; GCD : out Int_32) is XT : Int_32 := 1; YT : Int_32 := 0; procedure Recur (P, Q : in Int_32; -- a (i-1), a (i) X, Y : in Int_32; -- x (i), y (i) XP, YP : in out Int_32; -- x (i-1), y (i-1) GCD : out Int_32); procedure Recur (P, Q : in Int_32; X, Y : in Int_32; XP, YP : in out Int_32; GCD : out Int_32) is Quo : Int_32 := P / Q; -- q <-- |_ a (i-1) / a (i) _| XT : Int_32 := X; -- x (i) YT : Int_32 := Y; -- y (i) begin if P rem Q = 0 then -- while does not divide GCD := Q; XP := X; YP := Y; else Recur (Q, P - Q * Quo, XP - Quo * X, YP - Quo * Y, XT, YT, Quo); -- a (i) <== a (i) -- a (i+1) <-- a (i-1) - q*a (i) -- x (i+1) <-- x (i-1) - q*x (i) -- y (i+1) <-- y (i-1) - q*y (i) -- x (i) <== x (i) -- y (i) <== y (i) XP := XT; YP := YT; GCD := Quo; end if; end Recur; begin Recur (P, Q, 0, 1, XT, YT, GCD); X := XT; Y := YT; end Euclid; ------------ -- Euclid -- ------------ function Euclid (P, Q : Int_32) return Int_32 is X, Y, GCD : Int_32; begin Euclid (P, Q, X, Y, GCD); return X; end Euclid; ----------- -- Image -- ----------- function Image (Of_State : State) return String is begin return Int_32'Image (Of_State.X1) & ',' & Int_32'Image (Of_State.X2) & ',' & Int_32'Image (Of_State.P) & ',' & Int_32'Image (Of_State.Q) & ',' & Int_32'Image (Of_State.X); end Image; ------------ -- Random -- ------------ function Random (Gen : Generator) return Uniformly_Distributed is X1, X2, Temp : Int_32; begin Gen.Point.X1 := Square_Mod_N (Gen.Point.X1, Gen.Point.P); Gen.Point.X2 := Square_Mod_N (Gen.Point.X2, Gen.Point.Q); return Float ((Longer_Float (((Gen.Point.X2 - Gen.Point.X1) * Gen.Point.X) mod Gen.Point.Q) * Longer_Float (Gen.Point.P) + Longer_Float (Gen.Point.X1)) * Gen.Point.Scale); end Random; ----------- -- Reset -- ----------- procedure Reset (Gen : in Generator; Initiator : in Integer) is X1, X2, P, Q : Int_32; begin P := Initial_State.P; Q := Initial_State.Q; X1 := 2 + Int_32 (Initiator) rem (P - 3); X2 := 2 + Int_32 (Initiator) rem (Q - 3); for I in 1 .. 5 loop X1 := Square_Mod_N (X1, P); X2 := Square_Mod_N (X2, Q); end loop; -- eliminate effects of small Initiators. Gen.Point.all := (X1 => X1, X2 => X2, P => P, Q => Q, X => Euclid (P, Q), Scale => 1.0 / (Longer_Float (P) * Longer_Float (Q))); end Reset; ----------- -- Reset -- ----------- procedure Reset (Gen : Generator; From_State : State) is begin Gen.Point.all := From_State; end Reset; ----------- -- Reset -- ----------- procedure Reset (Gen : Generator) is use Ada.Calendar; Now : Time := Clock; X1, X2, P, Q : Int_32; begin P := Initial_State.P; Q := Initial_State.Q; X1 := Int_32 (Year (Now)) * 12 * 31 + Int_32 (Month (Now)) * 31 + Int_32 (Day (Now)); X2 := Int_32 (Seconds (Now) * Duration (1000.0)); -- X2 := Int_32 (1000 * Seconds (Now)); raises Constraint_Error X1 := 2 + X1 rem (P - 3); X2 := 2 + X2 rem (Q - 3); for I in 1 .. 5 loop X1 := Square_Mod_N (X1, P); X2 := Square_Mod_N (X2, Q); end loop; -- less justification here, but eliminate visible effects of same day -- starts. Gen.Point.all := (X1 => X1, X2 => X2, P => P, Q => Q, X => Euclid (P, Q), Scale => 1.0 / (Longer_Float (P) * Longer_Float (Q))); -- Why save the actual assignments to the end? To insure to the -- greatest extent possible that an exception won't leave the generator -- inconsistant. end Reset; ---------- -- Save -- ---------- procedure Save (Gen : in Generator; To_State : out State) is begin To_State := Gen.Point.all; end Save; ------------------ -- Square_Mod_N -- ------------------ -- don't mess with working code, but I want a function that returns X. -- Note rem is used below instead of mod where both arguments are known -- positive. This is faster on some hardware. function Square_Mod_N (X, N : Int_32) return Int_32 is subtype LF is Longer_Float; Temp : LF := LF (X) * LF (X); Div : Int_32 := Int_32 (Temp / LF (N)); begin Div := Int_32 (Temp - LF (Div) * LF (N)); if Div < 0 then return Div + N; else return Div; end if; end Square_Mod_N; --------------- -- Test_Math -- --------------- procedure Test_Math (N : in Int_32) is begin if Square_Mod_N (N - 1, N) /= 1 then raise Program_Error; end if; end Test_Math; ----------- -- Value -- ----------- function Value (Coded_State : String) return State is Start, Stop : Positive := Coded_State'First; Out_State : State; begin while Coded_State (Stop) /= ',' loop Stop := Stop + 1; end loop; Out_State.X1 := Int_32'Value (Coded_State (Start .. Stop - 1)); Start := Stop + 1; loop Stop := Stop + 1; exit when Coded_State (Stop) = ','; end loop; Out_State.X2 := Int_32'Value (Coded_State (Start .. Stop - 1)); Start := Stop + 1; loop Stop := Stop + 1; exit when Coded_State (Stop) = ','; end loop; Out_State.P := Int_32'Value (Coded_State (Start .. Stop - 1)); Out_State.Q := Int_32'Value (Coded_State (Stop + 1 .. Coded_State'Last)); Out_State.X := Euclid (Out_State.P, Out_State.Q); Out_State.Scale := 1.0 / (Longer_Float (Out_State.P) * Longer_Float (Out_State.Q)); -- now do SOME sanity checks. if Out_State.Q < 31 or else Out_State.P < 31 or else Out_State.X1 not in 2 .. Out_State.P - 1 or else Out_State.X2 not in 2 .. Out_State.Q - 1 then raise Constraint_Error; end if; Test_Math (Out_State.P); Test_Math (Out_State.Q); return Out_State; end Value; begin Test_Math (94_833_359); end Ada.Numerics.Float_Random;