########################################################################
##
## Copyright (C) 2020 The Octave Project Developers
##
## See the file COPYRIGHT.md in the top-level directory of this
## distribution or .
##
## This file is part of Octave.
##
## Octave is free software: you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING. If not, see
## .
##
########################################################################
## -*- texinfo -*-
## @deftypefn {} {} rng (@var{seed})
## @deftypefnx {} {} rng (@var{seed}, @var{generator})
## @deftypefnx {} {} rng ("shuffle")
## @deftypefnx {} {} rng ("shuffle", @var{generator})
## @deftypefnx {} {} rng ("default")
## @deftypefnx {} {@var{s} =} rng ()
## @deftypefnx {} {} rng (@var{s})
## @deftypefnx {} {@var{s} =} rng (@dots{})
## Set or query the seed of the random number generator used by rand and randn.
##
## Input scalar numeric value @qcode{seed} is used to initialize the state
## vector of the random number generator.
##
## The optional string @var{generator} specifies the type of random number
## generator to be used. Its value can be @qcode{"twister"}, @qcode{"v5uniform"}
## or @qcode{"v5normal"}. The @qcode{"twister"} keyword is described below.
## @qcode{"v5uniform"} and @qcode{"v5normal"} refer to older versions of Octave
## that used to use a different random number generator.
##
## The state or seed of the random number generator can be reset to a new random
## value using the @qcode{"shuffle"} keyword.
##
## The state or seed of the random number generator can be reset to its default
## values using the @qcode{"default"} keyword. The default values are to use the
## Mersenne Twister with a seed of 0.
##
## The optional return value @var{s} contains the state of the random number
## generator at the time the function is called (i.e. before it might be
## modified according to the input arguments). It is encoded as a structure
## variable with three fields: @qcode{"Type"}, @qcode{"Seed"} and
## @qcode{"State"}.
## Later, the random number generator can be restored to the state @var{s} using
## rng (s).
##
## By default, and with the @qcode{"twister"} option, pseudo-random sequences
## are computed using the Mersenne Twister with a period of @math{2^{19937}-1}
## (See @nospell{M. Matsumoto and T. Nishimura},
## @cite{Mersenne Twister: A 623-dimensionally equidistributed uniform
## pseudorandom number generator},
## @nospell{ACM} Trans.@: on Modeling and Computer Simulation Vol.@: 8, No.@: 1,
## pp.@: 3--30, January 1998,
## @url{http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html}).
## Do @strong{not} use for cryptography without securely hashing several
## returned values together, otherwise the generator state can be learned after
## reading 624 consecutive values.
##
## @seealso{rand, randn}
## @end deftypefn
function retval = rng (varargin)
if (nargin > 2 || nargout > 1)
print_usage ();
endif
## Store current settings of random number generator
## FIXME: there doesn't seem to be a way to query the type of generator
## currently used in Octave - assume "twister".
## FIXME: there doesn't seem to be a way to query the seed initialization
## value - use "Not applicable".
## FIXME: rand and randn use different generators - storing both states.
## For older Matlab generators (v4, v5), the settings are stored like this:
## struct ("Type","Legacy", "Seed", "Not applicable", "State",{[],[],...})
s = struct (...
"Type", "twister",... # generator name
"Seed", "Not applicable",... # seed initialization value
"State", {{rand("state"), randn("state")}}); # internal state of the generator
if (! nargin)
retval = s;
return;
endif
if (isscalar (varargin{1}) && isnumeric (varargin{1}) && isreal (varargin{1}) && varargin{1} >= 0)
s_rand = s_randn = varargin{1};
generator = check_generator (varargin(2:end));
elseif (ischar (varargin{1}) && strcmpi (varargin{1}, "shuffle"))
## Seed the random number generator based on the current time
s_rand = s_randn = "reset"; # or sum (1000*clock)
generator = check_generator (varargin(2:end));
elseif (ischar (varargin{1}) && strcmpi (varargin{1}, "default") && nargin == 1)
generator = "twister";
s_rand = s_randn = 0; # In Matlab, seed 0 corresponds to 5489
elseif (isstruct (varargin{1}) && isscalar (varargin{1}) && nargin == 1)
if (numfields (varargin{1}) != 3 || ! isfield (varargin{1}, "Type")
|| ! isfield (varargin{1}, "Seed") || ! isfield (varargin{1}, "State"))
error ("Input structure not compatible with the one returned by rng ()");
endif
## Only the internal state "State" and generator type "Type" are needed
generator = varargin{1}.Type;
if (iscell (varargin{1}.State))
[s_rand, s_randn] = deal (varargin{1}.State{:});
else
s_rand = s_randn = varargin{1}.State;
endif
else
print_usage ();
endif
## Set the type of random number generator and seed it
if (isempty (generator))
generator = s.Type;
endif
switch generator
case "twister"
rand ("state", s_rand);
randn ("state", s_randn);
case "legacy"
rand ("seed", s_rand);
randn ("seed", s_randn);
case "v5uniform"
rand ("seed", s_rand);
case "v5normal"
randn ("seed", s_randn);
otherwise
error ("Unknown type of random number generator");
endswitch
if (nargout > 0)
retval = s;
endif
endfunction
function gen = check_generator (val)
if (isempty (val))
gen = "";
return;
elseif (! iscellstr (val))
error ("Second input must be a type of random number generator");
endif
gen = tolower (char (val));
if (ismember (gen, {"simdtwister", "combrecursive", "philox", "threefry", "multfibonacci", "v4"}))
error ("This random number generator is not available in Octave");
elseif (! ismember (gen, {"twister", "v5uniform", "v5normal"}))
error ("This type of random number generator is unknown");
endif
endfunction
%!test
%! rng (42);
%! ru1 = rand ();
%! rn1 = randn ();
%! rng (42);
%! ru2 = rand ();
%! rn2 = randn ();
%! assert (ru2, ru1);
%! assert (rn2, rn1);
%! s1 = rng ();
%! s2 = rng (42);
%! assert (isequal (s1, s2));
%! ru1 = rand ();
%! rn1 = randn ();
%! s3 = rng (42);
%! ru2 = rand ();
%! rn2 = randn ();
%! assert (ru2, ru1);
%! assert (rn2, rn1);
%!test
%! rng (42, "twister");
%! ru1 = rand ();
%! rn1 = randn ();
%! rng (42, "twister");
%! ru2 = rand ();
%! rn2 = randn ();
%! assert (ru2, ru1);
%! assert (rn2, rn1);
%! s1 = rng ();
%! s2 = rng (42, "twister");
%! assert (isequal (s1, s2));
%! ru1 = rand ();
%! rn1 = randn ();
%! s3 = rng (42, "twister");
%! ru2 = rand ();
%! rn2 = randn ();
%! assert (ru2, ru1);
%! assert (rn2, rn1);
%!test
%! rng ("shuffle");
%! rng ("shuffle", "twister");
%! s = rng ("shuffle");
%! assert (! isequal (s, rng ("shuffle")));
%! s = rng ("shuffle", "twister");
%! assert (! isequal (s, rng ("shuffle", "twister")));
%!test
%! rng ("default");
%! ru1 = rand ();
%! rn1 = randn ();
%! rng (0, "twister");
%! ru2 = rand ();
%! rn2 = randn ();
%! assert (ru2, ru1);
%! assert (rn2, rn1);
%! rng (0, "twister");
%! s = rng ("default");
%! assert (isequal (s, rng ()));
%!test
%! s = rng ();
%! ru1 = rand ();
%! rn1 = randn ();
%! rng (s);
%! ru2 = rand ();
%! rn2 = randn ();
%! assert (ru2, ru1);
%! assert (rn2, rn1);
%! rng (42); rand (1,2); x = rand (1,2);
%! rng (42); rand (1,2); s = rng (); y = rand (1,2);
%! assert (x, y);
%! rng (s); z = rand (1,2);
%! assert (x, z);
%! s1 = rng ();
%! s2 = rng (s1);
%! assert (isequal (s1, s2));
## Test input validation
%!error rng (1, 2, 3)
%!error [a, b] = rng ()
%!error rng ({})
%!error rng ("unknown")
%!error rng (eye (2))
%!error rng (struct ("Seed", {1 ,2}))
%!error rng (-1)
%!error rng (struct ("Type",[],"State",[],"Seed",[]), 2)
%!error rng ("default", "twister")
%!error rng (0, struct ())
%!error rng (0, "unknown")
%!error rng ("shuffle", struct ())
%!error rng ("shuffle", "unknown")