Commit 5bec9717 by Arnaud Charlet

[multiple changes]

2010-06-22  Matthew Heaney  <heaney@adacore.com>

	* a-convec.adb, a-coinve.adb: Removed 64-bit types Int and UInt.

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

	* s-rannum.adb (Random_Float_Template): Replace with unbiased version
	that is able to produce all representable floating-point numbers in the
	unit interval. Remove template parameter Shift_Right, no longer used.
	* gnat_rm.texi: Document the period of the pseudo-random number
	generator under the description of its algorithm.
	* gcc-interface/Make-lang.in: Update dependencies.

From-SVN: r161202
parent 5087048c
2010-06-22 Matthew Heaney <heaney@adacore.com>
* a-convec.adb, a-coinve.adb: Removed 64-bit types Int and UInt.
2010-06-22 Paul Hilfinger <hilfinger@adacore.com>
* s-rannum.adb (Random_Float_Template): Replace with unbiased version
that is able to produce all representable floating-point numbers in the
unit interval. Remove template parameter Shift_Right, no longer used.
* gnat_rm.texi: Document the period of the pseudo-random number
generator under the description of its algorithm.
* gcc-interface/Make-lang.in: Update dependencies.
2010-06-22 Thomas Quinot <quinot@adacore.com> 2010-06-22 Thomas Quinot <quinot@adacore.com>
* exp_aggr.adb (Rewrite_Discriminant): Fix predicate used to identify * exp_aggr.adb (Rewrite_Discriminant): Fix predicate used to identify
......
...@@ -8869,7 +8869,8 @@ A.5.2(32). ...@@ -8869,7 +8869,8 @@ A.5.2(32).
@end cartouche @end cartouche
@noindent @noindent
The algorithm is the Mersenne Twister, as documented in the source file The algorithm is the Mersenne Twister, as documented in the source file
@file{s-rannum.adb}. @file{s-rannum.adb}. This version of the algorithm has a period of
2**19937-1.
@sp 1 @sp 1
@cartouche @cartouche
......
...@@ -191,17 +191,21 @@ package body System.Random_Numbers is ...@@ -191,17 +191,21 @@ package body System.Random_Numbers is
generic generic
type Unsigned is mod <>; type Unsigned is mod <>;
type Real is digits <>; type Real is digits <>;
with function Shift_Right (Value : Unsigned; Amount : Natural)
return Unsigned is <>;
with function Random (G : Generator) return Unsigned is <>; with function Random (G : Generator) return Unsigned is <>;
function Random_Float_Template (Gen : Generator) return Real; function Random_Float_Template (Gen : Generator) return Real;
pragma Inline (Random_Float_Template); pragma Inline (Random_Float_Template);
-- Template for a random-number generator implementation that delivers -- Template for a random-number generator implementation that delivers
-- values of type Real in the half-open range [0 .. 1), using values from -- values of type Real in the range [0 .. 1], using values from Gen,
-- Gen, assuming that Unsigned is large enough to hold the bits of -- assuming that Unsigned is large enough to hold the bits of a mantissa
-- a mantissa for type Real. -- for type Real.
function Random_Float_Template (Gen : Generator) return Real is function Random_Float_Template (Gen : Generator) return Real is
pragma Compile_Time_Error
(Unsigned'Last <= 2**(Real'Machine_Mantissa - 1),
"insufficiently large modular type used to hold mantissa");
begin
-- This code generates random floating-point numbers from unsigned -- This code generates random floating-point numbers from unsigned
-- integers. Assuming that Real'Machine_Radix = 2, it can deliver all -- integers. Assuming that Real'Machine_Radix = 2, it can deliver all
-- machine values of type Real (as implied by Real'Machine_Mantissa and -- machine values of type Real (as implied by Real'Machine_Mantissa and
...@@ -210,69 +214,118 @@ package body System.Random_Numbers is ...@@ -210,69 +214,118 @@ package body System.Random_Numbers is
-- integer>) / (<max random integer>+1). To do so, we first extract an -- integer>) / (<max random integer>+1). To do so, we first extract an
-- (M-1)-bit significand (where M is Real'Machine_Mantissa), and then -- (M-1)-bit significand (where M is Real'Machine_Mantissa), and then
-- decide on a normalized exponent by repeated coin flips, decrementing -- decide on a normalized exponent by repeated coin flips, decrementing
-- from 0 as long as we flip heads (1 bits). This yields the proper -- from 0 as long as we flip heads (1 bits). This process yields the
-- geometric distribution for the exponent: in a uniformly distributed -- proper geometric distribution for the exponent: in a uniformly
-- set of floating-point numbers, 1/2 of them will be in [0.5, 1), 1/4 -- distributed set of floating-point numbers, 1/2 of them will be in
-- will be in [0.25, 0.5), and so forth. If the process reaches -- (0.5, 1], 1/4 will be in (0.25, 0.5], and so forth. It makes a
-- Machine_Emin (an extremely rare event), it uses the selected mantissa -- further adjustment at binade boundaries (see comments below) to give
-- bits as an unnormalized fraction with Machine_Emin as exponent. -- the effect of selecting a uniformly distributed real deviate in
-- Otherwise, it adds a leading bit to the selected mantissa bits (thus -- [0..1] and then rounding to the nearest representable floating-point
-- giving a normalized fraction) and adjusts by the chosen exponent. The -- number. The algorithm attempts to be stingy with random integers. In
-- algorithm attempts to be stingy with random integers. In the worst -- the worst case, it can consume roughly -Real'Machine_Emin/32 32-bit
-- case, it can consume roughly -Real'Machine_Emin/32 32-bit integers, -- integers, but this case occurs with probability around
-- but this case occurs with probability 2**Machine_Emin, and the -- 2**Machine_Emin, and the expected number of calls to integer-valued
-- expected number of calls to integer-valued Random is 1. -- Random is 1. For another discussion of the issues addressed by this
-- process, see Allen Downey's unpublished paper at
-- http://allendowney.com/research/rand/downey07randfloat.pdf.
begin
if Real'Machine_Radix /= 2 then if Real'Machine_Radix /= 2 then
return Real'Machine
(Real (Unsigned'(Random (Gen))) * 2.0**(-Unsigned'Size));
else
declare declare
Val : constant Real := type Bit_Count is range 0 .. 4;
Real'Machine
(Real (Unsigned'(Random (Gen))) * 2.0**(-Unsigned'Size)); subtype T is Real'Base;
Trailing_Ones : constant array (Unsigned_32 range 0 .. 15)
of Bit_Count
:= (2#00000# => 0, 2#00001# => 1, 2#00010# => 0, 2#00011# => 2,
2#00100# => 0, 2#00101# => 1, 2#00110# => 0, 2#00111# => 3,
2#01000# => 0, 2#01001# => 1, 2#01010# => 0, 2#01011# => 2,
2#01100# => 0, 2#01101# => 1, 2#01110# => 0, 2#01111# => 4);
Pow_Tab : constant array (Bit_Count range 0 .. 3) of Real
:= (0 => 2.0**(0 - T'Machine_Mantissa),
1 => 2.0**(-1 - T'Machine_Mantissa),
2 => 2.0**(-2 - T'Machine_Mantissa),
3 => 2.0**(-3 - T'Machine_Mantissa));
Extra_Bits : constant Natural :=
(Unsigned'Size - T'Machine_Mantissa + 1);
-- Random bits left over after selecting mantissa
Mantissa : Unsigned;
X : Real; -- Scaled mantissa
R : Unsigned_32; -- Supply of random bits
R_Bits : Natural; -- Number of bits left in R
K : Bit_Count; -- Next decrement to exponent
begin begin
if Val < 1.0 then
return Real'Base (Val); Mantissa := Random (Gen) / 2**Extra_Bits;
R := Unsigned_32 (Mantissa mod 2**Extra_Bits);
R_Bits := Extra_Bits;
X := Real (2**(T'Machine_Mantissa - 1) + Mantissa); -- Exact
if Extra_Bits < 4 and then R < 2**Extra_Bits - 1 then
-- We got lucky and got a zero in our few extra bits
K := Trailing_Ones (R);
else else
return Real'Pred (1.0); Find_Zero : loop
-- R has R_Bits unprocessed random bits, a multiple of 4.
-- X needs to be halved for each trailing one bit. The
-- process stops as soon as a 0 bit is found. If R_Bits
-- becomes zero, reload R.
-- Process 4 bits at a time for speed: the two iterations
-- on average with three tests each was still too slow,
-- probably because the branches are not predictable.
-- This loop now will only execute once 94% of the cases,
-- doing more bits at a time will not help.
while R_Bits >= 4 loop
K := Trailing_Ones (R mod 16);
exit Find_Zero when K < 4; -- Exits 94% of the time
R_Bits := R_Bits - 4;
X := X / 16.0;
R := R / 16;
end loop;
-- Do not allow us to loop endlessly even in the (very
-- unlikely) case that Random (Gen) keeps yielding all ones.
exit Find_Zero when X = 0.0;
R := Random (Gen);
R_Bits := 32;
end loop Find_Zero;
end if; end if;
end;
else -- K has the count of trailing ones not reflected yet in X.
declare -- The following multiplication takes care of that, as well
Mant_Bits : constant Integer := Real'Machine_Mantissa - 1; -- as the correction to move the radix point to the left of
Mant_Mask : constant Unsigned := 2**Mant_Bits - 1; -- the mantissa. Doing it at the end avoids repeated rounding
Adjust32 : constant Integer := Real'Size - Unsigned_32'Size; -- errors in the exceedingly unlikely case of ever having
Leftover : constant Integer := -- a subnormal result.
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 X := X * Pow_Tab (K);
Rand_Bits := Unsigned_32 (Shift_Right (V, Adjust32));
Exp := 0; -- The smallest value in each binade is rounded to by 0.75 of
Bits_Left := Leftover; -- the span of real numbers as its next larger neighbor, and
Result := Real (Mant + 2**Mant_Bits) * 2.0**(-Mant_Bits - 1); -- 1.0 is rounded to by half of the span of real numbers as its
while Rand_Bits >= 2**31 loop -- next smaller neighbor. To account for this, when we encounter
if Exp = Real'Machine_Emin then -- the smallest number in a binade, we substitute the smallest
return Real (Mant) * 2.0**Real'Machine_Emin; -- value in the next larger binade with probability 1/2.
end if;
if Mantissa = 0 and then Unsigned_32'(Random (Gen)) mod 2 = 0 then
Result := Result * 0.5; X := 2.0 * X;
Exp := Exp - 1; end if;
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; return X;
end; end;
end if; end if;
end Random_Float_Template; end Random_Float_Template;
......
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