[m-rev.] RNGs
Mark Brown
mark at mercurylang.org
Tue Aug 13 17:32:33 AEST 2019
Hi all,
I've just made a pull request to add some RNGs to the standard library, as
well as a typeclass interface to them. Pull request is here:
https://github.com/Mercury-Language/mercury/pull/70
The diff is also attached, for anyone who would prefer to comment on
it via email.
Mark
-------------- next part --------------
commit f102ddaf4c128137b27ac6e85ac881faaf435c9b
Author: Mark Brown <mark at mercurylang.org>
Date: Tue Aug 13 16:42:38 2019 +1000
Some RNGs for the standard library, plus typeclass interfaces.
library/rng.m:
Top-level module containing typeclasses and generic routines,
as well as overall documentation.
library/rng.binfile.m:
A generator that reads from a file.
library/rng.marsaglia.m:
Fast, simple generator. Diehard results for default seed:
PASSED: 109
WEAK: 2
FAILED: 3
library/rng.tausworthe.m:
Combined Tausworthe generators. Diehard results for default seed:
generator T3
PASSED: 113
WEAK: 1
FAILED: 0
generator T4
PASSED: 112
WEAK: 2
FAILED: 0
library/library.m:
library/MODULES_DOC:
Add the new modules to the public interface.
library/float.m:
Add functions to convert from (u)int32/64 to float.
library/uint32.m:
Add function to convert from uint64 to uint32.
tests/hard_coded/Mmakefile:
tests/hard_coded/rng1.{m,exp}:
tests/hard_coded/rng2.{m,exp}:
Test the generators. The expected outputs are derived from
the output of the C reference implementations.
diff --git a/library/MODULES_DOC b/library/MODULES_DOC
index 08f1681..822524e 100644
--- a/library/MODULES_DOC
+++ b/library/MODULES_DOC
@@ -61,6 +61,10 @@ ranges.m
rational.m
rbtree.m
require.m
+rng.m
+rng.binfile.m
+rng.marsaglia.m
+rng.tausworthe.m
rtree.m
set.m
set_bbbtree.m
diff --git a/library/float.m b/library/float.m
index e19981c..4e2abfe 100644
--- a/library/float.m
+++ b/library/float.m
@@ -133,13 +133,25 @@
:- func from_int8(int8) = float.
% Convert a signed 16-bit integer into a float.
- % Always succeeds as all unsigned 8-bit integers have an exact
+ % Always succeeds as all signed 16-bit integers have an exact
% floating-point representation.
%
:- func from_int16(int16) = float.
+ % Convert a signed 32-bit integer into a float.
+ % Always succeeds as all signed 32-bit integers have an exact
+ % floating-point representation.
+ %
+:- func from_int32(int32) = float.
+
+ % Convert a signed 64-bit integer into a float.
+ % The behaviour when the integer exceeds the range of what can be
+ % exactly represented by a float is undefined.
+ %
+:- func from_int64(int64) = float.
+
% Convert an unsigned 8-bit integer into a float.
- % Always succeeds as all signed 16-bit integers have an exact
+ % Always succeeds as all unsigned 8-bit integers have an exact
% floating-point representation.
%
:- func from_uint8(uint8) = float.
@@ -150,6 +162,18 @@
%
:- func from_uint16(uint16) = float.
+ % Convert an unsigned 32-bit integer into a float.
+ % Always succeeds as all unsigned 32-bit integers have an exact
+ % floating-point representation.
+ %
+:- func from_uint32(uint32) = float.
+
+ % Convert an unsigned 64-bit integer into a float.
+ % The behaviour when the integer exceeds the range of what can be
+ % exactly represented by a float is undefined.
+ %
+:- func from_uint64(uint64) = float.
+
%---------------------------------------------------------------------------%
%
% Conversion to int.
@@ -558,6 +582,130 @@ X / Y = Z :-
").
%---------------------------------------------------------------------------%
+
+:- pragma foreign_proc("C",
+ from_int32(Int32Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe, will_not_modify_trail,
+ does_not_affect_liveness],
+"
+ FloatVal = Int32Val;
+").
+
+:- pragma foreign_proc("C#",
+ from_int32(Int32Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ FloatVal = (double) Int32Val;
+").
+
+:- pragma foreign_proc("Java",
+ from_int32(Int32Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ FloatVal = (double) Int32Val;
+").
+
+:- pragma foreign_proc("Erlang",
+ from_int32(Int32Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ FloatVal = float(Int32Val)
+").
+
+%---------------------------------------------------------------------------%
+
+:- pragma foreign_proc("C",
+ from_uint32(UInt32Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe, will_not_modify_trail,
+ does_not_affect_liveness],
+"
+ FloatVal = UInt32Val;
+").
+
+:- pragma foreign_proc("C#",
+ from_uint32(UInt32Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ FloatVal = (double) UInt32Val;
+").
+
+:- pragma foreign_proc("Java",
+ from_uint32(UInt32Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ FloatVal = (double) (UInt32Val & 0xffffffff);
+").
+
+:- pragma foreign_proc("Erlang",
+ from_uint32(UInt32Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ FloatVal = float(UInt32Val)
+").
+
+%---------------------------------------------------------------------------%
+
+:- pragma foreign_proc("C",
+ from_int64(Int64Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe, will_not_modify_trail,
+ does_not_affect_liveness],
+"
+ FloatVal = Int64Val;
+").
+
+:- pragma foreign_proc("C#",
+ from_int64(Int64Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ FloatVal = (double) Int64Val;
+").
+
+:- pragma foreign_proc("Java",
+ from_int64(Int64Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ FloatVal = (double) Int64Val;
+").
+
+:- pragma foreign_proc("Erlang",
+ from_int64(Int64Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ FloatVal = float(Int64Val)
+").
+
+%---------------------------------------------------------------------------%
+
+:- pragma foreign_proc("C",
+ from_uint64(UInt64Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe, will_not_modify_trail,
+ does_not_affect_liveness],
+"
+ FloatVal = UInt64Val;
+").
+
+:- pragma foreign_proc("C#",
+ from_uint64(UInt64Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ FloatVal = (double) UInt64Val;
+").
+
+:- pragma foreign_proc("Java",
+ from_uint64(UInt64Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ FloatVal = (double) UInt64Val;
+").
+
+:- pragma foreign_proc("Erlang",
+ from_uint64(UInt64Val::in) = (FloatVal::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ FloatVal = float(UInt64Val)
+").
+
+%---------------------------------------------------------------------------%
%
% Conversion to ints.
%
diff --git a/library/library.m b/library/library.m
index 4330097..9b79605 100644
--- a/library/library.m
+++ b/library/library.m
@@ -124,6 +124,10 @@
:- import_module rational.
:- import_module rbtree.
:- import_module require.
+:- import_module rng.
+:- import_module rng.binfile.
+:- import_module rng.marsaglia.
+:- import_module rng.tausworthe.
:- import_module robdd.
:- import_module rtree.
:- import_module set.
@@ -309,6 +313,10 @@ mercury_std_library_module("rational").
mercury_std_library_module("rbtree").
mercury_std_library_module("region_builtin").
mercury_std_library_module("require").
+mercury_std_library_module("rng").
+mercury_std_library_module("rng.binfile").
+mercury_std_library_module("rng.marsaglia").
+mercury_std_library_module("rng.tausworthe").
mercury_std_library_module("robdd").
mercury_std_library_module("rtree").
mercury_std_library_module("rtti_implementation").
diff --git a/library/rng.binfile.m b/library/rng.binfile.m
new file mode 100644
index 0000000..36c1966
--- /dev/null
+++ b/library/rng.binfile.m
@@ -0,0 +1,96 @@
+%---------------------------------------------------------------------------%
+% vim: ft=mercury ts=4 sts=4 sw=4 et
+%---------------------------------------------------------------------------%
+% Copyright (C) 2019 The Mercury team.
+% This file is distributed under the terms specified in COPYING.LIB.
+%---------------------------------------------------------------------------%
+%
+% File: rng.binfile.m
+% Main author: Mark Brown
+%
+% "Random" number generator that reads numbers from a binary file.
+%
+%---------------------------------------------------------------------------%
+%---------------------------------------------------------------------------%
+
+:- module rng.binfile.
+:- interface.
+
+:- import_module io.
+
+%---------------------------------------------------------------------------%
+
+:- type binfile.
+:- instance urng(binfile, io).
+
+ % Open a binfile generator from a filename. This should be closed
+ % when no longer needed.
+ %
+:- pred open(string, io.res(binfile), io, io).
+:- mode open(in, out, di, uo) is det.
+
+ % Close a binfile generator.
+ %
+:- pred close(binfile, io, io).
+:- mode close(in, di, uo) is det.
+
+%---------------------------------------------------------------------------%
+
+ % Generate a number between 0 and max_uint64. This reads 8 bytes
+ % at a time from the binfile and interprets them as an unsigned,
+ % big-endian integer.
+ %
+ % Throws an exception if the end-of-file is reached.
+ %
+:- pred rand(binfile, uint64, io, io).
+:- mode rand(in, out, di, uo) is det.
+
+%---------------------------------------------------------------------------%
+
+:- implementation.
+
+:- import_module require.
+:- import_module uint64.
+
+%---------------------------------------------------------------------------%
+
+:- type binfile
+ ---> binfile(binary_input_stream).
+
+:- instance urng(binfile, io) where [
+ pred(urandom/4) is rand,
+ ( urandom_max(_) = uint64.max_uint64 )
+].
+
+%---------------------------------------------------------------------------%
+
+open(Filename, Res, !IO) :-
+ io.open_binary_input(Filename, Res0, !IO),
+ (
+ Res0 = ok(Stream),
+ Res = ok(binfile(Stream))
+ ;
+ Res0 = error(E),
+ Res = error(E)
+ ).
+
+close(binfile(Stream), !IO) :-
+ io.close_binary_input(Stream, !IO).
+
+%---------------------------------------------------------------------------%
+
+rand(binfile(Stream), N, !IO) :-
+ io.read_binary_uint64_be(Stream, Res, !IO),
+ (
+ Res = ok(N)
+ ;
+ ( Res = eof
+ ; Res = incomplete(_)
+ ),
+ unexpected($pred, "end of file")
+ ;
+ Res = error(E),
+ unexpected($pred, io.error_message(E))
+ ).
+
+%---------------------------------------------------------------------------%
diff --git a/library/rng.m b/library/rng.m
new file mode 100644
index 0000000..24e8086
--- /dev/null
+++ b/library/rng.m
@@ -0,0 +1,345 @@
+%---------------------------------------------------------------------------%
+% vim: ft=mercury ts=4 sts=4 sw=4 et
+%---------------------------------------------------------------------------%
+% Copyright (C) 2019 The Mercury team.
+% This file is distributed under the terms specified in COPYING.LIB.
+%---------------------------------------------------------------------------%
+%
+% File: rng.m
+% Main author: Mark Brown
+%
+% This module provides an interface to several random number generators,
+% which can be found in the submodules.
+%
+% Two styles of the interface are provided, a ground style and a
+% unique style. Each has its own advantages and disadvantages:
+%
+% - Ground RNGs are easier to use; for example they can be easily
+% stored in larger data structures.
+% - Ground RNGs are easier to implement instances for.
+% - Unique RNGs are able to use destructive update, and therefore
+% are often able to operate more efficiently.
+% - Unique RNGs need to be explicitly duplicated (i.e., to produce
+% a new generator that will generate the same sequence of numbers).
+% This may be regarded as an advantage or a disadvantage.
+% - Some RNGs, for example the binfile generator that reads data from
+% a file, use the IO state and therefore must use the unique interface.
+%
+% Each RNG defined in the submodules is natively one of these two styles.
+% Conversion between the two styles can be done with make_urng/3 and
+% make_shared_rng/2, below, although this incurs additional overhead.
+%
+%
+% Example, ground style:
+%
+% main(!IO) :-
+% RNG0 = rng.marsaglia.init,
+% roll(RNG0, RNG1, !IO),
+% roll(RNG1, _, !IO).
+%
+% :- pred roll(RNG, RNG, io, io) <= rng(RNG).
+% :- mode roll(in, out, di, uo) is det.
+%
+% roll(!RNG, !IO) :-
+% random_int(1, 6, N, !RNG),
+% io.format("You rolled a %d\n", [i(N)], !IO).
+%
+%
+% Example, unique style:
+%
+% main(!IO) :-
+% rng.tausworthe.init_t3(RP, RS0),
+% roll(RP, RS0, RS1, !IO),
+% roll(RP, RS1, _, !IO).
+%
+% :- pred roll(RP, RS, RS, io, io) <= urng(RP, RS).
+% :- mode roll(in, di, uo, di, uo) is det. % note unique modes
+%
+% roll(RP, !RS, !IO) :-
+% urandom_int(RP, 1, 6, N, !RS),
+% io.format("You rolled a %d\n", [i(N)], !IO).
+%
+%
+% Example, converting style:
+%
+% main(!IO) :-
+% rng.tausworthe.init_t3(RP, RS),
+% RNG0 = make_shared_rng(RP, RS),
+% random_int(1, 6, N, RNG0, RNG1),
+% ...
+%
+% main(!IO) :-
+% RNG = rng.marsaglia.init,
+% make_urng(RNG, RP, RS0),
+% urandom_int(RP, 1, 6, N, RS0, RS1),
+% ...
+%
+%---------------------------------------------------------------------------%
+%---------------------------------------------------------------------------%
+
+:- module rng.
+:- interface.
+
+:- include_module binfile.
+:- include_module marsaglia.
+:- include_module tausworthe.
+
+%---------------------------------------------------------------------------%
+
+ % random_int(Start, Range, N, !RNG)
+ %
+ % Generate a random integer between Start and Start+Range-1 inclusive.
+ % Throws an exception if Range < 1 or Range > random_max.
+ %
+:- pred random_int(int, int, int, RNG, RNG) <= rng(RNG).
+:- mode random_int(in, in, out, in, out) is det.
+
+ % Generate a random float between 0.0 and 1.0, inclusive.
+ %
+:- pred random_float(float, RNG, RNG) <= rng(RNG).
+:- mode random_float(out, in, out) is det.
+
+ % Generate two random floats from a normal distribution with
+ % mean 0 and standard deviation 1, using the Box-Muller method.
+ %
+ % We generate two at a time for efficiency; they are independent of
+ % each other.
+ %
+:- pred random_gauss(float, float, RNG, RNG) <= rng(RNG).
+:- mode random_gauss(out, out, in, out) is det.
+
+%---------------------------------------------------------------------------%
+
+ % Interface to random number generators.
+ %
+:- typeclass rng(RNG) where [
+
+ % Generate a random integer between 0 and random_max, inclusive.
+ %
+ pred random(uint64, RNG, RNG),
+ mode random(out, in, out) is det,
+
+ % Return the largest integer that can be generated.
+ %
+ func random_max(RNG) = uint64
+].
+
+%---------------------------------------------------------------------------%
+%---------------------------------------------------------------------------%
+
+ % urandom_int(RP, Start, Range, N, !RS)
+ %
+ % Generate a random integer between Start and Start+Range-1 inclusive.
+ % Throws an exception if Range < 1 or Range > urandom_max.
+ %
+:- pred urandom_int(RP, int, int, int, RS, RS) <= urng(RP, RS).
+:- mode urandom_int(in, in, in, out, di, uo) is det.
+
+ % Generate a random float between 0.0 and 1.0, inclusive.
+ %
+:- pred urandom_float(RP, float, RS, RS) <= urng(RP, RS).
+:- mode urandom_float(in, out, di, uo) is det.
+
+ % Generate two random floats from a normal distribution with
+ % mean 0 and standard deviation 1, using the Box-Muller method.
+ %
+ % We generate two at a time for efficiency; they are independent of
+ % each other.
+ %
+:- pred urandom_gauss(RP, float, float, RS, RS) <= urng(RP, RS).
+:- mode urandom_gauss(in, out, out, di, uo) is det.
+
+%---------------------------------------------------------------------------%
+
+ % Interface to unique random number generators. Callers need to
+ % ensure they preserve the uniqueness of the random state, and in
+ % turn instances can use destructive update on it.
+ %
+:- typeclass urng(RP, RS) <= (RP -> RS) where [
+
+ % Generate a random integer between 0 and random_max, inclusive.
+ %
+ pred urandom(RP, uint64, RS, RS),
+ mode urandom(in, out, di, uo) is det,
+
+ % Return the largest integer that can be generated.
+ %
+ func urandom_max(RP) = uint64
+].
+
+:- typeclass urng_dup(RS) where [
+
+ % urandom_dup(!RS, !:RSdup)
+ %
+ % Create a duplicate random state that will generate the
+ % same sequence of integers.
+ %
+ pred urandom_dup(RS, RS, RS),
+ mode urandom_dup(di, uo, uo) is det
+].
+
+%---------------------------------------------------------------------------%
+%---------------------------------------------------------------------------%
+
+ % Convert any rng into a urng. This creates some additional overhead
+ % in the form of additional typeclass method calls.
+ %
+:- type urng_params(RNG).
+:- type urng_state(RNG).
+
+:- instance urng(urng_params(RNG), urng_state(RNG)) <= rng(RNG).
+:- instance urng_dup(urng_state(RNG)) <= rng(RNG).
+
+:- pred make_urng(RNG, urng_params(RNG), urng_state(RNG)) <= rng(RNG).
+:- mode make_urng(in, out, uo) is det.
+
+%---------------------------------------------------------------------------%
+
+ % Convert any urng into an rng. This duplicates the state every time
+ % a random number is generated, hence may use significantly more
+ % memory than if the unique version is used directly.
+ %
+:- type shared_rng(RP, RS).
+
+:- instance rng(shared_rng(RP, RS)) <= (urng(RP, RS), urng_dup(RS)).
+
+:- func make_shared_rng(RP, RS) = shared_rng(RP, RS).
+:- mode make_shared_rng(in, di) = out is det.
+
+%---------------------------------------------------------------------------%
+%---------------------------------------------------------------------------%
+
+:- implementation.
+
+:- import_module int.
+:- import_module float.
+:- import_module math.
+:- import_module uint64.
+
+%---------------------------------------------------------------------------%
+
+random_int(Start, Range0, N, !RNG) :-
+ Range = uint64.det_from_int(Range0),
+ random(N0, !RNG),
+ Max = random_max(!.RNG),
+ N1 = N0 // (Max // Range),
+ ( if N1 < Range then
+ N = Start + uint64.cast_to_int(N1)
+ else
+ random_int(Start, Range0, N, !RNG)
+ ).
+
+random_float(F, !RNG) :-
+ random(N, !RNG),
+ Max = random_max(!.RNG),
+ F = float.from_uint64(N) / float.from_uint64(Max).
+
+random_gauss(U, V, !RNG) :-
+ random_float(X, !RNG),
+ random_float(Y, !RNG),
+ ( if gauss(X, Y, U0, V0) then
+ U = U0,
+ V = V0
+ else
+ random_gauss(U, V, !RNG)
+ ).
+
+%---------------------------------------------------------------------------%
+
+urandom_int(RP, Start, Range0, N, !RS) :-
+ Range = uint64.det_from_int(Range0),
+ urandom(RP, N0, !RS),
+ Max = urandom_max(RP),
+ N1 = N0 // (Max // Range),
+ ( if N1 < Range then
+ N = Start + uint64.cast_to_int(N1)
+ else
+ urandom_int(RP, Start, Range0, N, !RS)
+ ).
+
+urandom_float(RP, F, !RS) :-
+ urandom(RP, N, !RS),
+ Max = urandom_max(RP),
+ F = float.from_uint64(N) / float.from_uint64(Max).
+
+urandom_gauss(RP, U, V, !RS) :-
+ urandom_float(RP, X, !RS),
+ urandom_float(RP, Y, !RS),
+ ( if gauss(X, Y, U0, V0) then
+ U = U0,
+ V = V0
+ else
+ urandom_gauss(RP, U, V, !RS)
+ ).
+
+%---------------------------------------------------------------------------%
+
+:- pred gauss(float, float, float, float).
+:- mode gauss(in, in, out, out) is semidet.
+
+gauss(X0, Y0, U, V) :-
+ X = 2.0 * X0 - 1.0,
+ Y = 2.0 * Y0 - 1.0,
+ S = X * X + Y * Y,
+ S > 0.0,
+ S < 1.0,
+ Fac = math.sqrt(-2.0 * math.ln(S) / S),
+ U = X * Fac,
+ V = Y * Fac.
+
+%---------------------------------------------------------------------------%
+%---------------------------------------------------------------------------%
+
+:- type urng_params(RNG)
+ ---> urng_params(
+ urng_max :: uint64
+ ).
+
+:- type urng_state(RNG)
+ ---> urng_state(
+ urng_rng :: RNG
+ ).
+
+:- instance urng(urng_params(RNG), urng_state(RNG)) <= rng(RNG) where [
+ ( urandom(_, N, RS0, RS) :-
+ RS0 = urng_state(RNG0),
+ random(N, RNG0, RNG),
+ RS = unsafe_promise_unique(urng_state(RNG))
+ ),
+ ( urandom_max(RP) = RP ^ urng_max )
+].
+
+:- instance urng_dup(urng_state(RNG)) <= rng(RNG) where [
+ ( urandom_dup(RS, RS1, RS2) :-
+ RS1 = unsafe_promise_unique(RS),
+ RS2 = unsafe_promise_unique(RS)
+ )
+].
+
+make_urng(RNG, RP, RS) :-
+ RP = urng_params(random_max(RNG)),
+ RS = unsafe_promise_unique(urng_state(RNG)).
+
+%---------------------------------------------------------------------------%
+
+:- type shared_rng(RP, RS)
+ ---> shared_rng(
+ shared_rng_params :: RP,
+ shared_rng_state :: RS
+ ).
+
+:- instance rng(shared_rng(RP, RS)) <= (urng(RP, RS), urng_dup(RS)) where [
+ ( random(N, RNG0, RNG) :-
+ RNG0 = shared_rng(RP, RS0),
+ RS1 = unsafe_promise_unique(RS0),
+ urandom_dup(RS1, _, RS2),
+ urandom(RP, N, RS2, RS),
+ RNG = shared_rng(RP, RS)
+ ),
+ ( random_max(RNG) = urandom_max(RNG ^ shared_rng_params) )
+].
+
+make_shared_rng(RP, RS) = shared_rng(RP, RS).
+
+%---------------------------------------------------------------------------%
+%---------------------------------------------------------------------------%
diff --git a/library/rng.marsaglia.m b/library/rng.marsaglia.m
new file mode 100644
index 0000000..8afac69
--- /dev/null
+++ b/library/rng.marsaglia.m
@@ -0,0 +1,100 @@
+%---------------------------------------------------------------------------%
+% vim: ft=mercury ts=4 sts=4 sw=4 et
+%---------------------------------------------------------------------------%
+% Copyright (C) 2019 The Mercury team.
+% This file is distributed under the terms specified in COPYING.LIB.
+%---------------------------------------------------------------------------%
+%
+% File: rng.marsaglia.m
+% Main author: Mark Brown
+%
+% Very fast concatenation of two 16-bit MWC generators.
+%
+% http://gcrhoads.byethost4.com/Code/Random/marsaglia.c
+%
+% "Algorithm recommended by Marsaglia."
+%
+%---------------------------------------------------------------------------%
+%---------------------------------------------------------------------------%
+
+:- module rng.marsaglia.
+:- interface.
+
+%---------------------------------------------------------------------------%
+
+:- type marsaglia.
+
+:- instance rng(marsaglia).
+
+ % Initialise a marsaglia RNG with the default seed.
+ %
+:- func init = marsaglia.
+
+ % Initialise a marsaglia RNG with the given seed.
+ %
+:- func seed(uint32, uint32) = marsaglia.
+
+%---------------------------------------------------------------------------%
+
+ % Generate a random number between 0 and max_uint32.
+ %
+:- pred rand(uint32, marsaglia, marsaglia).
+:- mode rand(out, in, out) is det.
+
+%---------------------------------------------------------------------------%
+
+:- implementation.
+
+:- import_module uint32.
+
+%---------------------------------------------------------------------------%
+
+:- type marsaglia
+ ---> marsaglia(uint64).
+
+:- instance rng(marsaglia) where [
+ ( random(N, !RNG) :-
+ rand(N0, !RNG),
+ N = uint32.cast_to_uint64(N0)
+ ),
+ ( random_max(_) = uint32.cast_to_uint64(uint32.max_uint32) )
+].
+
+%---------------------------------------------------------------------------%
+
+init = seed(0u32, 0u32).
+
+seed(SX0, SY0) = RNG :-
+ SX = ( if SX0 = 0u32 then 521288629u32 else SX0 ),
+ SY = ( if SY0 = 0u32 then 362436069u32 else SY0 ),
+ RNG = marsaglia(pack_uint64(SX, SY)).
+
+%---------------------------------------------------------------------------%
+
+rand(N, RNG0, RNG) :-
+ RNG0 = marsaglia(S0),
+ unpack_uint64(S0, SX0, SY0),
+ A = 18000u32,
+ B = 30903u32,
+ M = 0xffffu32,
+ SX = A * (SX0 /\ M) + (SX0 >> 16),
+ SY = B * (SY0 /\ M) + (SY0 >> 16),
+ N = (SX << 16) + (SY /\ M),
+ S = pack_uint64(SX, SY),
+ RNG = marsaglia(S).
+
+%---------------------------------------------------------------------------%
+
+:- func pack_uint64(uint32, uint32) = uint64.
+
+pack_uint64(Hi, Lo) =
+ (uint32.cast_to_uint64(Hi) << 32) + uint32.cast_to_uint64(Lo).
+
+:- pred unpack_uint64(uint64, uint32, uint32).
+:- mode unpack_uint64(in, out, out) is det.
+
+unpack_uint64(S, Hi, Lo) :-
+ Hi = uint32.cast_from_uint64(S >> 32),
+ Lo = uint32.cast_from_uint64(S /\ 0xffffffffu64).
+
+%---------------------------------------------------------------------------%
diff --git a/library/rng.tausworthe.m b/library/rng.tausworthe.m
new file mode 100644
index 0000000..38aa059
--- /dev/null
+++ b/library/rng.tausworthe.m
@@ -0,0 +1,302 @@
+%---------------------------------------------------------------------------%
+% vim: ft=mercury ts=4 sts=4 sw=4 et
+%---------------------------------------------------------------------------%
+% Copyright (C) 2019 The Mercury team.
+% This file is distributed under the terms specified in COPYING.LIB.
+%---------------------------------------------------------------------------%
+%
+% File: rng.tausworthe.m
+% Main author: Mark Brown
+%
+% Combined Tausworthe-type generators. See:
+%
+% Pierre L'Ecuyer, "Maximally Equidistributed Combined Tausworthe Generators",
+% Mathematics of Computation, vol. 65, no. 213 (1996)
+% Pierre L'Ecuyer, "Tables of Maximally-Equidistributed Combined LFSR
+% Generators", Mathematics of Computation, vol. 68, no. 225 (1999)
+%
+% http://gcrhoads.byethost4.com/Code/Random/tausworth.c
+% http://gcrhoads.byethost4.com/Code/Random/tausworth4.c
+%
+%---------------------------------------------------------------------------%
+%---------------------------------------------------------------------------%
+
+:- module rng.tausworthe.
+:- interface.
+
+:- import_module maybe.
+
+%---------------------------------------------------------------------------%
+
+:- type params.
+:- type state.
+
+:- instance urng(params, state).
+:- instance urng_dup(state).
+
+ % Initialise a 3-combo tausworthe RNG with the default seed
+ % and parameters.
+ %
+:- pred init_t3(params, state).
+:- mode init_t3(out, uo) is det.
+
+ % Initialise a 4-combo tausworthe RNG with the default seed
+ % and parameters.
+ %
+:- pred init_t4(params, state).
+:- mode init_t4(out, uo) is det.
+
+ % Initialise a 3-combo tausworthe RNG with the given seed.
+ % If given, the first argument selects from one of two sets of
+ % parameters, depending on its value modulo 2.
+ %
+:- pred seed_t3(maybe(int), uint32, uint32, uint32, params, state).
+:- mode seed_t3(in, in, in, in, out, uo) is det.
+
+ % Initialise a 4-combo tausworthe RNG with the given seed.
+ % If given, the first argument selects from one of 62 sets of
+ % parameters, depending on its value modulo 62.
+ %
+:- pred seed_t4(maybe(int), uint32, uint32, uint32, uint32, params, state).
+:- mode seed_t4(in, in, in, in, in, out, uo) is det.
+
+%---------------------------------------------------------------------------%
+
+ % Generate a random number between 0 and max_uint32. Throws an
+ % exception if the params and state are not the same size (i.e.,
+ % both 3-combo or both 4-combo).
+ %
+:- pred rand(params, uint32, state, state).
+:- mode rand(in, out, di, uo) is det.
+
+ % Duplicate a tausworthe RNG state.
+ %
+:- pred dup(state, state, state).
+:- mode dup(di, uo, uo) is det.
+
+%---------------------------------------------------------------------------%
+
+:- implementation.
+
+:- import_module array.
+:- import_module int.
+:- import_module list.
+:- import_module require.
+:- import_module uint32.
+
+%---------------------------------------------------------------------------%
+
+:- type params
+ ---> params(
+ qs :: array(int),
+ ps :: array(int),
+ shft :: array(int),
+ mask :: array(uint32)
+ ).
+
+:- type state
+ ---> state(
+ seed :: array(uint32)
+ ).
+
+:- instance urng(params, state) where [
+ ( urandom(RP, N, !RS) :-
+ rand(RP, N0, !RS),
+ N = uint32.cast_to_uint64(N0)
+ ),
+ ( urandom_max(_) = uint32.cast_to_uint64(uint32.max_uint32) )
+].
+
+:- instance urng_dup(state) where [
+ pred(urandom_dup/3) is dup
+].
+
+dup(S, S1, S2) :-
+ S = state(A),
+ S1 = unsafe_promise_unique(S),
+ S2 = unsafe_promise_unique(state(array.copy(A))).
+
+%---------------------------------------------------------------------------%
+
+:- pred seed(array(int), array(int), array(uint32), params, state).
+:- mode seed(in, in, array_di, out, uo) is det.
+
+seed(Qs, Ps, Seed0, RP, RS) :-
+ Size = array.size(Seed0),
+ Ks = array([31, 29, 28, 25]),
+ Ds = array([390451501u32, 613566701u32, 858993401u32, 943651322u32]),
+ Shft0 = array.init(Size, 0),
+ Mask0 = array.init(Size, 0u32),
+ seed_2(0, Size, Ks, Ps, Ds, Shft0, Shft, Mask0, Mask, Seed0, Seed),
+ RP = params(Qs, Ps, Shft, Mask),
+ RS0 = unsafe_promise_unique(state(Seed)),
+ rand(RP, _, RS0, RS).
+
+:- pred seed_2(int, int, array(int), array(int), array(uint32),
+ array(int), array(int), array(uint32), array(uint32),
+ array(uint32), array(uint32)).
+:- mode seed_2(in, in, in, in, in,
+ array_di, array_uo, array_di, array_uo, array_di, array_uo) is det.
+
+seed_2(I, Size, Ks, Ps, Ds, !Shft, !Mask, !Seed) :-
+ ( if I < Size then
+ K = array.lookup(Ks, I),
+ P = array.lookup(Ps, I),
+ S = array.lookup(!.Seed, I),
+ J = 32 - K,
+ array.set(I, K - P, !Shft),
+ array.set(I, uint32.max_uint32 << J, !Mask),
+ ( if S > (1u32 << J) then
+ true
+ else
+ D = array.lookup(Ds, I),
+ array.set(I, D, !Seed)
+ ),
+ seed_2(I + 1, Size, Ks, Ps, Ds, !Shft, !Mask, !Seed)
+ else
+ true
+ ).
+
+%---------------------------------------------------------------------------%
+
+rand(RP, N, RS0, RS) :-
+ RS0 = state(Seed0),
+ Size = array.size(Seed0),
+ rand_2(RP, 0, Size, 0u32, N, Seed0, Seed),
+ RS = unsafe_promise_unique(state(Seed)).
+
+:- pred rand_2(params, int, int, uint32, uint32, array(uint32), array(uint32)).
+:- mode rand_2(in, in, in, in, out, array_di, array_uo) is det.
+
+rand_2(RP, I, Size, N0, N, !Seed) :-
+ ( if I < Size then
+ Q = array.lookup(RP ^ qs, I),
+ P = array.lookup(RP ^ ps, I),
+ Shft = array.lookup(RP ^ shft, I),
+ Mask = array.lookup(RP ^ mask, I),
+ S0 = array.lookup(!.Seed, I),
+ B = ((S0 << Q) `xor` S0) >> Shft,
+ S = ((S0 /\ Mask) << P) `xor` B,
+ array.set(I, S, !Seed),
+ N1 = N0 `xor` S,
+ rand_2(RP, I + 1, Size, N1, N, !Seed)
+ else
+ N = N0
+ ).
+
+%---------------------------------------------------------------------------%
+%---------------------------------------------------------------------------%
+
+init_t3(RP, RS) :-
+ seed_t3(no, 0u32, 0u32, 0u32, RP, RS).
+
+seed_t3(MZ, A, B, C, RP, RS) :-
+ (
+ MZ = yes(Z)
+ ;
+ MZ = no,
+ Z = 0
+ ),
+ ( if params_t3(Z mod 2, Q1, Q2, Q3, P1, P2, P3) then
+ Qs = array([Q1, Q2, Q3]),
+ Ps = array([P1, P2, P3])
+ else
+ unexpected($pred, "unexpected failure")
+ ),
+ Seed = array([A, B, C]),
+ seed(Qs, Ps, Seed, RP, RS).
+
+:- pred params_t3(int, int, int, int, int, int, int).
+:- mode params_t3(in, out, out, out, out, out, out) is semidet.
+
+params_t3(0, 13, 2, 3, 12, 4, 17).
+params_t3(1, 3, 2, 13, 20, 16, 7).
+
+%---------------------------------------------------------------------------%
+
+init_t4(RP, RS) :-
+ seed_t4(no, 0u32, 0u32, 0u32, 0u32, RP, RS).
+
+seed_t4(MZ, A, B, C, D, RP, RS) :-
+ (
+ MZ = yes(Z)
+ ;
+ MZ = no,
+ Z = 58
+ ),
+ ( if params_t4(Z mod 62, P1, P2, P3, P4) then
+ Qs = array([6, 2, 13, 3]),
+ Ps = array([P1, P2, P3, P4])
+ else
+ unexpected($pred, "unexpected failure")
+ ),
+ Seed = array([A, B, C, D]),
+ seed(Qs, Ps, Seed, RP, RS).
+
+:- pred params_t4(int, int, int, int, int).
+:- mode params_t4(in, out, out, out, out) is semidet.
+
+params_t4(0, 18, 2, 7, 13).
+params_t4(1, 13, 3, 4, 9).
+params_t4(2, 24, 3, 11, 12).
+params_t4(3, 10, 4, 2, 6).
+params_t4(4, 16, 4, 2, 12).
+params_t4(5, 11, 5, 4, 3).
+params_t4(6, 17, 5, 4, 6).
+params_t4(7, 12, 5, 11, 9).
+params_t4(8, 23, 5, 11, 12).
+params_t4(9, 23, 6, 7, 8).
+params_t4(10, 14, 8, 2, 9).
+params_t4(11, 22, 8, 7, 4).
+params_t4(12, 21, 8, 11, 4).
+params_t4(13, 10, 9, 8, 2).
+params_t4(14, 22, 9, 11, 9).
+params_t4(15, 3, 10, 4, 15).
+params_t4(16, 24, 10, 7, 8).
+params_t4(17, 21, 10, 8, 4).
+params_t4(18, 12, 10, 8, 15).
+params_t4(19, 17, 10, 11, 6).
+params_t4(20, 3, 11, 4, 12).
+params_t4(21, 9, 11, 4, 13).
+params_t4(22, 9, 11, 7, 4).
+params_t4(23, 11, 12, 4, 10).
+params_t4(24, 20, 12, 7, 15).
+params_t4(25, 17, 12, 11, 11).
+params_t4(26, 21, 13, 4, 14).
+params_t4(27, 11, 14, 8, 7).
+params_t4(28, 6, 14, 8, 13).
+params_t4(29, 20, 15, 7, 13).
+params_t4(30, 12, 16, 2, 10).
+params_t4(31, 4, 16, 8, 3).
+params_t4(32, 22, 17, 4, 6).
+params_t4(33, 21, 17, 4, 13).
+params_t4(34, 20, 17, 7, 8).
+params_t4(35, 19, 17, 11, 6).
+params_t4(36, 4, 17, 11, 7).
+params_t4(37, 12, 17, 11, 15).
+params_t4(38, 15, 18, 4, 9).
+params_t4(39, 17, 18, 4, 15).
+params_t4(40, 12, 18, 7, 4).
+params_t4(41, 15, 18, 8, 11).
+params_t4(42, 6, 18, 11, 13).
+params_t4(43, 8, 19, 2, 9).
+params_t4(44, 13, 19, 4, 2).
+params_t4(45, 5, 19, 8, 3).
+params_t4(46, 6, 19, 8, 11).
+params_t4(47, 24, 19, 11, 5).
+params_t4(48, 6, 20, 2, 10).
+params_t4(49, 13, 20, 4, 10).
+params_t4(50, 24, 21, 2, 7).
+params_t4(51, 14, 21, 8, 13).
+params_t4(52, 10, 22, 8, 13).
+params_t4(53, 7, 22, 8, 14).
+params_t4(54, 15, 23, 8, 5).
+params_t4(55, 9, 23, 11, 4).
+params_t4(56, 20, 24, 4, 8).
+params_t4(57, 16, 24, 4, 14).
+params_t4(58, 20, 24, 4, 14).
+params_t4(59, 23, 24, 7, 3).
+params_t4(60, 14, 24, 8, 10).
+params_t4(61, 16, 24, 11, 12).
+
+%---------------------------------------------------------------------------%
diff --git a/library/uint32.m b/library/uint32.m
index 1c16900..11beb36 100644
--- a/library/uint32.m
+++ b/library/uint32.m
@@ -75,7 +75,7 @@
%---------------------------------------------------------------------------%
%
-% Conversion to uint64.
+% Conversion to/from uint64.
%
% cast_to_uint64(U32) = U64:
@@ -86,6 +86,14 @@
%
:- func cast_to_uint64(uint32) = uint64.
+ % cast_from_uint64(U64) = U32:
+ %
+ % Convert a uint64 to a uint32.
+ % Always succeeds, but will yield a result that is mathematically equal
+ % to I only if I is in [0, 2^32 - 1].
+ %
+:- func cast_from_uint64(uint64) = uint32.
+
%---------------------------------------------------------------------------%
%
% Change of signedness.
@@ -514,6 +522,35 @@ cast_to_uint64(_) = _ :-
%---------------------------------------------------------------------------%
+:- pragma no_determinism_warning(cast_from_uint64/1).
+
+:- pragma foreign_proc("C",
+ cast_from_uint64(U64::in) = (U32::out),
+ [will_not_call_mercury, promise_pure, thread_safe, will_not_modify_trail,
+ does_not_affect_liveness],
+"
+ U32 = (uint32_t) U64;
+").
+
+:- pragma foreign_proc("C#",
+ cast_from_uint64(U64::in) = (U32::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ U32 = (uint) U64;
+").
+
+:- pragma foreign_proc("Java",
+ cast_from_uint64(U64::in) = (U32::out),
+ [will_not_call_mercury, promise_pure, thread_safe],
+"
+ U32 = (int) U64;
+").
+
+cast_from_uint64(_) = _ :-
+ sorry($module, "uint32.cast_from_uint64/1 NYI for Erlang").
+
+%---------------------------------------------------------------------------%
+
:- pragma no_determinism_warning(cast_from_int32/1).
:- pragma foreign_proc("C",
diff --git a/tests/hard_coded/Mmakefile b/tests/hard_coded/Mmakefile
index 58395dc..ccbd972 100644
--- a/tests/hard_coded/Mmakefile
+++ b/tests/hard_coded/Mmakefile
@@ -321,6 +321,8 @@ ORDINARY_PROGS = \
rev_arith \
reverse_arith \
rnd \
+ rng1 \
+ rng2 \
rtree_test \
rtti_strings \
sectag_bits \
diff --git a/tests/hard_coded/rng1.exp b/tests/hard_coded/rng1.exp
new file mode 100644
index 0000000..228e6fd
--- /dev/null
+++ b/tests/hard_coded/rng1.exp
@@ -0,0 +1,65 @@
+marsaglia:
+1168299085
+520487819
+1761612921
+3632618539
+610669668
+2136514290
+3850311835
+2494138816
+3923280858
+1280618954
+309986706
+924303156
+2252542156
+1444019197
+2955985350
+1185139548
+3579107875
+3047601897
+1651990379
+2165617597
+
+tausworthe3:
+364603069
+528378279
+1153580643
+643237034
+3988596671
+1788716332
+626833507
+3768515118
+3526246283
+979916873
+497809124
+3522765921
+1904307014
+4035450154
+758388753
+2195520256
+1345056435
+1718236369
+823666345
+2531321601
+
+tausworthe4:
+3298080016
+1006674250
+784842863
+3826950035
+1766034713
+2314274634
+2461174380
+1680209578
+3954198082
+1441070313
+3013911521
+3001839125
+563675899
+2431136453
+632203520
+1481012674
+3251476639
+4143656215
+2141916911
+1746317775
diff --git a/tests/hard_coded/rng1.m b/tests/hard_coded/rng1.m
new file mode 100644
index 0000000..0664669
--- /dev/null
+++ b/tests/hard_coded/rng1.m
@@ -0,0 +1,42 @@
+%---------------------------------------------------------------------------%
+% vim: ts=4 sw=4 sts=4 et ft=mercury
+%---------------------------------------------------------------------------%
+
+:- module rng1.
+:- interface.
+:- import_module io.
+
+:- pred main(io::di, io::uo) is det.
+
+:- implementation.
+
+:- import_module int.
+:- import_module rng.
+:- import_module rng.marsaglia.
+:- import_module rng.tausworthe.
+
+main(!IO) :-
+ io.write_string("marsaglia:\n", !IO),
+ make_urng(marsaglia.init, RP1, RS1),
+ test(20, RP1, RS1, _, !IO),
+
+ io.write_string("\ntausworthe3:\n", !IO),
+ tausworthe.init_t3(RP2, RS2),
+ test(20, RP2, RS2, _, !IO),
+
+ io.write_string("\ntausworthe4:\n", !IO),
+ tausworthe.init_t4(RP3, RS3),
+ test(20, RP3, RS3, _, !IO).
+
+:- pred test(int, RP, RS, RS, io, io) <= urng(RP, RS).
+:- mode test(in, in, di, uo, di, uo) is det.
+
+test(Count, RP, !RS, !IO) :-
+ ( if Count > 0 then
+ urandom(RP, N, !RS),
+ io.write_uint64(N, !IO),
+ io.nl(!IO),
+ test(Count - 1, RP, !RS, !IO)
+ else
+ true
+ ).
diff --git a/tests/hard_coded/rng2.exp b/tests/hard_coded/rng2.exp
new file mode 100644
index 0000000..228e6fd
--- /dev/null
+++ b/tests/hard_coded/rng2.exp
@@ -0,0 +1,65 @@
+marsaglia:
+1168299085
+520487819
+1761612921
+3632618539
+610669668
+2136514290
+3850311835
+2494138816
+3923280858
+1280618954
+309986706
+924303156
+2252542156
+1444019197
+2955985350
+1185139548
+3579107875
+3047601897
+1651990379
+2165617597
+
+tausworthe3:
+364603069
+528378279
+1153580643
+643237034
+3988596671
+1788716332
+626833507
+3768515118
+3526246283
+979916873
+497809124
+3522765921
+1904307014
+4035450154
+758388753
+2195520256
+1345056435
+1718236369
+823666345
+2531321601
+
+tausworthe4:
+3298080016
+1006674250
+784842863
+3826950035
+1766034713
+2314274634
+2461174380
+1680209578
+3954198082
+1441070313
+3013911521
+3001839125
+563675899
+2431136453
+632203520
+1481012674
+3251476639
+4143656215
+2141916911
+1746317775
diff --git a/tests/hard_coded/rng2.m b/tests/hard_coded/rng2.m
new file mode 100644
index 0000000..c9e36da
--- /dev/null
+++ b/tests/hard_coded/rng2.m
@@ -0,0 +1,44 @@
+%---------------------------------------------------------------------------%
+% vim: ts=4 sw=4 sts=4 et ft=mercury
+%---------------------------------------------------------------------------%
+
+:- module rng2.
+:- interface.
+:- import_module io.
+
+:- pred main(io::di, io::uo) is det.
+
+:- implementation.
+
+:- import_module int.
+:- import_module rng.
+:- import_module rng.marsaglia.
+:- import_module rng.tausworthe.
+
+main(!IO) :-
+ io.write_string("marsaglia:\n", !IO),
+ RNG1 = marsaglia.init,
+ test(20, RNG1, _, !IO),
+
+ io.write_string("\ntausworthe3:\n", !IO),
+ tausworthe.init_t3(RP2, RS2),
+ RNG2 = make_shared_rng(RP2, RS2),
+ test(20, RNG2, _, !IO),
+
+ io.write_string("\ntausworthe4:\n", !IO),
+ tausworthe.init_t4(RP3, RS3),
+ RNG3 = make_shared_rng(RP3, RS3),
+ test(20, RNG3, _, !IO).
+
+:- pred test(int, RNG, RNG, io, io) <= rng(RNG).
+:- mode test(in, in, out, di, uo) is det.
+
+test(Count, !RNG, !IO) :-
+ ( if Count > 0 then
+ random(N, !RNG),
+ io.write_uint64(N, !IO),
+ io.nl(!IO),
+ test(Count - 1, !RNG, !IO)
+ else
+ true
+ ).
More information about the reviews
mailing list