|
#############################################################################
##
## This file is part of GAP, a system for computational discrete algebra.
## This file's authors include Frank Lübeck, Max Neunhöffer.
##
## Copyright of GAP belongs to its developers, whose names are too numerous
## to list here. Please refer to the COPYRIGHT file for details.
##
## SPDX-License-Identifier: GPL-2.0-or-later
##
## This file implements the basic operations for some types of random
## sources.
##
###########################################################################
## Generic methods for random sources.
##
# The generic initializer of a random source creates a dummy object of the
# right type and then calls 'Init'.
InstallMethod(RandomSource, [IsOperation, IsObject], function(rep, seed)
local res;
res := Objectify(NewType(RandomSourcesFamily, rep), rec());
return Init(res, seed);
end);
InstallMethod(RandomSource, [IsOperation], function(rep)
return RandomSource(rep, 1);
end);
# A generic Reset if no seed is given (then seed=1 is chosen).
InstallMethod(Reset, [IsRandomSource], function(rs)
return Reset(rs, 1);
end);
# Generic fallback methods, such that it is sufficient to install Random for
# lists or for pairs of integers.
InstallMethod(Random, [IsRandomSource, IsInt, IsInt], function(rs, a, b)
local d, x, r, y;
d := b - a;
if d < 0 then
return fail;
elif a = b then
return a;
else
x := LogInt( d, 2 ) + 1;
r := 0;
while 0 < x do
y := Minimum( 10, x );
x := x - y;
r := r * 2 ^ y + Random( rs, [ 0 .. (2 ^ y - 1) ] );
od;
if d < r then
return Random( rs, a, b );
else
return a + r;
fi;
fi;
end);
InstallMethod(Random, [IsRandomSource, IsList and IsDenseList],
function(rs, list)
return list[Random(rs, 1, Length(list))];
end);
# A print method.
InstallMethod(PrintObj, [IsRandomSource], function(rs)
local cat;
cat := Difference(CategoriesOfObject(rs), [ "IsRandomSource" ]);
Print("<RandomSource in ", JoinStringsWithSeparator(cat, " and "), ">");
end);
############################################################################
## The classical GAP random generator as independent random sources.
##
# We need to compute modulo 2^28 repeatedly. If we just write 2^28 into the
# code, though, GAP needs to evaluate it to 268435456 each time it executes
# the function. To avoid this, we put it into a constant. As an additional
# trick, we actually store -2^28, which gives identical results, but has the
# benefit of being an immediate integer even on 32 bit systems.
BIND_CONSTANT("R_228", -2^28);
InstallMethod(Init, [IsGAPRandomSource, IsObject], function(rs, seed)
local R_N, R_X, i;
if seed = 1 then
rs!.R_N := 45;
rs!.R_X := [ 66318732, 86395905, 22233618, 21989103, 237245480,
264566285, 240037038, 264902875, 9274660, 180361945, 94688010, 24032135,
106293216, 27264613, 126456102, 243761907, 80312412, 2522186, 59575208,
70682510, 228947516, 173992210, 175178224, 250788150, 73030390, 210575942,
128491926, 194508966, 201311350, 63569414, 185485910, 62786150, 213986102,
88913350, 94904086, 252860454, 247700982, 233113990, 75685846, 196780518,
74570934, 7958751, 130274620, 247708693, 183364378, 82600777, 28385464,
184547675, 20423483, 75041763, 235736203, 54265107, 49075195, 100648387,
114539755 ];
elif IsInt(seed) then
R_N := 1;
R_X := [ seed mod R_228 ];
for i in [2..55] do
R_X[i] := (1664525 * R_X[i-1] + 1) mod R_228;
od;
for i in [1..99] do
R_N := R_N mod 55 + 1;
R_X[R_N] := (R_X[R_N] + R_X[(R_N+30) mod 55+1]) mod R_228;
od;
rs!.R_N := R_N;
rs!.R_X := R_X;
else
rs!.R_N := seed[1];
rs!.R_X := ShallowCopy(seed[2]);
fi;
return rs;
end);
InstallMethod(State, [IsGAPRandomSource], function(rs)
return [rs!.R_N, ShallowCopy(rs!.R_X)];
end);
InstallMethod(Reset, [IsGAPRandomSource, IsObject], function(rs, seed)
local old;
old := State(rs);
Init(rs, seed);
return old;
end);
InstallMethod(Random, [IsGAPRandomSource, IsList and IsDenseList],
function(rs, list)
local rx, rn;
if Length(list) < 2^28 then
rx := rs!.R_X;
rn := rs!.R_N mod 55 + 1;
rs!.R_N := rn;
rx[rn] := (rx[rn] + rx[(rn+30) mod 55+1]) mod R_228;
return list[ QUO_INT( rx[rn] * LEN_LIST(list), -R_228 ) + 1 ];
else
return list[Random(rs, 1, Length(list))];
fi;
end);
############################################################################
## We provide the "classical" GAP random generator via a random source.
##
if IsHPCGAP then
BIND_GLOBAL("RANDOM_SEED_COUNTER", FixedAtomicList(1, 0));
BIND_GLOBAL("GET_RANDOM_SEED_COUNTER", {} ->
ATOMIC_ADDITION(RANDOM_SEED_COUNTER, 1, 1) );
# HACK to enforce backwards compatibility: when reading this file,
# we are on the main thread, and want to initialize GlobalRandomSource
# with seed 1, without modifying the RANDOM_SEED_COUNTER, to get
# behavior identical with that of prior GAP versions, and also identical
# with regular GAP (on a single thread).
BIND_GLOBAL("GlobalRandomSource", RandomSource(IsGAPRandomSource, 1));
BindThreadLocalConstructor("GlobalRandomSource", {} ->
RandomSource(IsGAPRandomSource, GET_RANDOM_SEED_COUNTER()));
else
BindGlobal("GlobalRandomSource", RandomSource(IsGAPRandomSource, 1));
fi;
##############################################################################
## Random source using the Mersenne twister kernel functions.
##
InstallMethod(Init, [IsMersenneTwister, IsObject], function(rs, seed)
local st, endianseed, endiansys, perm, tmp, i;
if IsPlistRep(seed) and IsString(seed[1]) and Length(seed[1]) = 2504 then
st := ShallowCopy(seed[1]);
# maybe adjust endianness if seed comes from different machine:
endianseed := st{[2501..2504]};
endiansys := GlobalMersenneTwister!.state{[2501..2504]};
if endianseed <> endiansys then
perm := List(endiansys, c-> Position(endianseed, c));
tmp := "";
for i in [0..625] do
tmp{[4*i+1..4*i+4]} := st{4*i+perm};
od;
st := tmp;
fi;
rs!.state := st;
else
if not IsString(seed) then
seed := ShallowCopy(String(seed));
# padding such that length is positive and divisible by 4
while Length(seed) = 0 or Length(seed) mod 4 <> 0 do
Add(seed, CHAR_INT(0));
od;
fi;
rs!.state := InitRandomMT(seed);
fi;
return rs;
end);
InstallMethod(State, [IsMersenneTwister], function(rs)
return [ShallowCopy(rs!.state)];
end);
InstallMethod(Reset, [IsMersenneTwister, IsObject], function(rs, seed)
local old;
old := State(rs);
Init(rs, seed);
return old;
end);
# We do not need a 'Random' method for 'IsMersenneTwister' and a dense list,
# the default method for 'IsRandomSource' is enough.
InstallMethod(Random, [IsMersenneTwister, IsInt, IsInt], function(rs, a, b)
local d, nrbits, res;
d := b-a+1;
if d <= 0 then
return fail;
fi;
# Here we could return 'a' in the case 'd = 1'.
# However, this would change the sequence of random numbers
# w.r.t. earlier GAP versions.
nrbits := Log2Int(d) + 1;
repeat
res := RandomIntegerMT(rs!.state, nrbits);
until res < d;
return res + a;
end);
# One global Mersenne twister random source, can be used to overwrite
# the library Random(list) and Random(a,b) methods.
if IsHPCGAP then
BindThreadLocalConstructor("GlobalMersenneTwister", {} ->
RandomSource(IsMersenneTwister, String(GET_RANDOM_SEED_COUNTER())));
else
BindGlobal("GlobalMersenneTwister", RandomSource(IsMersenneTwister, "1"));
fi;
# default random method for lists and pairs of integers using the Mersenne
# twister
InstallMethod( Random, "for a dense internal list",
[ IsList and IsDenseList and IsInternalRep ], 100, function(l)
return l[Random(GlobalMersenneTwister, 1, Length(l))];
end );
InstallMethod( Random,
"for two integers",
IsIdenticalObj,
[ IsInt,
IsInt ],
0,
function(low, high)
return Random(GlobalMersenneTwister, low, high);
end );
(function()
local func;
func := function(installType)
return function(args...)
local filterpos, i, func, info;
# Check we understand arguments
# Second value must be an info string
if not IsString(args[2]) then
ErrorNoReturn("Second argument must be an info string");
fi;
# Info strings always tend to begin 'for ', and here we want
# to be able to edit it, so we check.
if args[2]{[1..23]} <> "for a random source and" then
ErrorNoReturn("Info string must begin 'for a random source and'");
fi;
# Filters must start with 'IsRandomSource'
for i in [1..Length(args)] do
if IsList(args[i]) and args[i][1] = IsRandomSource then
filterpos := i;
fi;
od;
if not IsBound(filterpos) then
ErrorNoReturn("Must use a list of filters beginning 'IsRandomSource'");
fi;
# Last argument must be the actual method
if not IsFunction(Last(args)) then
ErrorNoReturn("Argument list must end with the method");
fi;
# Install
CallFuncList(installType, args);
# Install random, wrapping random source argument
# Remove 'IsRandomSource' from the filter list
Remove(args[filterpos], 1);
# Correct info string by removing 'a random source and'
info := "for";
APPEND_LIST(info, args[2]{[24..Length(args[2])]});
args[2] := info;
func := Remove(args);
if Length(args[filterpos]) = 1 then
Add(args, x -> func(GlobalMersenneTwister,x));
elif Length(args[filterpos]) = 2 then
Add(args, {x,y} -> func(GlobalMersenneTwister,x,y));
else
Error("Only 2 or 3 argument methods supported");
fi;
CallFuncList(installType, args);
end;
end;
InstallGlobalFunction("InstallMethodWithRandomSource", func(InstallMethod));
InstallGlobalFunction("InstallOtherMethodWithRandomSource", func(InstallOtherMethod));
end)();
# This method must rank below Random(SomeRandomSource, IsList)
# for any random source SomeRandomSource, to avoid an infinite loop.
InstallMethodWithRandomSource( Random, "for a random source and a (finite) collection",
[ IsRandomSource, IsCollection and IsFinite ],
{} -> -RankFilter(IsCollection and IsFinite),
{rs, C} -> RandomList(rs, Enumerator( C ) ) );
[ Dauer der Verarbeitung: 0.31 Sekunden
(vorverarbeitet)
]
|
