Commit 41195c94 by Paul Hilfinger Committed by Arnaud Charlet

2010-06-22 Paul Hilfinger <hilfinger@adacore.com>

	* 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
parent 07309d58
2010-06-22 Paul Hilfinger <hilfinger@adacore.com>
* 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 <dismukes@adacore.com> 2010-06-22 Gary Dismukes <dismukes@adacore.com>
* sem_ch5.adb (Analyze_Assignment): Revise test for illegal assignment * sem_ch5.adb (Analyze_Assignment): Revise test for illegal assignment
......
# Makefile.rtl for GNU Ada Compiler (GNAT). # 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. #This file is part of GCC.
...@@ -359,6 +359,8 @@ GNATRTL_NONTASKING_OBJS= \ ...@@ -359,6 +359,8 @@ GNATRTL_NONTASKING_OBJS= \
g-io$(objext) \ g-io$(objext) \
g-io_aux$(objext) \ g-io_aux$(objext) \
g-locfil$(objext) \ g-locfil$(objext) \
g-mbdira$(objext) \
g-mbflra$(objext) \
g-md5$(objext) \ g-md5$(objext) \
g-memdum$(objext) \ g-memdum$(objext) \
g-moreex$(objext) \ g-moreex$(objext) \
......
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
-- -- -- --
-- B o d y -- -- 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 -- -- 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- -- -- terms of the GNU General Public License as published by the Free Soft- --
...@@ -29,9 +29,7 @@ ...@@ -29,9 +29,7 @@
-- -- -- --
------------------------------------------------------------------------------ ------------------------------------------------------------------------------
with Ada.Calendar; with System.Random_Numbers; use System.Random_Numbers;
with Interfaces; use Interfaces;
package body Ada.Numerics.Discrete_Random is package body Ada.Numerics.Discrete_Random is
...@@ -49,249 +47,55 @@ 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 -- 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. -- Generator is a limited type and will thus always be passed by reference.
type Pointer is access all State; subtype Rep_Generator is System.Random_Numbers.Generator;
subtype Rep_State is System.Random_Numbers.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; function Rep_Random is
pragma Inline (Square_Mod_N); new Random_Discrete (Result_Subtype, Result_Subtype'First);
-- Computes X**2 mod N avoiding intermediate overflow
----------- function Random (Gen : Generator) return Result_Subtype is
-- Image --
-----------
function Image (Of_State : State) return String is
begin begin
return Int'Image (Of_State.X1) & return Rep_Random (Gen.Rep);
',' &
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; end Random;
----------- procedure Reset (Gen : Generator;
-- Reset -- Initiator : Integer) is
----------- G : Rep_Generator renames Gen.Rep'Unrestricted_Access.all;
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 begin
X1 := Int (Calendar.Year (Now)) * 12 * 31 + Reset (G, Initiator);
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; end Reset;
----------- procedure Reset (Gen : Generator) is
-- Reset -- G : Rep_Generator renames Gen.Rep'Unrestricted_Access.all;
-----------
procedure Reset (Gen : Generator; From_State : State) is
Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access;
begin begin
Genp.all := From_State; Reset (G);
end Reset; end Reset;
---------- procedure Save (Gen : Generator;
-- Save -- To_State : out State) is
----------
procedure Save (Gen : Generator; To_State : out State) is
begin begin
To_State := Gen.Gen_State; Save (Gen.Rep, State (To_State));
end Save; end Save;
------------------ procedure Reset (Gen : Generator;
-- Square_Mod_N -- From_State : State) is
------------------ G : Rep_Generator renames Gen.Rep'Unrestricted_Access.all;
function Square_Mod_N (X, N : Int) return Int is
begin begin
return Int ((Integer_64 (X) ** 2) mod (Integer_64 (N))); Reset (G, From_State);
end Square_Mod_N; end Reset;
----------- function Image (Of_State : State) return String is
-- Value -- begin
----------- return Image (Rep_State (Of_State));
end Image;
function Value (Coded_State : String) return State is function Value (Coded_State : String) return State is
Last : constant Natural := Coded_State'Last; G : Generator;
Start : Positive := Coded_State'First; S : Rep_State;
Stop : Positive := Coded_State'First;
Outs : State;
begin begin
while Stop <= Last and then Coded_State (Stop) /= ',' loop Reset (G.Rep, Coded_State);
Stop := Stop + 1; System.Random_Numbers.Save (G.Rep, S);
end loop; return State (S);
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 Value;
end Ada.Numerics.Discrete_Random; end Ada.Numerics.Discrete_Random;
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
-- -- -- --
-- S p e c -- -- 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 -- -- This specification is derived from the Ada Reference Manual for use with --
-- GNAT. The copyright notice above, and the license provisions that follow -- -- GNAT. The copyright notice above, and the license provisions that follow --
...@@ -33,39 +33,24 @@ ...@@ -33,39 +33,24 @@
-- -- -- --
------------------------------------------------------------------------------ ------------------------------------------------------------------------------
-- Note: the implementation used in this package was contributed by Robert -- Note: the implementation used in this package is a version of the
-- Eachus. It is based on the work of L. Blum, M. Blum, and M. Shub, SIAM -- Mersenne Twister. See s-rannum.adb for details and references.
-- 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.
with Interfaces; with System.Random_Numbers;
generic generic
type Result_Subtype is (<>); type Result_Subtype is (<>);
package Ada.Numerics.Discrete_Random 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 -- Basic facilities
type Generator is limited private; type Generator is limited private;
function Random (Gen : Generator) return Result_Subtype; function Random (Gen : Generator) return Result_Subtype;
procedure Reset (Gen : Generator);
procedure Reset (Gen : Generator; Initiator : Integer); procedure Reset (Gen : Generator; Initiator : Integer);
procedure Reset (Gen : Generator);
-- Advanced facilities -- Advanced facilities
...@@ -74,41 +59,17 @@ package Ada.Numerics.Discrete_Random is ...@@ -74,41 +59,17 @@ package Ada.Numerics.Discrete_Random is
procedure Save (Gen : Generator; To_State : out State); procedure Save (Gen : Generator; To_State : out State);
procedure Reset (Gen : Generator; From_State : 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 Image (Of_State : State) return String;
function Value (Coded_State : String) return State; function Value (Coded_State : String) return State;
private 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 type Generator is limited record
Gen_State : State; Rep : System.Random_Numbers.Generator;
end record; end record;
type State is new System.Random_Numbers.State;
end Ada.Numerics.Discrete_Random; end Ada.Numerics.Discrete_Random;
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
-- -- -- --
-- B o d y -- -- 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 -- -- 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- -- -- terms of the GNU General Public License as published by the Free Soft- --
...@@ -29,7 +29,9 @@ ...@@ -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 package body Ada.Numerics.Float_Random is
...@@ -47,105 +49,16 @@ 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 -- 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. -- Generator is a limited type and will thus always be passed by reference.
type Pointer is access all State; subtype Rep_Generator is System.Random_Numbers.Generator;
subtype Rep_State is System.Random_Numbers.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 -- -- Random --
------------ ------------
function Random (Gen : Generator) return Uniformly_Distributed is function Random (Gen : Generator) return Uniformly_Distributed is
Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access;
begin begin
Genp.X1 := Square_Mod_N (Genp.X1, Genp.P); return Random (Gen.Rep);
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; end Random;
----------- -----------
...@@ -155,157 +68,56 @@ package body Ada.Numerics.Float_Random is ...@@ -155,157 +68,56 @@ package body Ada.Numerics.Float_Random is
-- Version that works from given initiator value -- Version that works from given initiator value
procedure Reset (Gen : Generator; Initiator : Integer) is procedure Reset (Gen : Generator; Initiator : Integer) is
Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; G : Rep_Generator renames Gen.Rep'Unrestricted_Access.all;
X1, X2 : Int;
begin begin
X1 := 2 + Int (Initiator) mod (K1 - 3); Reset (G, Integer_32 (Initiator));
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; end Reset;
-- Version that works from calendar -- Version that works from calendar
procedure Reset (Gen : Generator) is procedure Reset (Gen : Generator) is
Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access; G : Rep_Generator renames Gen.Rep'Unrestricted_Access.all;
Now : constant Calendar.Time := Calendar.Clock;
X1, X2 : Int;
begin begin
X1 := Int (Calendar.Year (Now)) * 12 * 31 + Reset (G);
Int (Calendar.Month (Now)) * 31 + end Reset;
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 := -- Version that works from specific saved state
(X1 => X1,
X2 => X2,
P => K1,
Q => K2,
X => 1,
Scl => Scal);
procedure Reset (Gen : Generator; From_State : State) is
G : Rep_Generator renames Gen.Rep'Unrestricted_Access.all;
begin
Reset (G, From_State);
end Reset; end Reset;
---------- ----------
-- Save -- -- Save --
---------- ----------
procedure Save (Gen : Generator; To_State : out State) is procedure Save (Gen : Generator; To_State : out State) is
begin begin
To_State := Gen.Gen_State; Save (Gen.Rep, State (To_State));
end Save; end Save;
------------------ -----------
-- Square_Mod_N -- -- Image --
------------------ -----------
function Square_Mod_N (X, N : Int) return Int is
Temp : constant Flt := Flt (X) * Flt (X);
Div : Int;
function Image (Of_State : State) return String is
begin begin
Div := Int (Temp / Flt (N)); return Image (Rep_State (Of_State));
Div := Int (Temp - Flt (Div) * Flt (N)); end Image;
if Div < 0 then
return Div + N;
else
return Div;
end if;
end Square_Mod_N;
----------- -----------
-- Value -- -- Value --
----------- -----------
function Value (Coded_State : String) return State is function Value (Coded_State : String) return State is
Last : constant Natural := Coded_State'Last; G : Generator;
Start : Positive := Coded_State'First; S : Rep_State;
Stop : Positive := Coded_State'First;
Outs : State;
begin begin
while Stop <= Last and then Coded_State (Stop) /= ',' loop Reset (G.Rep, Coded_State);
Stop := Stop + 1; System.Random_Numbers.Save (G.Rep, S);
end loop; return State (S);
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 Value;
end Ada.Numerics.Float_Random; end Ada.Numerics.Float_Random;
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
-- -- -- --
-- S p e c -- -- 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 -- -- This specification is derived from the Ada Reference Manual for use with --
-- GNAT. The copyright notice above, and the license provisions that follow -- -- GNAT. The copyright notice above, and the license provisions that follow --
...@@ -33,17 +33,10 @@ ...@@ -33,17 +33,10 @@
-- -- -- --
------------------------------------------------------------------------------ ------------------------------------------------------------------------------
-- Note: the implementation used in this package was contributed by -- Note: the implementation used in this package is a version of the
-- Robert Eachus. It is based on the work of L. Blum, M. Blum, and -- Mersenne Twister. See s-rannum.adb for details and references.
-- 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.
with Interfaces; with System.Random_Numbers;
package Ada.Numerics.Float_Random is package Ada.Numerics.Float_Random is
...@@ -65,35 +58,17 @@ 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 Save (Gen : Generator; To_State : out State);
procedure Reset (Gen : Generator; From_State : 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 Image (Of_State : State) return String;
function Value (Coded_State : String) return State; function Value (Coded_State : String) return State;
private 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 type Generator is limited record
Gen_State : State; Rep : System.Random_Numbers.Generator;
end record; end record;
type State is new System.Random_Numbers.State;
end Ada.Numerics.Float_Random; end Ada.Numerics.Float_Random;
------------------------------------------------------------------------------
-- --
-- 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 --
-- <http://www.gnu.org/licenses/>. --
-- --
-- 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;
------------------------------------------------------------------------------
-- --
-- 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 --
-- <http://www.gnu.org/licenses/>. --
-- --
-- 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;
------------------------------------------------------------------------------
-- --
-- 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 --
-- <http://www.gnu.org/licenses/>. --
-- --
-- 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;
...@@ -8848,7 +8848,7 @@ floating-point. ...@@ -8848,7 +8848,7 @@ floating-point.
@code{Numerics.Float_Random.Max_Image_Width}. See A.5.2(27). @code{Numerics.Float_Random.Max_Image_Width}. See A.5.2(27).
@end cartouche @end cartouche
@noindent @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 @sp 1
@cartouche @cartouche
...@@ -8857,7 +8857,7 @@ Maximum image width is 649, see library file @file{a-numran.ads}. ...@@ -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). @code{Numerics.Discrete_Random.Max_Image_Width}. See A.5.2(27).
@end cartouche @end cartouche
@noindent @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 @sp 1
@cartouche @cartouche
...@@ -8866,8 +8866,8 @@ Maximum image width is 80, see library file @file{a-nudira.ads}. ...@@ -8866,8 +8866,8 @@ Maximum image width is 80, see library file @file{a-nudira.ads}.
A.5.2(32). A.5.2(32).
@end cartouche @end cartouche
@noindent @noindent
The algorithm is documented in the source files @file{a-numran.ads} and The algorithm is the Mersenne Twister, as documented in the source file
@file{a-numran.adb}. @file{s-rannum.adb}.
@sp 1 @sp 1
@cartouche @cartouche
...@@ -8876,7 +8876,9 @@ The algorithm is documented in the source files @file{a-numran.ads} and ...@@ -8876,7 +8876,9 @@ The algorithm is documented in the source files @file{a-numran.ads} and
state. See A.5.2(38). state. See A.5.2(38).
@end cartouche @end cartouche
@noindent @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 @sp 1
@cartouche @cartouche
...@@ -11839,12 +11841,12 @@ This is a predefined instantiation of ...@@ -11839,12 +11841,12 @@ This is a predefined instantiation of
build the type @code{Complex} and @code{Imaginary}. build the type @code{Complex} and @code{Imaginary}.
@item Ada.Numerics.Discrete_Random @item Ada.Numerics.Discrete_Random
This package provides a random number generator suitable for generating This generic package provides a random number generator suitable for generating
random integer values from a specified range. uniformly distributed values of a specified discrete subtype.
@item Ada.Numerics.Float_Random @item Ada.Numerics.Float_Random
This package provides a random number generator suitable for generating 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 @item Ada.Numerics.Generic_Complex_Elementary_Functions
This is a generic version of the package that provides the This is a generic version of the package that provides the
......
...@@ -250,6 +250,8 @@ package body Impunit is ...@@ -250,6 +250,8 @@ package body Impunit is
"g-io ", -- GNAT.IO "g-io ", -- GNAT.IO
"g-io_aux", -- GNAT.IO_Aux "g-io_aux", -- GNAT.IO_Aux
"g-locfil", -- GNAT.Lock_Files "g-locfil", -- GNAT.Lock_Files
"g-mbdira", -- GNAT.MBBS_Discrete_Random
"g-mbflra", -- GNAT.MBBS_Float_Random
"g-md5 ", -- GNAT.MD5 "g-md5 ", -- GNAT.MD5
"g-memdum", -- GNAT.Memory_Dump "g-memdum", -- GNAT.Memory_Dump
"g-moreex", -- GNAT.Most_Recent_Exception "g-moreex", -- GNAT.Most_Recent_Exception
......
...@@ -188,35 +188,104 @@ package body System.Random_Numbers is ...@@ -188,35 +188,104 @@ package body System.Random_Numbers is
return Y; return Y;
end Random; end Random;
function Random (Gen : Generator) return Float is generic
type Unsigned is mod <>;
-- Note: The application of Float'Machine (...) is necessary to avoid type Real is digits <>;
-- returning extra significand bits. Without it, the function's value with function Shift_Right (Value : Unsigned; Amount : Natural)
-- will change if it is spilled, for example, causing return Unsigned is <>;
-- gratuitous nondeterminism. 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(<random integer>) / (<max random integer>+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 begin
if Result < 1.0 then if Real'Machine_Radix /= 2 then
return Result; 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 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 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; end Random;
function Random (Gen : Generator) return Long_Float is function Random (Gen : Generator) return Long_Float is
Result : constant Long_Float := function F is new Random_Float_Template (Unsigned_64, Long_Float);
Long_Float'Machine ((Long_Float (Unsigned_32'(Random (Gen)))
* 2.0 ** (-32))
+ (Long_Float (Unsigned_32'(Random (Gen))) * 2.0 ** (-64)));
begin begin
if Result < 1.0 then return F (Gen);
return Result;
else
return Long_Float'Adjacent (1.0, 0.0);
end if;
end Random; end Random;
function Random (Gen : Generator) return Unsigned_64 is function Random (Gen : Generator) return Unsigned_64 is
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment