From 41195c94583a5f2da9178234bbcdacee27db69ff Mon Sep 17 00:00:00 2001 From: Paul Hilfinger Date: Tue, 22 Jun 2010 15:24:10 +0000 Subject: [PATCH] 2010-06-22 Paul Hilfinger * a-nudira.adb, a-nudira.ads, a-nuflra.adb, a-nuflra.ads, gnat_rm.texi, impunit.adb, Makefile.rtl, s-rannum.adb (Random_Float_Template, Random): New method of creating uniform floating-point variables that allow the creation of all machine values in [0 .. 1). * g-mbdira.adb, g-mbflra.adb, g-mbdira.ads, g-mbflra.ads: New file. From-SVN: r161191 --- gcc/ada/ChangeLog | 10 ++ gcc/ada/Makefile.rtl | 4 +- gcc/ada/a-nudira.adb | 260 +++++------------------------------ gcc/ada/a-nudira.ads | 57 ++------ gcc/ada/a-nuflra.adb | 252 +++++----------------------------- gcc/ada/a-nuflra.ads | 41 ++---- gcc/ada/g-mbdira.adb | 299 +++++++++++++++++++++++++++++++++++++++++ gcc/ada/g-mbdira.ads | 118 ++++++++++++++++ gcc/ada/g-mbflra.adb | 313 +++++++++++++++++++++++++++++++++++++++++++ gcc/ada/g-mbflra.ads | 103 ++++++++++++++ gcc/ada/gnat_rm.texi | 18 +-- gcc/ada/impunit.adb | 2 + gcc/ada/s-rannum.adb | 111 ++++++++++++--- 13 files changed, 1029 insertions(+), 559 deletions(-) create mode 100644 gcc/ada/g-mbdira.adb create mode 100644 gcc/ada/g-mbdira.ads create mode 100644 gcc/ada/g-mbflra.adb create mode 100644 gcc/ada/g-mbflra.ads diff --git a/gcc/ada/ChangeLog b/gcc/ada/ChangeLog index 3bb621317f7..700fa324ac2 100644 --- a/gcc/ada/ChangeLog +++ b/gcc/ada/ChangeLog @@ -1,3 +1,13 @@ +2010-06-22 Paul Hilfinger + + * a-nudira.adb, a-nudira.ads, a-nuflra.adb, a-nuflra.ads, + gnat_rm.texi, impunit.adb, Makefile.rtl, s-rannum.adb + (Random_Float_Template, Random): New method of creating + uniform floating-point variables that allow the creation of all machine + values in [0 .. 1). + + * g-mbdira.adb, g-mbflra.adb, g-mbdira.ads, g-mbflra.ads: New file. + 2010-06-22 Gary Dismukes * sem_ch5.adb (Analyze_Assignment): Revise test for illegal assignment diff --git a/gcc/ada/Makefile.rtl b/gcc/ada/Makefile.rtl index f101a52e025..c130ad31a83 100644 --- a/gcc/ada/Makefile.rtl +++ b/gcc/ada/Makefile.rtl @@ -1,5 +1,5 @@ # Makefile.rtl for GNU Ada Compiler (GNAT). -# Copyright (C) 2003-2008, Free Software Foundation, Inc. +# Copyright (C) 2003-2010, Free Software Foundation, Inc. #This file is part of GCC. @@ -359,6 +359,8 @@ GNATRTL_NONTASKING_OBJS= \ g-io$(objext) \ g-io_aux$(objext) \ g-locfil$(objext) \ + g-mbdira$(objext) \ + g-mbflra$(objext) \ g-md5$(objext) \ g-memdum$(objext) \ g-moreex$(objext) \ diff --git a/gcc/ada/a-nudira.adb b/gcc/ada/a-nudira.adb index 87abcd8f100..e17945c07a2 100644 --- a/gcc/ada/a-nudira.adb +++ b/gcc/ada/a-nudira.adb @@ -6,7 +6,7 @@ -- -- -- B o d y -- -- -- --- Copyright (C) 1992-2009, Free Software Foundation, Inc. -- +-- Copyright (C) 1992-2010, Free Software Foundation, Inc. -- -- -- -- GNAT is free software; you can redistribute it and/or modify it under -- -- terms of the GNU General Public License as published by the Free Soft- -- @@ -29,9 +29,7 @@ -- -- ------------------------------------------------------------------------------ -with Ada.Calendar; - -with Interfaces; use Interfaces; +with System.Random_Numbers; use System.Random_Numbers; package body Ada.Numerics.Discrete_Random is @@ -49,249 +47,55 @@ package body Ada.Numerics.Discrete_Random is -- get a pointer to the state in the passed Generator. This works because -- Generator is a limited type and will thus always be passed by reference. - type Pointer is access all State; - - Fits_In_32_Bits : constant Boolean := - Rst'Size < 31 - or else (Rst'Size = 31 - and then Rst'Pos (Rst'First) < 0); - -- This is set True if we do not need more than 32 bits in the result. If - -- we need 64-bits, we will only use the meaningful 48 bits of any 64-bit - -- number generated, since if more than 48 bits are required, we split the - -- computation into two separate parts, since the algorithm does not behave - -- above 48 bits. - - -- The way this expression works is that obviously if the size is 31 bits, - -- it fits in 32 bits. In the 32-bit case, it fits in 32-bit signed if the - -- range has negative values. It is too conservative in the case that the - -- programmer has set a size greater than the default, e.g. a size of 33 - -- for an integer type with a range of 1..10, but an over-conservative - -- result is OK. The important thing is that the value is only True if - -- we know the result will fit in 32-bits signed. If the value is False - -- when it could be True, the behavior will be correct, just a bit less - -- efficient than it could have been in some unusual cases. - -- - -- One might assume that we could get a more accurate result by testing - -- the lower and upper bounds of the type Rst against the bounds of 32-bit - -- Integer. However, there is no easy way to do that. Why? Because in the - -- relatively rare case where this expresion has to be evaluated at run - -- time rather than compile time (when the bounds are dynamic), we need a - -- type to use for the computation. But the possible range of upper bound - -- values for Rst (remembering the possibility of 64-bit modular types) is - -- from -2**63 to 2**64-1, and no run-time type has a big enough range. - - ----------------------- - -- Local Subprograms -- - ----------------------- + subtype Rep_Generator is System.Random_Numbers.Generator; + subtype Rep_State is System.Random_Numbers.State; - function Square_Mod_N (X, N : Int) return Int; - pragma Inline (Square_Mod_N); - -- Computes X**2 mod N avoiding intermediate overflow + function Rep_Random is + new Random_Discrete (Result_Subtype, Result_Subtype'First); - ----------- - -- Image -- - ----------- - - function Image (Of_State : State) return String is + function Random (Gen : Generator) return Result_Subtype is begin - return Int'Image (Of_State.X1) & - ',' & - Int'Image (Of_State.X2) & - ',' & - Int'Image (Of_State.Q); - end Image; - - ------------ - -- Random -- - ------------ - - function Random (Gen : Generator) return Rst is - Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; - Temp : Int; - TF : Flt; - - begin - -- Check for flat range here, since we are typically run with checks - -- off, note that in practice, this condition will usually be static - -- so we will not actually generate any code for the normal case. - - if Rst'Last < Rst'First then - raise Constraint_Error; - end if; - - -- Continue with computation if non-flat range - - Genp.X1 := Square_Mod_N (Genp.X1, Genp.P); - Genp.X2 := Square_Mod_N (Genp.X2, Genp.Q); - Temp := Genp.X2 - Genp.X1; - - -- Following duplication is not an error, it is a loop unwinding! - - if Temp < 0 then - Temp := Temp + Genp.Q; - end if; - - if Temp < 0 then - Temp := Temp + Genp.Q; - end if; - - TF := Offs + (Flt (Temp) * Flt (Genp.P) + Flt (Genp.X1)) * Genp.Scl; - - -- 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. - - if TF >= Flt (Rst'Pos (Rst'Last)) + 0.5 then - return Rst'First; - - elsif not Fits_In_32_Bits then - return Rst'Val (Interfaces.Integer_64 (TF)); - - else - return Rst'Val (Int (TF)); - end if; + return Rep_Random (Gen.Rep); end Random; - ----------- - -- Reset -- - ----------- - - procedure Reset (Gen : Generator; Initiator : Integer) is - Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; - X1, X2 : Int; - - begin - X1 := 2 + Int (Initiator) mod (K1 - 3); - X2 := 2 + Int (Initiator) mod (K2 - 3); - - for J in 1 .. 5 loop - X1 := Square_Mod_N (X1, K1); - X2 := Square_Mod_N (X2, K2); - end loop; - - -- Eliminate effects of small Initiators - - Genp.all := - (X1 => X1, - X2 => X2, - P => K1, - Q => K2, - FP => K1F, - Scl => Scal); - end Reset; - - ----------- - -- Reset -- - ----------- - - procedure Reset (Gen : Generator) is - Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; - Now : constant Calendar.Time := Calendar.Clock; - X1 : Int; - X2 : Int; - + procedure Reset (Gen : Generator; + Initiator : Integer) is + G : Rep_Generator renames Gen.Rep'Unrestricted_Access.all; begin - X1 := Int (Calendar.Year (Now)) * 12 * 31 + - Int (Calendar.Month (Now) * 31) + - Int (Calendar.Day (Now)); - - X2 := Int (Calendar.Seconds (Now) * Duration (1000.0)); - - X1 := 2 + X1 mod (K1 - 3); - X2 := 2 + X2 mod (K2 - 3); - - -- Eliminate visible effects of same day starts - - for J in 1 .. 5 loop - X1 := Square_Mod_N (X1, K1); - X2 := Square_Mod_N (X2, K2); - end loop; - - Genp.all := - (X1 => X1, - X2 => X2, - P => K1, - Q => K2, - FP => K1F, - Scl => Scal); - + Reset (G, Initiator); end Reset; - ----------- - -- Reset -- - ----------- - - procedure Reset (Gen : Generator; From_State : State) is - Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; + procedure Reset (Gen : Generator) is + G : Rep_Generator renames Gen.Rep'Unrestricted_Access.all; begin - Genp.all := From_State; + Reset (G); end Reset; - ---------- - -- Save -- - ---------- - - procedure Save (Gen : Generator; To_State : out State) is + procedure Save (Gen : Generator; + To_State : out State) is begin - To_State := Gen.Gen_State; + Save (Gen.Rep, State (To_State)); end Save; - ------------------ - -- Square_Mod_N -- - ------------------ - - function Square_Mod_N (X, N : Int) return Int is + procedure Reset (Gen : Generator; + From_State : State) is + G : Rep_Generator renames Gen.Rep'Unrestricted_Access.all; begin - return Int ((Integer_64 (X) ** 2) mod (Integer_64 (N))); - end Square_Mod_N; + Reset (G, From_State); + end Reset; - ----------- - -- Value -- - ----------- + function Image (Of_State : State) return String is + begin + return Image (Rep_State (Of_State)); + end Image; function Value (Coded_State : String) return State is - Last : constant Natural := Coded_State'Last; - Start : Positive := Coded_State'First; - Stop : Positive := Coded_State'First; - Outs : State; - + G : Generator; + S : Rep_State; begin - while Stop <= Last and then Coded_State (Stop) /= ',' loop - Stop := Stop + 1; - end loop; - - if Stop > Last then - raise Constraint_Error; - end if; - - Outs.X1 := Int'Value (Coded_State (Start .. Stop - 1)); - Start := Stop + 1; - - loop - Stop := Stop + 1; - exit when Stop > Last or else Coded_State (Stop) = ','; - end loop; - - if Stop > Last then - raise Constraint_Error; - end if; - - Outs.X2 := Int'Value (Coded_State (Start .. Stop - 1)); - Outs.Q := Int'Value (Coded_State (Stop + 1 .. Last)); - Outs.P := Outs.Q * 2 + 1; - Outs.FP := Flt (Outs.P); - Outs.Scl := (RstL - RstF + 1.0) / (Flt (Outs.P) * Flt (Outs.Q)); - - -- Now do *some* sanity checks - - if Outs.Q < 31 - or else Outs.X1 not in 2 .. Outs.P - 1 - or else Outs.X2 not in 2 .. Outs.Q - 1 - then - raise Constraint_Error; - end if; - - return Outs; + Reset (G.Rep, Coded_State); + System.Random_Numbers.Save (G.Rep, S); + return State (S); end Value; end Ada.Numerics.Discrete_Random; diff --git a/gcc/ada/a-nudira.ads b/gcc/ada/a-nudira.ads index 425aa6f9bc9..03ce48b38b4 100644 --- a/gcc/ada/a-nudira.ads +++ b/gcc/ada/a-nudira.ads @@ -6,7 +6,7 @@ -- -- -- S p e c -- -- -- --- Copyright (C) 1992-2009, Free Software Foundation, Inc. -- +-- Copyright (C) 1992-2010, Free Software Foundation, Inc. -- -- -- -- This specification is derived from the Ada Reference Manual for use with -- -- GNAT. The copyright notice above, and the license provisions that follow -- @@ -33,39 +33,24 @@ -- -- ------------------------------------------------------------------------------ --- Note: the implementation used in this package was contributed by Robert --- Eachus. It is based on the work of L. Blum, M. Blum, and M. Shub, SIAM --- Journal of Computing, Vol 15. No 2, May 1986. The particular choices for P --- and Q chosen here guarantee a period of 562,085,314,430,582 (about 2**49), --- and the generated sequence has excellent randomness properties. For further --- details, see the paper "Fast Generation of Trustworthy Random Numbers", by --- Robert Eachus, which describes both the algorithm and the efficient --- implementation approach used here. +-- Note: the implementation used in this package is a version of the +-- Mersenne Twister. See s-rannum.adb for details and references. -with Interfaces; +with System.Random_Numbers; generic type Result_Subtype is (<>); package Ada.Numerics.Discrete_Random is - -- The algorithm used here is reliable from a required statistical point of - -- view only up to 48 bits. We try to behave reasonably in the case of - -- larger types, but we can't guarantee the required properties. So - -- generate a warning for these (slightly) dubious cases. - - pragma Compile_Time_Warning - (Result_Subtype'Size > 48, - "statistical properties not guaranteed for size > 48"); - -- Basic facilities type Generator is limited private; function Random (Gen : Generator) return Result_Subtype; - procedure Reset (Gen : Generator); procedure Reset (Gen : Generator; Initiator : Integer); + procedure Reset (Gen : Generator); -- Advanced facilities @@ -74,41 +59,17 @@ package Ada.Numerics.Discrete_Random is procedure Save (Gen : Generator; To_State : out State); procedure Reset (Gen : Generator; From_State : State); - Max_Image_Width : constant := 80; + Max_Image_Width : constant := System.Random_Numbers.Max_Image_Width; function Image (Of_State : State) return String; function Value (Coded_State : String) return State; private - subtype Int is Interfaces.Integer_32; - subtype Rst is Result_Subtype; - - -- We prefer to use 14 digits for Flt, but some targets are more limited - - type Flt is digits Positive'Min (14, Long_Long_Float'Digits); - - RstF : constant Flt := Flt (Rst'Pos (Rst'First)); - RstL : constant Flt := Flt (Rst'Pos (Rst'Last)); - - Offs : constant Flt := RstF - 0.5; - - K1 : constant := 94_833_359; - K1F : constant := 94_833_359.0; - K2 : constant := 47_416_679; - K2F : constant := 47_416_679.0; - Scal : constant Flt := (RstL - RstF + 1.0) / (K1F * K2F); - - type State is record - X1 : Int := Int (2999 ** 2); - X2 : Int := Int (1439 ** 2); - P : Int := K1; - Q : Int := K2; - FP : Flt := K1F; - Scl : Flt := Scal; - end record; type Generator is limited record - Gen_State : State; + Rep : System.Random_Numbers.Generator; end record; + type State is new System.Random_Numbers.State; + end Ada.Numerics.Discrete_Random; diff --git a/gcc/ada/a-nuflra.adb b/gcc/ada/a-nuflra.adb index 7e6323b8e8d..e58ff9247c2 100644 --- a/gcc/ada/a-nuflra.adb +++ b/gcc/ada/a-nuflra.adb @@ -6,7 +6,7 @@ -- -- -- B o d y -- -- -- --- Copyright (C) 1992-2009, Free Software Foundation, Inc. -- +-- Copyright (C) 1992-2010, Free Software Foundation, Inc. -- -- -- -- GNAT is free software; you can redistribute it and/or modify it under -- -- terms of the GNU General Public License as published by the Free Soft- -- @@ -29,7 +29,9 @@ -- -- ------------------------------------------------------------------------------ -with Ada.Calendar; +with Interfaces; use Interfaces; + +with System.Random_Numbers; use System.Random_Numbers; package body Ada.Numerics.Float_Random is @@ -47,105 +49,16 @@ package body Ada.Numerics.Float_Random is -- get a pointer to the state in the passed Generator. This works because -- Generator is a limited type and will thus always be passed by reference. - type Pointer is access all State; - - ----------------------- - -- Local Subprograms -- - ----------------------- - - procedure Euclid (P, Q : Int; X, Y : out Int; GCD : out Int); - - function Euclid (P, Q : Int) return Int; - - function Square_Mod_N (X, N : Int) return Int; - - ------------ - -- Euclid -- - ------------ - - procedure Euclid (P, Q : Int; X, Y : out Int; GCD : out Int) is - - XT : Int := 1; - YT : Int := 0; - - procedure Recur - (P, Q : Int; -- a (i-1), a (i) - X, Y : Int; -- x (i), y (i) - XP, YP : in out Int; -- x (i-1), y (i-1) - GCD : out Int); - - procedure Recur - (P, Q : Int; - X, Y : Int; - XP, YP : in out Int; - GCD : out Int) - is - Quo : Int := P / Q; -- q <-- |_ a (i-1) / a (i) _| - XT : Int := X; -- x (i) - YT : Int := 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; - - -- Start of processing for Euclid - - begin - Recur (P, Q, 0, 1, XT, YT, GCD); - X := XT; - Y := YT; - end Euclid; - - function Euclid (P, Q : Int) return Int is - X, Y, GCD : Int; - pragma Unreferenced (Y, GCD); - begin - Euclid (P, Q, X, Y, GCD); - return X; - end Euclid; - - ----------- - -- Image -- - ----------- - - function Image (Of_State : State) return String is - begin - return Int'Image (Of_State.X1) & ',' & Int'Image (Of_State.X2) - & ',' & - Int'Image (Of_State.P) & ',' & Int'Image (Of_State.Q); - end Image; + subtype Rep_Generator is System.Random_Numbers.Generator; + subtype Rep_State is System.Random_Numbers.State; ------------ -- Random -- ------------ - function Random (Gen : Generator) return Uniformly_Distributed is - Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; - + function Random (Gen : Generator) return Uniformly_Distributed is begin - Genp.X1 := Square_Mod_N (Genp.X1, Genp.P); - Genp.X2 := Square_Mod_N (Genp.X2, Genp.Q); - return - Float ((Flt (((Genp.X2 - Genp.X1) * Genp.X) - mod Genp.Q) * Flt (Genp.P) - + Flt (Genp.X1)) * Genp.Scl); + return Random (Gen.Rep); end Random; ----------- @@ -155,157 +68,56 @@ package body Ada.Numerics.Float_Random is -- Version that works from given initiator value procedure Reset (Gen : Generator; Initiator : Integer) is - Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; - X1, X2 : Int; - + G : Rep_Generator renames Gen.Rep'Unrestricted_Access.all; begin - X1 := 2 + Int (Initiator) mod (K1 - 3); - X2 := 2 + Int (Initiator) mod (K2 - 3); - - -- Eliminate effects of small initiators - - for J in 1 .. 5 loop - X1 := Square_Mod_N (X1, K1); - X2 := Square_Mod_N (X2, K2); - end loop; - - Genp.all := - (X1 => X1, - X2 => X2, - P => K1, - Q => K2, - X => 1, - Scl => Scal); - end Reset; - - -- Version that works from specific saved state - - procedure Reset (Gen : Generator; From_State : State) is - Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; - - begin - Genp.all := From_State; + Reset (G, Integer_32 (Initiator)); end Reset; -- Version that works from calendar procedure Reset (Gen : Generator) is - Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; - Now : constant Calendar.Time := Calendar.Clock; - X1, X2 : Int; - + G : Rep_Generator renames Gen.Rep'Unrestricted_Access.all; begin - X1 := Int (Calendar.Year (Now)) * 12 * 31 + - Int (Calendar.Month (Now)) * 31 + - Int (Calendar.Day (Now)); - - X2 := Int (Calendar.Seconds (Now) * Duration (1000.0)); - - X1 := 2 + X1 mod (K1 - 3); - X2 := 2 + X2 mod (K2 - 3); - - -- Eliminate visible effects of same day starts - - for J in 1 .. 5 loop - X1 := Square_Mod_N (X1, K1); - X2 := Square_Mod_N (X2, K2); - end loop; + Reset (G); + end Reset; - Genp.all := - (X1 => X1, - X2 => X2, - P => K1, - Q => K2, - X => 1, - Scl => Scal); + -- Version that works from specific saved state + procedure Reset (Gen : Generator; From_State : State) is + G : Rep_Generator renames Gen.Rep'Unrestricted_Access.all; + begin + Reset (G, From_State); end Reset; ---------- -- Save -- ---------- - procedure Save (Gen : Generator; To_State : out State) is + procedure Save (Gen : Generator; To_State : out State) is begin - To_State := Gen.Gen_State; + Save (Gen.Rep, State (To_State)); end Save; - ------------------ - -- Square_Mod_N -- - ------------------ - - function Square_Mod_N (X, N : Int) return Int is - Temp : constant Flt := Flt (X) * Flt (X); - Div : Int; + ----------- + -- Image -- + ----------- + function Image (Of_State : State) return String is begin - Div := Int (Temp / Flt (N)); - Div := Int (Temp - Flt (Div) * Flt (N)); - - if Div < 0 then - return Div + N; - else - return Div; - end if; - end Square_Mod_N; + return Image (Rep_State (Of_State)); + end Image; ----------- -- Value -- ----------- function Value (Coded_State : String) return State is - Last : constant Natural := Coded_State'Last; - Start : Positive := Coded_State'First; - Stop : Positive := Coded_State'First; - Outs : State; - + G : Generator; + S : Rep_State; begin - while Stop <= Last and then Coded_State (Stop) /= ',' loop - Stop := Stop + 1; - end loop; - - if Stop > Last then - raise Constraint_Error; - end if; - - Outs.X1 := Int'Value (Coded_State (Start .. Stop - 1)); - Start := Stop + 1; - - loop - Stop := Stop + 1; - exit when Stop > Last or else Coded_State (Stop) = ','; - end loop; - - if Stop > Last then - raise Constraint_Error; - end if; - - Outs.X2 := Int'Value (Coded_State (Start .. Stop - 1)); - Start := Stop + 1; - - loop - Stop := Stop + 1; - exit when Stop > Last or else Coded_State (Stop) = ','; - end loop; - - if Stop > Last then - raise Constraint_Error; - end if; - - Outs.P := Int'Value (Coded_State (Start .. Stop - 1)); - Outs.Q := Int'Value (Coded_State (Stop + 1 .. Last)); - Outs.X := Euclid (Outs.P, Outs.Q); - Outs.Scl := 1.0 / (Flt (Outs.P) * Flt (Outs.Q)); - - -- Now do *some* sanity checks - - if Outs.Q < 31 or else Outs.P < 31 - or else Outs.X1 not in 2 .. Outs.P - 1 - or else Outs.X2 not in 2 .. Outs.Q - 1 - then - raise Constraint_Error; - end if; - - return Outs; + Reset (G.Rep, Coded_State); + System.Random_Numbers.Save (G.Rep, S); + return State (S); end Value; + end Ada.Numerics.Float_Random; diff --git a/gcc/ada/a-nuflra.ads b/gcc/ada/a-nuflra.ads index e81842e23d8..9f8308121bb 100644 --- a/gcc/ada/a-nuflra.ads +++ b/gcc/ada/a-nuflra.ads @@ -6,7 +6,7 @@ -- -- -- S p e c -- -- -- --- Copyright (C) 1992-2009, Free Software Foundation, Inc. -- +-- Copyright (C) 1992-2010, Free Software Foundation, Inc. -- -- -- -- This specification is derived from the Ada Reference Manual for use with -- -- GNAT. The copyright notice above, and the license provisions that follow -- @@ -33,17 +33,10 @@ -- -- ------------------------------------------------------------------------------ --- Note: the implementation used in this package was contributed by --- Robert Eachus. It is based on the work of L. Blum, M. Blum, and --- M. Shub, SIAM Journal of Computing, Vol 15. No 2, May 1986. The --- particular choices for P and Q chosen here guarantee a period of --- 562,085,314,430,582 (about 2**49), and the generated sequence has --- excellent randomness properties. For further details, see the --- paper "Fast Generation of Trustworthy Random Numbers", by Robert --- Eachus, which describes both the algorithm and the efficient --- implementation approach used here. +-- Note: the implementation used in this package is a version of the +-- Mersenne Twister. See s-rannum.adb for details and references. -with Interfaces; +with System.Random_Numbers; package Ada.Numerics.Float_Random is @@ -65,35 +58,17 @@ package Ada.Numerics.Float_Random is procedure Save (Gen : Generator; To_State : out State); procedure Reset (Gen : Generator; From_State : State); - Max_Image_Width : constant := 80; + Max_Image_Width : constant := System.Random_Numbers.Max_Image_Width; function Image (Of_State : State) return String; function Value (Coded_State : String) return State; private - type Int is new Interfaces.Integer_32; - - -- We prefer to use 14 digits for Flt, but some targets are more limited - - type Flt is digits Positive'Min (14, Long_Long_Float'Digits); - - K1 : constant := 94_833_359; - K1F : constant := 94_833_359.0; - K2 : constant := 47_416_679; - K2F : constant := 47_416_679.0; - Scal : constant := 1.0 / (K1F * K2F); - - type State is record - X1 : Int := 2999 ** 2; -- Square mod p - X2 : Int := 1439 ** 2; -- Square mod q - P : Int := K1; - Q : Int := K2; - X : Int := 1; - Scl : Flt := Scal; - end record; type Generator is limited record - Gen_State : State; + Rep : System.Random_Numbers.Generator; end record; + type State is new System.Random_Numbers.State; + end Ada.Numerics.Float_Random; diff --git a/gcc/ada/g-mbdira.adb b/gcc/ada/g-mbdira.adb new file mode 100644 index 00000000000..20cb746d339 --- /dev/null +++ b/gcc/ada/g-mbdira.adb @@ -0,0 +1,299 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- G N A T . M B S S _ D I S C R E T E _ R A N D O M -- +-- -- +-- B o d y -- +-- -- +-- Copyright (C) 1992-2010, Free Software Foundation, Inc. -- +-- -- +-- GNAT is free software; you can redistribute it and/or modify it under -- +-- terms of the GNU General Public License as published by the Free Soft- -- +-- ware Foundation; either version 3, or (at your option) any later ver- -- +-- sion. GNAT is distributed in the hope that it will be useful, but WITH- -- +-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -- +-- or FITNESS FOR A PARTICULAR PURPOSE. -- +-- -- +-- As a special exception under Section 7 of GPL version 3, you are granted -- +-- additional permissions described in the GCC Runtime Library Exception, -- +-- version 3.1, as published by the Free Software Foundation. -- +-- -- +-- You should have received a copy of the GNU General Public License and -- +-- a copy of the GCC Runtime Library Exception along with this program; -- +-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see -- +-- . -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Calendar; + +with Interfaces; use Interfaces; + +package body GNAT.MBBS_Discrete_Random is + + ------------------------- + -- Implementation Note -- + ------------------------- + + -- The design of this spec is very awkward, as a result of Ada 95 not + -- permitting in-out parameters for function formals (most naturally + -- Generator values would be passed this way). In pure Ada 95, the only + -- solution is to use the heap and pointers, and, to avoid memory leaks, + -- controlled types. + + -- This is awfully heavy, so what we do is to use Unrestricted_Access to + -- get a pointer to the state in the passed Generator. This works because + -- Generator is a limited type and will thus always be passed by reference. + + package Calendar renames Ada.Calendar; + + type Pointer is access all State; + + Fits_In_32_Bits : constant Boolean := + Rst'Size < 31 + or else (Rst'Size = 31 + and then Rst'Pos (Rst'First) < 0); + -- This is set True if we do not need more than 32 bits in the result. If + -- we need 64-bits, we will only use the meaningful 48 bits of any 64-bit + -- number generated, since if more than 48 bits are required, we split the + -- computation into two separate parts, since the algorithm does not behave + -- above 48 bits. + + -- The way this expression works is that obviously if the size is 31 bits, + -- it fits in 32 bits. In the 32-bit case, it fits in 32-bit signed if the + -- range has negative values. It is too conservative in the case that the + -- programmer has set a size greater than the default, e.g. a size of 33 + -- for an integer type with a range of 1..10, but an over-conservative + -- result is OK. The important thing is that the value is only True if + -- we know the result will fit in 32-bits signed. If the value is False + -- when it could be True, the behavior will be correct, just a bit less + -- efficient than it could have been in some unusual cases. + -- + -- One might assume that we could get a more accurate result by testing + -- the lower and upper bounds of the type Rst against the bounds of 32-bit + -- Integer. However, there is no easy way to do that. Why? Because in the + -- relatively rare case where this expresion has to be evaluated at run + -- time rather than compile time (when the bounds are dynamic), we need a + -- type to use for the computation. But the possible range of upper bound + -- values for Rst (remembering the possibility of 64-bit modular types) is + -- from -2**63 to 2**64-1, and no run-time type has a big enough range. + + ----------------------- + -- Local Subprograms -- + ----------------------- + + function Square_Mod_N (X, N : Int) return Int; + pragma Inline (Square_Mod_N); + -- Computes X**2 mod N avoiding intermediate overflow + + ----------- + -- Image -- + ----------- + + function Image (Of_State : State) return String is + begin + return Int'Image (Of_State.X1) & + ',' & + Int'Image (Of_State.X2) & + ',' & + Int'Image (Of_State.Q); + end Image; + + ------------ + -- Random -- + ------------ + + function Random (Gen : Generator) return Rst is + Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; + Temp : Int; + TF : Flt; + + begin + -- Check for flat range here, since we are typically run with checks + -- off, note that in practice, this condition will usually be static + -- so we will not actually generate any code for the normal case. + + if Rst'Last < Rst'First then + raise Constraint_Error; + end if; + + -- Continue with computation if non-flat range + + Genp.X1 := Square_Mod_N (Genp.X1, Genp.P); + Genp.X2 := Square_Mod_N (Genp.X2, Genp.Q); + Temp := Genp.X2 - Genp.X1; + + -- Following duplication is not an error, it is a loop unwinding! + + if Temp < 0 then + Temp := Temp + Genp.Q; + end if; + + if Temp < 0 then + Temp := Temp + Genp.Q; + end if; + + TF := Offs + (Flt (Temp) * Flt (Genp.P) + Flt (Genp.X1)) * Genp.Scl; + + -- 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. + + if TF >= Flt (Rst'Pos (Rst'Last)) + 0.5 then + return Rst'First; + + elsif not Fits_In_32_Bits then + return Rst'Val (Interfaces.Integer_64 (TF)); + + else + return Rst'Val (Int (TF)); + end if; + end Random; + + ----------- + -- Reset -- + ----------- + + procedure Reset (Gen : Generator; Initiator : Integer) is + Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; + X1, X2 : Int; + + begin + X1 := 2 + Int (Initiator) mod (K1 - 3); + X2 := 2 + Int (Initiator) mod (K2 - 3); + + for J in 1 .. 5 loop + X1 := Square_Mod_N (X1, K1); + X2 := Square_Mod_N (X2, K2); + end loop; + + -- Eliminate effects of small Initiators + + Genp.all := + (X1 => X1, + X2 => X2, + P => K1, + Q => K2, + FP => K1F, + Scl => Scal); + end Reset; + + ----------- + -- Reset -- + ----------- + + procedure Reset (Gen : Generator) is + Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; + Now : constant Calendar.Time := Calendar.Clock; + X1 : Int; + X2 : Int; + + begin + X1 := Int (Calendar.Year (Now)) * 12 * 31 + + Int (Calendar.Month (Now) * 31) + + Int (Calendar.Day (Now)); + + X2 := Int (Calendar.Seconds (Now) * Duration (1000.0)); + + X1 := 2 + X1 mod (K1 - 3); + X2 := 2 + X2 mod (K2 - 3); + + -- Eliminate visible effects of same day starts + + for J in 1 .. 5 loop + X1 := Square_Mod_N (X1, K1); + X2 := Square_Mod_N (X2, K2); + end loop; + + Genp.all := + (X1 => X1, + X2 => X2, + P => K1, + Q => K2, + FP => K1F, + Scl => Scal); + + end Reset; + + ----------- + -- Reset -- + ----------- + + procedure Reset (Gen : Generator; From_State : State) is + Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; + begin + Genp.all := From_State; + end Reset; + + ---------- + -- Save -- + ---------- + + procedure Save (Gen : Generator; To_State : out State) is + begin + To_State := Gen.Gen_State; + end Save; + + ------------------ + -- Square_Mod_N -- + ------------------ + + function Square_Mod_N (X, N : Int) return Int is + begin + return Int ((Integer_64 (X) ** 2) mod (Integer_64 (N))); + end Square_Mod_N; + + ----------- + -- Value -- + ----------- + + function Value (Coded_State : String) return State is + Last : constant Natural := Coded_State'Last; + Start : Positive := Coded_State'First; + Stop : Positive := Coded_State'First; + Outs : State; + + begin + while Stop <= Last and then Coded_State (Stop) /= ',' loop + Stop := Stop + 1; + end loop; + + if Stop > Last then + raise Constraint_Error; + end if; + + Outs.X1 := Int'Value (Coded_State (Start .. Stop - 1)); + Start := Stop + 1; + + loop + Stop := Stop + 1; + exit when Stop > Last or else Coded_State (Stop) = ','; + end loop; + + if Stop > Last then + raise Constraint_Error; + end if; + + Outs.X2 := Int'Value (Coded_State (Start .. Stop - 1)); + Outs.Q := Int'Value (Coded_State (Stop + 1 .. Last)); + Outs.P := Outs.Q * 2 + 1; + Outs.FP := Flt (Outs.P); + Outs.Scl := (RstL - RstF + 1.0) / (Flt (Outs.P) * Flt (Outs.Q)); + + -- Now do *some* sanity checks + + if Outs.Q < 31 + or else Outs.X1 not in 2 .. Outs.P - 1 + or else Outs.X2 not in 2 .. Outs.Q - 1 + then + raise Constraint_Error; + end if; + + return Outs; + end Value; + +end GNAT.MBBS_Discrete_Random; diff --git a/gcc/ada/g-mbdira.ads b/gcc/ada/g-mbdira.ads new file mode 100644 index 00000000000..5c614e4ec42 --- /dev/null +++ b/gcc/ada/g-mbdira.ads @@ -0,0 +1,118 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- G N A T . M B S S _ D I S C R E T E _ R A N D O M -- +-- -- +-- S p e c -- +-- -- +-- Copyright (C) 1992-2010, Free Software Foundation, Inc. -- +-- -- +-- This specification is derived from the Ada Reference Manual for use with -- +-- GNAT. The copyright notice above, and the license provisions that follow -- +-- apply solely to the contents of the part following the private keyword. -- +-- -- +-- GNAT is free software; you can redistribute it and/or modify it under -- +-- terms of the GNU General Public License as published by the Free Soft- -- +-- ware Foundation; either version 3, or (at your option) any later ver- -- +-- sion. GNAT is distributed in the hope that it will be useful, but WITH- -- +-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -- +-- or FITNESS FOR A PARTICULAR PURPOSE. -- +-- -- +-- As a special exception under Section 7 of GPL version 3, you are granted -- +-- additional permissions described in the GCC Runtime Library Exception, -- +-- version 3.1, as published by the Free Software Foundation. -- +-- -- +-- You should have received a copy of the GNU General Public License and -- +-- a copy of the GCC Runtime Library Exception along with this program; -- +-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see -- +-- . -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +-- The implementation used in this package was contributed by Robert +-- Eachus. It is based on the work of L. Blum, M. Blum, and M. Shub, SIAM +-- Journal of Computing, Vol 15. No 2, May 1986. The particular choices for P +-- and Q chosen here guarantee a period of 562,085,314,430,582 (about 2**49), +-- and the generated sequence has excellent randomness properties. For further +-- details, see the paper "Fast Generation of Trustworthy Random Numbers", by +-- Robert Eachus, which describes both the algorithm and the efficient +-- implementation approach used here. + +-- Formerly, this package was Ada.Numerics.Discrete_Random. It is retained +-- here in part to allow users to reconstruct number sequences generated +-- by previous versions. + +with Interfaces; + +generic + type Result_Subtype is (<>); + +package GNAT.MBBS_Discrete_Random is + + -- The algorithm used here is reliable from a required statistical point of + -- view only up to 48 bits. We try to behave reasonably in the case of + -- larger types, but we can't guarantee the required properties. So + -- generate a warning for these (slightly) dubious cases. + + pragma Compile_Time_Warning + (Result_Subtype'Size > 48, + "statistical properties not guaranteed for size > 48"); + + -- Basic facilities + + type Generator is limited private; + + function Random (Gen : Generator) return Result_Subtype; + + procedure Reset (Gen : Generator); + procedure Reset (Gen : Generator; Initiator : Integer); + + -- Advanced facilities + + type State is private; + + procedure Save (Gen : Generator; To_State : out State); + procedure Reset (Gen : Generator; From_State : State); + + Max_Image_Width : constant := 80; + + function Image (Of_State : State) return String; + function Value (Coded_State : String) return State; + +private + subtype Int is Interfaces.Integer_32; + subtype Rst is Result_Subtype; + + -- We prefer to use 14 digits for Flt, but some targets are more limited + + type Flt is digits Positive'Min (14, Long_Long_Float'Digits); + + RstF : constant Flt := Flt (Rst'Pos (Rst'First)); + RstL : constant Flt := Flt (Rst'Pos (Rst'Last)); + + Offs : constant Flt := RstF - 0.5; + + K1 : constant := 94_833_359; + K1F : constant := 94_833_359.0; + K2 : constant := 47_416_679; + K2F : constant := 47_416_679.0; + Scal : constant Flt := (RstL - RstF + 1.0) / (K1F * K2F); + + type State is record + X1 : Int := Int (2999 ** 2); + X2 : Int := Int (1439 ** 2); + P : Int := K1; + Q : Int := K2; + FP : Flt := K1F; + Scl : Flt := Scal; + end record; + + type Generator is limited record + Gen_State : State; + end record; + +end GNAT.MBBS_Discrete_Random; diff --git a/gcc/ada/g-mbflra.adb b/gcc/ada/g-mbflra.adb new file mode 100644 index 00000000000..cf455707f64 --- /dev/null +++ b/gcc/ada/g-mbflra.adb @@ -0,0 +1,313 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- G N A T . M B S S _ F L O A T _ R A N D O M -- +-- -- +-- B o d y -- +-- -- +-- Copyright (C) 1992-2010, Free Software Foundation, Inc. -- +-- -- +-- GNAT is free software; you can redistribute it and/or modify it under -- +-- terms of the GNU General Public License as published by the Free Soft- -- +-- ware Foundation; either version 3, or (at your option) any later ver- -- +-- sion. GNAT is distributed in the hope that it will be useful, but WITH- -- +-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -- +-- or FITNESS FOR A PARTICULAR PURPOSE. -- +-- -- +-- As a special exception under Section 7 of GPL version 3, you are granted -- +-- additional permissions described in the GCC Runtime Library Exception, -- +-- version 3.1, as published by the Free Software Foundation. -- +-- -- +-- You should have received a copy of the GNU General Public License and -- +-- a copy of the GCC Runtime Library Exception along with this program; -- +-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see -- +-- . -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +with Ada.Calendar; + +package body GNAT.MBBS_Float_Random is + + ------------------------- + -- Implementation Note -- + ------------------------- + + -- The design of this spec is very awkward, as a result of Ada 95 not + -- permitting in-out parameters for function formals (most naturally + -- Generator values would be passed this way). In pure Ada 95, the only + -- solution is to use the heap and pointers, and, to avoid memory leaks, + -- controlled types. + + -- This is awfully heavy, so what we do is to use Unrestricted_Access to + -- get a pointer to the state in the passed Generator. This works because + -- Generator is a limited type and will thus always be passed by reference. + + package Calendar renames Ada.Calendar; + + type Pointer is access all State; + + ----------------------- + -- Local Subprograms -- + ----------------------- + + procedure Euclid (P, Q : Int; X, Y : out Int; GCD : out Int); + + function Euclid (P, Q : Int) return Int; + + function Square_Mod_N (X, N : Int) return Int; + + ------------ + -- Euclid -- + ------------ + + procedure Euclid (P, Q : Int; X, Y : out Int; GCD : out Int) is + + XT : Int := 1; + YT : Int := 0; + + procedure Recur + (P, Q : Int; -- a (i-1), a (i) + X, Y : Int; -- x (i), y (i) + XP, YP : in out Int; -- x (i-1), y (i-1) + GCD : out Int); + + procedure Recur + (P, Q : Int; + X, Y : Int; + XP, YP : in out Int; + GCD : out Int) + is + Quo : Int := P / Q; -- q <-- |_ a (i-1) / a (i) _| + XT : Int := X; -- x (i) + YT : Int := 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; + + -- Start of processing for Euclid + + begin + Recur (P, Q, 0, 1, XT, YT, GCD); + X := XT; + Y := YT; + end Euclid; + + function Euclid (P, Q : Int) return Int is + X, Y, GCD : Int; + pragma Unreferenced (Y, GCD); + begin + Euclid (P, Q, X, Y, GCD); + return X; + end Euclid; + + ----------- + -- Image -- + ----------- + + function Image (Of_State : State) return String is + begin + return Int'Image (Of_State.X1) & ',' & Int'Image (Of_State.X2) + & ',' & + Int'Image (Of_State.P) & ',' & Int'Image (Of_State.Q); + end Image; + + ------------ + -- Random -- + ------------ + + function Random (Gen : Generator) return Uniformly_Distributed is + Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; + + begin + Genp.X1 := Square_Mod_N (Genp.X1, Genp.P); + Genp.X2 := Square_Mod_N (Genp.X2, Genp.Q); + return + Float ((Flt (((Genp.X2 - Genp.X1) * Genp.X) + mod Genp.Q) * Flt (Genp.P) + + Flt (Genp.X1)) * Genp.Scl); + end Random; + + ----------- + -- Reset -- + ----------- + + -- Version that works from given initiator value + + procedure Reset (Gen : Generator; Initiator : Integer) is + Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; + X1, X2 : Int; + + begin + X1 := 2 + Int (Initiator) mod (K1 - 3); + X2 := 2 + Int (Initiator) mod (K2 - 3); + + -- Eliminate effects of small initiators + + for J in 1 .. 5 loop + X1 := Square_Mod_N (X1, K1); + X2 := Square_Mod_N (X2, K2); + end loop; + + Genp.all := + (X1 => X1, + X2 => X2, + P => K1, + Q => K2, + X => 1, + Scl => Scal); + end Reset; + + -- Version that works from specific saved state + + procedure Reset (Gen : Generator; From_State : State) is + Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; + + begin + Genp.all := From_State; + end Reset; + + -- Version that works from calendar + + procedure Reset (Gen : Generator) is + Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; + Now : constant Calendar.Time := Calendar.Clock; + X1, X2 : Int; + + begin + X1 := Int (Calendar.Year (Now)) * 12 * 31 + + Int (Calendar.Month (Now)) * 31 + + Int (Calendar.Day (Now)); + + X2 := Int (Calendar.Seconds (Now) * Duration (1000.0)); + + X1 := 2 + X1 mod (K1 - 3); + X2 := 2 + X2 mod (K2 - 3); + + -- Eliminate visible effects of same day starts + + for J in 1 .. 5 loop + X1 := Square_Mod_N (X1, K1); + X2 := Square_Mod_N (X2, K2); + end loop; + + Genp.all := + (X1 => X1, + X2 => X2, + P => K1, + Q => K2, + X => 1, + Scl => Scal); + + end Reset; + + ---------- + -- Save -- + ---------- + + procedure Save (Gen : Generator; To_State : out State) is + begin + To_State := Gen.Gen_State; + end Save; + + ------------------ + -- Square_Mod_N -- + ------------------ + + function Square_Mod_N (X, N : Int) return Int is + Temp : constant Flt := Flt (X) * Flt (X); + Div : Int; + + begin + Div := Int (Temp / Flt (N)); + Div := Int (Temp - Flt (Div) * Flt (N)); + + if Div < 0 then + return Div + N; + else + return Div; + end if; + end Square_Mod_N; + + ----------- + -- Value -- + ----------- + + function Value (Coded_State : String) return State is + Last : constant Natural := Coded_State'Last; + Start : Positive := Coded_State'First; + Stop : Positive := Coded_State'First; + Outs : State; + + begin + while Stop <= Last and then Coded_State (Stop) /= ',' loop + Stop := Stop + 1; + end loop; + + if Stop > Last then + raise Constraint_Error; + end if; + + Outs.X1 := Int'Value (Coded_State (Start .. Stop - 1)); + Start := Stop + 1; + + loop + Stop := Stop + 1; + exit when Stop > Last or else Coded_State (Stop) = ','; + end loop; + + if Stop > Last then + raise Constraint_Error; + end if; + + Outs.X2 := Int'Value (Coded_State (Start .. Stop - 1)); + Start := Stop + 1; + + loop + Stop := Stop + 1; + exit when Stop > Last or else Coded_State (Stop) = ','; + end loop; + + if Stop > Last then + raise Constraint_Error; + end if; + + Outs.P := Int'Value (Coded_State (Start .. Stop - 1)); + Outs.Q := Int'Value (Coded_State (Stop + 1 .. Last)); + Outs.X := Euclid (Outs.P, Outs.Q); + Outs.Scl := 1.0 / (Flt (Outs.P) * Flt (Outs.Q)); + + -- Now do *some* sanity checks + + if Outs.Q < 31 or else Outs.P < 31 + or else Outs.X1 not in 2 .. Outs.P - 1 + or else Outs.X2 not in 2 .. Outs.Q - 1 + then + raise Constraint_Error; + end if; + + return Outs; + end Value; +end GNAT.MBBS_Float_Random; diff --git a/gcc/ada/g-mbflra.ads b/gcc/ada/g-mbflra.ads new file mode 100644 index 00000000000..f9ad5af4238 --- /dev/null +++ b/gcc/ada/g-mbflra.ads @@ -0,0 +1,103 @@ +------------------------------------------------------------------------------ +-- -- +-- GNAT RUN-TIME COMPONENTS -- +-- -- +-- G N A T . M B S S _ F L O A T _ R A N D O M -- +-- -- +-- S p e c -- +-- -- +-- Copyright (C) 1992-2010, Free Software Foundation, Inc. -- +-- -- +-- This specification is derived from the Ada Reference Manual for use with -- +-- GNAT. The copyright notice above, and the license provisions that follow -- +-- apply solely to the contents of the part following the private keyword. -- +-- -- +-- GNAT is free software; you can redistribute it and/or modify it under -- +-- terms of the GNU General Public License as published by the Free Soft- -- +-- ware Foundation; either version 3, or (at your option) any later ver- -- +-- sion. GNAT is distributed in the hope that it will be useful, but WITH- -- +-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -- +-- or FITNESS FOR A PARTICULAR PURPOSE. -- +-- -- +-- As a special exception under Section 7 of GPL version 3, you are granted -- +-- additional permissions described in the GCC Runtime Library Exception, -- +-- version 3.1, as published by the Free Software Foundation. -- +-- -- +-- You should have received a copy of the GNU General Public License and -- +-- a copy of the GCC Runtime Library Exception along with this program; -- +-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see -- +-- . -- +-- -- +-- GNAT was originally developed by the GNAT team at New York University. -- +-- Extensive contributions were provided by Ada Core Technologies Inc. -- +-- -- +------------------------------------------------------------------------------ + +-- The implementation used in this package was contributed by +-- Robert Eachus. It is based on the work of L. Blum, M. Blum, and +-- M. Shub, SIAM Journal of Computing, Vol 15. No 2, May 1986. The +-- particular choices for P and Q chosen here guarantee a period of +-- 562,085,314,430,582 (about 2**49), and the generated sequence has +-- excellent randomness properties. For further details, see the +-- paper "Fast Generation of Trustworthy Random Numbers", by Robert +-- Eachus, which describes both the algorithm and the efficient +-- implementation approach used here. + +-- Formerly, this package was Ada.Numerics.Float_Random. It is retained +-- here in part to allow users to reconstruct number sequences generated +-- by previous versions. + +with Interfaces; + +package GNAT.MBBS_Float_Random is + + -- Basic facilities + + type Generator is limited private; + + subtype Uniformly_Distributed is Float range 0.0 .. 1.0; + + function Random (Gen : Generator) return Uniformly_Distributed; + + procedure Reset (Gen : Generator); + procedure Reset (Gen : Generator; Initiator : Integer); + + -- Advanced facilities + + type State is private; + + procedure Save (Gen : Generator; To_State : out State); + procedure Reset (Gen : Generator; From_State : State); + + Max_Image_Width : constant := 80; + + function Image (Of_State : State) return String; + function Value (Coded_State : String) return State; + +private + type Int is new Interfaces.Integer_32; + + -- We prefer to use 14 digits for Flt, but some targets are more limited + + type Flt is digits Positive'Min (14, Long_Long_Float'Digits); + + K1 : constant := 94_833_359; + K1F : constant := 94_833_359.0; + K2 : constant := 47_416_679; + K2F : constant := 47_416_679.0; + Scal : constant := 1.0 / (K1F * K2F); + + type State is record + X1 : Int := 2999 ** 2; -- Square mod p + X2 : Int := 1439 ** 2; -- Square mod q + P : Int := K1; + Q : Int := K2; + X : Int := 1; + Scl : Flt := Scal; + end record; + + type Generator is limited record + Gen_State : State; + end record; + +end GNAT.MBBS_Float_Random; diff --git a/gcc/ada/gnat_rm.texi b/gcc/ada/gnat_rm.texi index 3d3adae14c3..f3089d44f6a 100644 --- a/gcc/ada/gnat_rm.texi +++ b/gcc/ada/gnat_rm.texi @@ -8848,7 +8848,7 @@ floating-point. @code{Numerics.Float_Random.Max_Image_Width}. See A.5.2(27). @end cartouche @noindent -Maximum image width is 649, see library file @file{a-numran.ads}. +Maximum image width is 6864, see library file @file{s-rannum.ads}. @sp 1 @cartouche @@ -8857,7 +8857,7 @@ Maximum image width is 649, see library file @file{a-numran.ads}. @code{Numerics.Discrete_Random.Max_Image_Width}. See A.5.2(27). @end cartouche @noindent -Maximum image width is 80, see library file @file{a-nudira.ads}. +Maximum image width is 6864, see library file @file{s-rannum.ads}. @sp 1 @cartouche @@ -8866,8 +8866,8 @@ Maximum image width is 80, see library file @file{a-nudira.ads}. A.5.2(32). @end cartouche @noindent -The algorithm is documented in the source files @file{a-numran.ads} and -@file{a-numran.adb}. +The algorithm is the Mersenne Twister, as documented in the source file +@file{s-rannum.adb}. @sp 1 @cartouche @@ -8876,7 +8876,9 @@ The algorithm is documented in the source files @file{a-numran.ads} and state. See A.5.2(38). @end cartouche @noindent -See the documentation contained in the file @file{a-numran.adb}. +The value returned by the Image function is the concatenation of +the fixed-width decimal representations of the 624 32-bit integers +of the state vector. @sp 1 @cartouche @@ -11839,12 +11841,12 @@ This is a predefined instantiation of build the type @code{Complex} and @code{Imaginary}. @item Ada.Numerics.Discrete_Random -This package provides a random number generator suitable for generating -random integer values from a specified range. +This generic package provides a random number generator suitable for generating +uniformly distributed values of a specified discrete subtype. @item Ada.Numerics.Float_Random This package provides a random number generator suitable for generating -uniformly distributed floating point values. +uniformly distributed floating point values in the unit interval. @item Ada.Numerics.Generic_Complex_Elementary_Functions This is a generic version of the package that provides the diff --git a/gcc/ada/impunit.adb b/gcc/ada/impunit.adb index 83710311cd2..47e5b16311e 100644 --- a/gcc/ada/impunit.adb +++ b/gcc/ada/impunit.adb @@ -250,6 +250,8 @@ package body Impunit is "g-io ", -- GNAT.IO "g-io_aux", -- GNAT.IO_Aux "g-locfil", -- GNAT.Lock_Files + "g-mbdira", -- GNAT.MBBS_Discrete_Random + "g-mbflra", -- GNAT.MBBS_Float_Random "g-md5 ", -- GNAT.MD5 "g-memdum", -- GNAT.Memory_Dump "g-moreex", -- GNAT.Most_Recent_Exception diff --git a/gcc/ada/s-rannum.adb b/gcc/ada/s-rannum.adb index c3865a69a33..c161b6ebbe0 100644 --- a/gcc/ada/s-rannum.adb +++ b/gcc/ada/s-rannum.adb @@ -188,35 +188,104 @@ package body System.Random_Numbers is return Y; end Random; - function Random (Gen : Generator) return Float is - - -- Note: The application of Float'Machine (...) is necessary to avoid - -- returning extra significand bits. Without it, the function's value - -- will change if it is spilled, for example, causing - -- gratuitous nondeterminism. + generic + type Unsigned is mod <>; + type Real is digits <>; + with function Shift_Right (Value : Unsigned; Amount : Natural) + return Unsigned is <>; + with function Random (G : Generator) return Unsigned is <>; + function Random_Float_Template (Gen : Generator) return Real; + pragma Inline (Random_Float_Template); + -- Template for a random-number generator implementation that delivers + -- values of type Real in the half-open range [0 .. 1), using values from + -- Gen, assuming that Unsigned is large enough to hold the bits of + -- a mantissa for type Real. + + function Random_Float_Template (Gen : Generator) return Real is + -- This code generates random floating-point numbers from unsigned + -- integers. Assuming that Real'Machine_Radix = 2, it can deliver + -- all machine values of type Real (at least as implied by + -- Real'Machine_Mantissa and Real'Machine_Emin), which is not true + -- of the standard method (to which we fall back for non-binary + -- radix): computing Real() / (+1). + -- To do so, we first extract an (M-1)-bit significand (where M + -- is Real'Machine_Mantissa), and then decide on a normalized + -- exponent by repeated coin flips, decrementing from 0 as long as + -- we flip heads (1 bits). This yields the proper geometric + -- distribution for the exponent: in a uniformly distributed set of + -- floating-point numbers, 1/2 of them will be in [0.5, 1), 1/4 will + -- be in [0.25, 0.5), and so forth. If the process reaches + -- Machine_Emin (an extremely rare event), it uses the selected + -- mantissa bits as an unnormalized fraction with Machine_Emin as + -- exponent. Otherwise, it adds a leading bit to the selected + -- mantissa bits (thus giving a normalized fraction) and adjusts by + -- the chosen exponent. The algorithm attempts to be stingy with + -- random integers. In the worst case, it can consume roughly + -- -Real'Machine_Emin/32 32-bit integers, but this case occurs with + -- probability 2**Machine_Emin, and the expected number of calls to + -- integer-valued Random is 1. - Result : constant Float := - Float'Machine - (Float (Unsigned_32'(Random (Gen))) * 2.0 ** (-32)); begin - if Result < 1.0 then - return Result; + if Real'Machine_Radix /= 2 then + declare + Val : constant Real := Real'Machine + (Real (Unsigned'(Random (Gen))) * 2.0**(-Unsigned'Size)); + begin + if Val < 1.0 then + return Real'Base (Val); + else + return Real'Pred (1.0); + end if; + end; else - return Float'Adjacent (1.0, 0.0); + declare + Mant_Bits : constant Integer := Real'Machine_Mantissa - 1; + Mant_Mask : constant Unsigned := 2**Mant_Bits - 1; + Adjust32 : constant Integer := Real'Size - Unsigned_32'Size; + Leftover : constant Integer := + Unsigned'Size - Real'Machine_Mantissa + 1; + + V : constant Unsigned := Random (Gen); + Mant : constant Unsigned := V and Mant_Mask; + Rand_Bits : Unsigned_32; + Exp : Integer; + Bits_Left : Integer; + Result : Real; + begin + Rand_Bits := Unsigned_32 (Shift_Right (V, Adjust32)); + Exp := 0; + Bits_Left := Leftover; + Result := Real (Mant + 2**Mant_Bits) * 2.0**(-Mant_Bits - 1); + while Rand_Bits >= 2**31 loop + if Exp = Real'Machine_Emin then + return Real (Mant) * 2.0**Real'Machine_Emin; + end if; + + Result := Result * 0.5; + Exp := Exp - 1; + Rand_Bits := 2 * Rand_Bits; + Bits_Left := Bits_Left - 1; + + if Bits_Left = 0 then + Bits_Left := 32; + Rand_Bits := Random (Gen); + end if; + end loop; + return Result; + end; end if; + end Random_Float_Template; + + function Random (Gen : Generator) return Float is + function F is new Random_Float_Template (Unsigned_32, Float); + begin + return F (Gen); end Random; function Random (Gen : Generator) return Long_Float is - Result : constant Long_Float := - Long_Float'Machine ((Long_Float (Unsigned_32'(Random (Gen))) - * 2.0 ** (-32)) - + (Long_Float (Unsigned_32'(Random (Gen))) * 2.0 ** (-64))); + function F is new Random_Float_Template (Unsigned_64, Long_Float); begin - if Result < 1.0 then - return Result; - else - return Long_Float'Adjacent (1.0, 0.0); - end if; + return F (Gen); end Random; function Random (Gen : Generator) return Unsigned_64 is -- 2.30.2