diff --git a/changelog/std-random-MersenneTwisterEngine.dd b/changelog/std-random-MersenneTwisterEngine.dd new file mode 100644 index 00000000000..22cf1d6bbcb --- /dev/null +++ b/changelog/std-random-MersenneTwisterEngine.dd @@ -0,0 +1,16 @@ +`MersenneTwisterEngine` has been updated so that its template signature +matches the C++11 standard (adding two new template parameters, an extra +tempering parameter `d` and the initialization multiplier `f`). + +For anyone using the standard template instantiation `Mt19937` this will +have no noticeable effect. However, this will be a breaking change for +anyone using the `MersenneTwisterEngine` template directly. + +The internal implementation has been reworked to use Ilya Yaroshenko's +highly optimized algorithm from `mir.random`. This should have a very +noticeable positive effect for anyone who cares about generating a lot +of random numbers quickly. + +A new `Mt19937_64` template instantiation has been added, corresponding +to the standard 64-bit implementation of the algorithm (MT19937-64). +This fixes $(LINK https://issues.dlang.org/show_bug.cgi?id=10900). diff --git a/std/random.d b/std/random.d index 47d586a04a6..f3b1b455ce2 100644 --- a/std/random.d +++ b/std/random.d @@ -44,6 +44,7 @@ License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0). Authors: $(HTTP erdani.org, Andrei Alexandrescu) Masahiro Nakagawa (Xorshift random generator) $(HTTP braingam.es, Joseph Rushton Wakeling) (Algorithm D for random sampling) + Ilya Yaroshenko (Mersenne Twister implementation, adapted from $(HTTPS github.com/libmir/mir-random, mir.random)) Credits: The entire random number library architecture is derived from the excellent $(HTTP open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2461.pdf, C++0X) random number facility proposed by Jens Maurer and contributed to by @@ -527,9 +528,9 @@ alias MinstdRand = LinearCongruentialEngine!(uint, 48271, 0, 2147483647); The $(LUCKY Mersenne Twister) generator. */ struct MersenneTwisterEngine(UIntType, size_t w, size_t n, size_t m, size_t r, - UIntType a, size_t u, size_t s, + UIntType a, size_t u, UIntType d, size_t s, UIntType b, size_t t, - UIntType c, size_t l) + UIntType c, size_t l, UIntType f) if (isUnsigned!UIntType) { static assert(0 < w && w <= UIntType.sizeof * 8); @@ -537,6 +538,7 @@ struct MersenneTwisterEngine(UIntType, size_t w, size_t n, size_t m, size_t r, static assert(0 <= r && 0 <= u && 0 <= s && 0 <= t && 0 <= l); static assert(r <= w && u <= w && s <= w && t <= w && l <= w); static assert(0 <= a && 0 <= b && 0 <= c); + static assert(n <= sizediff_t.max); ///Mark this as a Rng enum bool isUniformRandom = true; @@ -549,21 +551,79 @@ Parameters for the generator. enum size_t shiftSize = m; /// ditto enum size_t maskBits = r; /// ditto enum UIntType xorMask = a; /// ditto - enum UIntType temperingU = u; /// ditto + enum size_t temperingU = u; /// ditto + enum UIntType temperingD = d; /// ditto enum size_t temperingS = s; /// ditto enum UIntType temperingB = b; /// ditto enum size_t temperingT = t; /// ditto enum UIntType temperingC = c; /// ditto enum size_t temperingL = l; /// ditto + enum UIntType initializationMultiplier = f; /// ditto /// Smallest generated value (0). enum UIntType min = 0; /// Largest generated value. enum UIntType max = UIntType.max >> (UIntType.sizeof * 8u - w); - static assert(a <= max && b <= max && c <= max); + // note, `max` also serves as a bitmask for the lowest `w` bits + static assert(a <= max && b <= max && c <= max && f <= max); + /// The default seed value. enum UIntType defaultSeed = 5489u; + // Bitmasks used in the 'twist' part of the algorithm + private enum UIntType lowerMask = (cast(UIntType) 1u << r) - 1, + upperMask = (~lowerMask) & this.max; + + /* + Collection of all state variables + used by the generator + */ + private struct State + { + /* + State array of the generator. This + is iterated through backwards (from + last element to first), providing a + few extra compiler optimizations by + comparison to the forward iteration + used in most implementations. + */ + UIntType[n] data; + + /* + Cached copy of most recently updated + element of `data` state array, ready + to be tempered to generate next + `front` value + */ + UIntType z; + + /* + Most recently generated random variate + */ + UIntType front; + + /* + Index of the entry in the `data` + state array that will be twisted + in the next `popFront()` call + */ + size_t index; + } + + /* + State variables used by the generator; + initialized to values equivalent to + explicitly seeding the generator with + `defaultSeed` + */ + private State state = defaultState(); + /* NOTE: the above is a workaround to ensure + backwards compatibility with the original + implementation, which permitted implicit + construction. With `@disable this();` + it would not be necessary. */ + /** Constructs a MersenneTwisterEngine object. */ @@ -572,36 +632,60 @@ Parameters for the generator. seed(value); } + /** + Generates the default initial state for a Mersenne + Twister; equivalent to the internal state obtained + by calling `seed(defaultSeed)` + */ + private static State defaultState() @safe pure nothrow @nogc + { + if (!__ctfe) assert(false); + State mtState; + seedImpl(defaultSeed, mtState); + return mtState; + } + /** Seeds a MersenneTwisterEngine object. Note: - This seed function gives 2^32 starting points. To allow the RNG to be started in any one of its - internal states use the seed overload taking an InputRange. + This seed function gives 2^w starting points (the lowest w bits of + the value provided will be used). To allow the RNG to be started + in any one of its internal states use the seed overload taking an + InputRange. */ void seed()(UIntType value = defaultSeed) @safe pure nothrow @nogc { - static if (w == UIntType.sizeof * 8) - { - mt[0] = value; - } - else + this.seedImpl(value, this.state); + } + + /** + Implementation of the seeding mechanism, which + can be used with an arbitrary `State` instance + */ + private static void seedImpl(UIntType value, ref State mtState) + { + mtState.data[$ - 1] = value; + static if (this.max != UIntType.max) { - static assert(max + 1 > 0); - mt[0] = value % (max + 1); + mtState.data[$ - 1] &= this.max; } - for (mti = 1; mti < n; ++mti) + + foreach_reverse (size_t i, ref e; mtState.data[0 .. $ - 1]) { - mt[mti] = - cast(UIntType) - (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> (w - 2))) + mti); - /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ - /* In the previous versions, MSBs of the seed affect */ - /* only MSBs of the array mt[]. */ - /* 2002/01/09 modified by Makoto Matsumoto */ - //mt[mti] &= ResultType.max; - /* for >32 bit machines */ + e = f * (mtState.data[i + 1] ^ (mtState.data[i + 1] >> (w - 2))) + cast(UIntType)(n - (i + 1)); + static if (this.max != UIntType.max) + { + e &= this.max; + } } - popFront(); + + mtState.index = n - 1; + + /* double popFront() to guarantee both `mtState.z` + and `mtState.front` are derived from the newly + set values in `mtState.data` */ + MersenneTwisterEngine.popFrontImpl(mtState); + MersenneTwisterEngine.popFrontImpl(mtState); } /** @@ -612,14 +696,26 @@ Parameters for the generator. The number of elements required is the 'n' template parameter of the MersenneTwisterEngine struct. */ void seed(T)(T range) if (isInputRange!T && is(Unqual!(ElementType!T) == UIntType)) + { + this.seedImpl(range, this.state); + } + + /** + Implementation of the range-based seeding mechanism, + which can be used with an arbitrary `State` instance + */ + private static void seedImpl(T)(T range, ref State mtState) + if (isInputRange!T && is(Unqual!(ElementType!T) == UIntType)) { size_t j; for (j = 0; j < n && !range.empty; ++j, range.popFront()) { - mt[j] = range.front; + sizediff_t idx = n - j - 1; + mtState.data[idx] = range.front; } - mti = n; + mtState.index = n - 1; + if (range.empty && j < n) { import core.internal.string : UnsignedStringBuf, unsignedToTempString; @@ -630,17 +726,11 @@ Parameters for the generator. throw new Exception(s); } - popFront(); - } - - /// - @safe unittest - { - import std.algorithm.iteration : map; - import std.range : repeat; - - Mt19937 gen; - gen.seed(map!((a) => unpredictableSeed)(repeat(0))); + /* double popFront() to guarantee both `mtState.z` + and `mtState.front` are derived from the newly + set values in `mtState.data` */ + MersenneTwisterEngine.popFrontImpl(mtState); + MersenneTwisterEngine.popFrontImpl(mtState); } /** @@ -648,59 +738,71 @@ Parameters for the generator. */ void popFront() @safe pure nothrow @nogc { - if (mti == size_t.max) seed(); - enum UIntType - upperMask = ~((cast(UIntType) 1u << - (UIntType.sizeof * 8 - (w - r))) - 1), - lowerMask = (cast(UIntType) 1u << r) - 1; - static immutable UIntType[2] mag01 = [0x0UL, a]; + this.popFrontImpl(this.state); + } - ulong y = void; + /* + Internal implementation of `popFront()`, which + can be used with an arbitrary `State` instance + */ + private static void popFrontImpl(ref State mtState) + { + /* This function blends two nominally independent + processes: (i) calculation of the next random + variate `mtState.front` from the cached previous + `data` entry `z`, and (ii) updating the value + of `data[index]` and `mtState.z` and advancing + the `index` value to the next in sequence. - if (mti >= n) - { - /* generate N words at one time */ + By interweaving the steps involved in these + procedures, rather than performing each of + them separately in sequence, the variables + are kept 'hot' in CPU registers, allowing + for significantly faster performance. */ + sizediff_t index = mtState.index; + sizediff_t next = index - 1; + if (next < 0) + next = n - 1; + auto z = mtState.z; + sizediff_t conj = index - m; + if (conj < 0) + conj = index - m + n; - int kk = 0; - const limit1 = n - m; - for (; kk < limit1; ++kk) - { - y = (mt[kk] & upperMask)|(mt[kk + 1] & lowerMask); - mt[kk] = cast(UIntType) (mt[kk + m] ^ (y >> 1) - ^ mag01[cast(UIntType) y & 0x1U]); - } - const limit2 = n - 1; - for (; kk < limit2; ++kk) - { - y = (mt[kk] & upperMask)|(mt[kk + 1] & lowerMask); - mt[kk] = cast(UIntType) (mt[kk + (m -n)] ^ (y >> 1) - ^ mag01[cast(UIntType) y & 0x1U]); - } - y = (mt[n -1] & upperMask)|(mt[0] & lowerMask); - mt[n - 1] = cast(UIntType) (mt[m - 1] ^ (y >> 1) - ^ mag01[cast(UIntType) y & 0x1U]); - - mti = 0; + static if (d == UIntType.max) + { + z ^= (z >> u); + } + else + { + z ^= (z >> u) & d; } - y = mt[mti++]; - - /* Tempering */ - y ^= (y >> temperingU); - y ^= (y << temperingS) & temperingB; - y ^= (y << temperingT) & temperingC; - y ^= (y >> temperingL); + auto q = mtState.data[index] & upperMask; + auto p = mtState.data[next] & lowerMask; + z ^= (z << s) & b; + auto y = q | p; + auto x = y >> 1; + z ^= (z << t) & c; + if (y & 1) + x ^= a; + auto e = mtState.data[conj] ^ x; + z ^= (z >> l); + mtState.z = mtState.data[index] = e; + mtState.index = next; - _y = cast(UIntType) y; + /* technically we should take the lowest `w` + bits here, but if the tempering bitmasks + `b` and `c` are set correctly, this should + be unnecessary */ + mtState.front = z; } /** Returns the current random value. */ - @property UIntType front() @safe pure nothrow @nogc + @property UIntType front() @safe const pure nothrow @nogc { - if (mti == size_t.max) seed(); - return _y; + return this.state.front; } /// @@ -713,10 +815,6 @@ Parameters for the generator. Always $(D false). */ enum bool empty = false; - - private UIntType[n] mt; - private size_t mti = size_t.max; /* means mt is not initialized */ - UIntType _y = UIntType.max; } /** @@ -728,9 +826,9 @@ generation unless memory is severely restricted, in which case a $(D LinearCongruentialEngine) would be the generator of choice. */ alias Mt19937 = MersenneTwisterEngine!(uint, 32, 624, 397, 31, - 0x9908b0df, 11, 7, + 0x9908b0df, 11, 0xffffffff, 7, 0x9d2c5680, 15, - 0xefc60000, 18); + 0xefc60000, 18, 1812433253); /// @safe unittest @@ -753,8 +851,60 @@ alias Mt19937 = MersenneTwisterEngine!(uint, 32, 624, 397, 31, static assert(isSeedable!(Mt19937, uint)); static assert(isSeedable!(Mt19937, typeof(map!((a) => unpredictableSeed)(repeat(0))))); Mt19937 gen; + assert(gen.front == 3499211612); popFrontN(gen, 9999); assert(gen.front == 4123659995); + try { gen.seed(iota(624u)); } catch (Exception) { assert(false); } + assert(gen.front == 3708921088u); + popFrontN(gen, 9999); + assert(gen.front == 165737292u); +} + +/** +A $(D MersenneTwisterEngine) instantiated with the parameters of the +original engine $(HTTP en.wikipedia.org/wiki/Mersenne_Twister, +MT19937-64), generating uniformly-distributed 64-bit numbers with a +period of 2 to the power of 19937. +*/ +alias Mt19937_64 = MersenneTwisterEngine!(ulong, 64, 312, 156, 31, + 0xb5026f5aa96619e9, 29, 0x5555555555555555, 17, + 0x71d67fffeda60000, 37, + 0xfff7eee000000000, 43, 6364136223846793005); + +/// +@safe unittest +{ + // Seed with a constant + auto gen = Mt19937_64(12345); + auto n = gen.front; // same for each run + // Seed with an unpredictable value + gen.seed(unpredictableSeed); + n = gen.front; // different across runs +} + +@safe nothrow unittest +{ + import std.algorithm; + import std.range; + static assert(isUniformRNG!Mt19937_64); + static assert(isUniformRNG!(Mt19937_64, ulong)); + static assert(isSeedable!Mt19937_64); + static assert(isSeedable!(Mt19937_64, ulong)); + // FIXME: this test demonstrates viably that Mt19937_64 + // is seedable with an infinite range of `ulong` values + // but it's a poor example of how to actually seed the + // generator, since it can't cover the full range of + // possible seed values. Ideally we need a 64-bit + // unpredictable seed to complement the 32-bit one! + static assert(isSeedable!(Mt19937_64, typeof(map!((a) => (cast(ulong)unpredictableSeed))(repeat(0))))); + Mt19937_64 gen; + assert(gen.front == 14514284786278117030uL); + popFrontN(gen, 9999); + assert(gen.front == 9981545732273789042uL); + try { gen.seed(iota(312uL)); } catch (Exception) { assert(false); } + assert(gen.front == 14660652410669508483uL); + popFrontN(gen, 9999); + assert(gen.front == 15956361063660440239uL); } @safe unittest @@ -792,7 +942,7 @@ alias Mt19937 = MersenneTwisterEngine!(uint, 32, 624, 397, 31, { import std.range; // Check .save works - foreach (Type; std.meta.AliasSeq!(Mt19937)) + foreach (Type; std.meta.AliasSeq!(Mt19937, Mt19937_64)) { auto gen1 = Type(unpredictableSeed); auto gen2 = gen1.save; @@ -806,12 +956,24 @@ alias Mt19937 = MersenneTwisterEngine!(uint, 32, 624, 397, 31, @safe pure nothrow unittest //11690 { alias MT(UIntType, uint w) = MersenneTwisterEngine!(UIntType, w, 624, 397, 31, - 0x9908b0df, 11, 7, + 0x9908b0df, 11, 0xffffffff, 7, 0x9d2c5680, 15, - 0xefc60000, 18); + 0xefc60000, 18, 1812433253); - foreach (R; std.meta.AliasSeq!(MT!(uint, 32), MT!(ulong, 32), MT!(ulong, 48), MT!(ulong, 64))) + ulong[] expectedFirstValue = [3499211612uL, 3499211612uL, + 171143175841277uL, 1145028863177033374uL]; + + ulong[] expected10kValue = [4123659995uL, 4123659995uL, + 51991688252792uL, 3031481165133029945uL]; + + foreach (i, R; std.meta.AliasSeq!(MT!(uint, 32), MT!(ulong, 32), MT!(ulong, 48), MT!(ulong, 64))) + { auto a = R(); + a.seed(a.defaultSeed); // checks that some alternative paths in `seed` are utilized + assert(a.front == expectedFirstValue[i]); + a.popFrontN(9999); + assert(a.front == expected10kValue[i]); + } }