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-nudira.adb < prev next >

Wrap

Text File | 1995-10-17 | 8.2 KB | 256 lines

------------------------------------------------------------------------------ -- -- -- GNAT RUNTIME COMPONENTS -- -- -- -- A D A . N U M E R I C S . D I S C R E T E _ R A N D O M -- -- -- -- B o d y -- -- -- -- $Revision: 1.9 $ -- -- -- -- 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.Discrete_Random is ----------------------- -- Local Subprograms -- ----------------------- function Square_Mod_N (X, N : Int_32) return Int_32; pragma Inline (Square_Mod_N); procedure Test_Math (N : Int_32); ------------ -- Create -- ------------ function Create return Pointer is begin return new State' (X1 => 2999 ** 2, X2 => 1439 ** 2, P => 94_833_359, Q => 47_416_679, FP => 94_833_359.0, Offset => Floating (Result_Subtype'Pos (Result_Subtype'First)) + 0.5, Scale => (Floating (Result_Subtype'Pos (Result_Subtype'Last)) - Floating (Result_Subtype'Pos (Result_Subtype'First)) + 1.0) / (94_833_359.0 * 47_416_679.0)); end Create; ----------- -- 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); end Image; ------------ -- Random -- ------------ function Random (Gen : Generator) return Result_Subtype is 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); Temp := Gen.Point.X2 - Gen.Point.X1; if Temp < 0 then Temp := Temp + Gen.Point.Q; end if; if Temp < 0 then Temp := Temp + Gen.Point.Q; end if; -- duplication is not an error, it is a loop unwinding. Temp := Int_32 (Gen.Point.Offset + (Floating (Temp) * Floating (Gen.Point.P) + Floating (Gen.Point.X1)) * Gen.Point.Scale); if Temp > Int_32 (Result_Subtype'Pos (Result_Subtype'Last)) then return Result_Subtype'First; else return Result_Subtype'Val (Temp); end if; -- Pathological, but there do exist cases where the rounding implicit -- in calculating the scale factor will cause rounding to 'Last + 1. -- In those cases, returning 'First results in the least bias. end Random; ----------- -- Reset -- ----------- procedure Reset (Gen : Generator; Initiator : Integer) is X1, X2, P, Q : Int_32; begin P := 94_833_359; Q := 47_416_679; 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, FP => Floating (P), Offset => Floating (Result_Subtype'Pos (Result_Subtype'First)) + 0.5, Scale => (Floating (Result_Subtype'Pos (Result_Subtype'Last)) - Floating (Result_Subtype'Pos (Result_Subtype'First)) + 1.0) / (Floating (P) * Floating (Q))); end Reset; ----------- -- Reset -- ----------- procedure Reset (Gen : Generator) is use Ada.Calendar; Now : Time := Clock; X1, X2, P, Q : Int_32; begin P := 94_833_359; Q := 47_416_679; X1 := Int_32 (Year (Now)) * 12 * 31 + Int_32 (Month (Now) * 31) + Int_32 (Day (Now)); X2 := Int_32 (Seconds (Now) * Duration (1000.0)); 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, FP => Floating (P), Offset => Floating (Result_Subtype'Pos (Result_Subtype'First)) + 0.5, Scale => (Floating (Result_Subtype'Pos (Result_Subtype'Last)) - Floating (Result_Subtype'Pos (Result_Subtype'First)) + 1.0) / (Floating (P) * Floating (Q))); -- Why save the actual assignments to the end? To insure to the -- greatest extent possible that an exception won't leave the generator -- inconsistent. end Reset; ----------- -- Reset -- ----------- procedure Reset (Gen : Generator; From_State : State) is begin Gen.Point.all := From_State; end Reset; ---------- -- Save -- ---------- procedure Save (Gen : Generator; To_State : out State) is begin To_State := Gen.Point.all; end Save; ------------------ -- Square_Mod_N -- ------------------ function Square_Mod_N (X, N : Int_32) return Int_32 is Temp : Floating := Floating (X) * Floating (X); Div : Int_32 := Int_32 (Temp / Floating (N)); begin Div := Int_32 (Temp - Floating (Div) * Floating (N)); if Div < 0 then return Div + N; else return Div; end if; end Square_Mod_N; --------------- -- Test_Math -- --------------- procedure Test_Math (N : 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)); Out_State.Q := Int_32'Value (Coded_State (Stop + 1 .. Coded_State'Last)); Out_State.P := Out_State.Q * 2 + 1; Out_State.FP := Floating (Out_State.P); Out_State.Scale := (Floating (Result_Subtype'Pos (Result_Subtype'Last)) - Floating (Result_Subtype'Pos (Result_Subtype'First)) + 1.0) / (Floating (Out_State.P) * Floating (Out_State.Q)); -- now do SOME sanity checks. if Out_State.Q < 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); return Out_State; end Value; begin Test_Math (94_833_359); end Ada.Numerics.Discrete_Random;